我在python中有以下函数,我不知道如何用矢量化的形式来表示。
对我来说,cov是一个numpy形状数组(2,2),mu是具有形状(2,)的平均向量,xtp是形状(~50000,2)。
我知道scipy提供scipy.stats.multivariate_normal.pdf,但我正在努力学习如何编写有效的矢量化代码。拜托

def mvnpdf(xtp, mu, cov):
    temp = np.zeros(xtp.shape[0])
    i = 0
    length = xtp.shape[0]
    const = 1 / ( ((2* np.pi)**(len(mu)/2)) * (np.linalg.det(cov)**(1/2)) )
    inv = np.linalg.inv(cov)
    while i < length:
        x = xtp[i]-mu
        exponent = (-1/2) * (x.dot(inv).dot(x))
        temp[i] =  (const * np.exp(exponent))
        i+=1
    return temp

最佳答案

矢量化唯一棘手的部分是double.dot。让我们隔离一下:

x = xtp - mu  # move this out of the loop
ddot = [i.dot(inv).dot(i) for i in x]
temp = const * np.exp(-0.5 * ddot)

把它放到你的代码中,看看它是否会产生同样的结果。
有几种“矢量化”adot的方法。我喜欢先尝试的是einsum。在我的测试中,这相当于:
ddot = np.einsum('ij,jk,ik->i',x,inv,x)

我建议你试试看它是否有效,并加快速度。在一个交互式shell中使用更小的数组(而不是大约50000个数组)来处理这些计算。
我在测试
In [225]: x
Out[225]:
array([[  0.,   2.],
       [  1.,   3.],
       ...
       [  7.,   9.],
       [  8.,  10.],
       [  9.,  11.]])
In [226]: inv
Out[226]:
array([[ 1.,  0.],
       [ 0.,  1.]])

既然这是一个学习练习,我就把细节留给你。
使用(2,2)时,计算速度可能比使用covdet函数更快。但迭代是时间消耗者。

关于python - 如何以向量化形式编写此numpy代码?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/32533392/

10-09 05:35
查看更多