我只想加速我用numpy编写的数值算法。一个关键的部分是计算对数似然函数(两个截断正态CDF之间的差)。我的函数非常慢(每个循环31.9ms),每次迭代都需要运行2000次。在
我试过用西皮的”标准cdf“功能”而不是“ecfc”。但速度较慢。我也尝试了Numba软件包中的“@jit”。但它也比原始代码慢。在
我想也许我需要用Cython。但我对C几乎一无所知。我试着从Cython for numpy users网页上学习Cython,但这对我来说真的很难。在
有人能帮我用Cython重写代码吗?或者建议我怎样写得更快?在
import numpy as np
from scipy.special import erfc
# The bloody function for calculating the difference between two truncated normal CDFs
def my_loglikelihood2(x,b,c,z):
log_likelihood=np.zeros(np.shape(z)[0])
log_likelihood[x==1]=np.log(0.5*erfc(-(c[1]-np.dot(z[x==1,:],b)) / np.sqrt(2.)) - 0.5*erfc(-(c[0]-np.dot(z[x==1],b)) / np.sqrt(2.)))
log_likelihood[x==2]=np.log(0.5*erfc(-(c[2]-np.dot(z[x==2,:],b)) / np.sqrt(2.)) - 0.5*erfc(-(c[1]-np.dot(z[x==2],b)) / np.sqrt(2.)))
log_likelihood[x==3]=np.log(0.5*erfc(-(c[3]-np.dot(z[x==3,:],b)) / np.sqrt(2.)) - 0.5*erfc(-(c[2]-np.dot(z[x==3],b)) / np.sqrt(2.)))
return log_likelihood
# generate random values
x=np.random.randint(low=1, high=4, size=50000)
b=np.random.normal(0,1,70)
c=np.array([-999,-1,1,999],dtype='f')
z=np.random.multivariate_normal(np.zeros(70), np.eye(70), 50000)
%timeit my_loglikelihood2(x,b,c,z)
# 10 loops, best of 3: 31.9 ms per loop :(
更新1-根据@jackvdp的建议。已接种4.5倍。但我仍在寻找更快的代码:
^{pr2}$更新2-根据@DSM的建议。更换美国运输部(z,b)通过zdotb=z.dot(b)。提高1.5ms
def my_loglikelihood2(x,b,low_c,up_c,z):
up_c=up_cutoff(x,c)
low_c=low_cutoff(x,c)
zdotb = z.dot(b)
return np.log(0.5*erfc(-(up_c-zdotb) / np.sqrt(2.)) - 0.5*erfc(-(low_c-zdotb) / np.sqrt(2.)))
%timeit my_loglikelihood2(x,b,low_c,up_c,z)
100 loops, best of 3: 5.02 ms per loop
如果您的代码由于Python中的循环而变得很慢,那么将它移植到Cython可以看到很大的改进。但是您的示例只调用现有的numpy/scipy函数六次。在
它主要是对
np.log, erfc, np.dot, np.sqrt
的调用。我不确定erfc
,但其他人已经使用编译代码。Cython不碰那些。在我们可以检查
erfc
。在但最好的办法是用更大的数组调用此代码。在
相关问题 更多 >
编程相关推荐