scipy.stats.multivariate_normal.pdf 与使用 numpy 编写的相同函数有何不同？

Question

我需要在脚本中使用多元正态分布。我注意到我的版本给出了与 scipy 方法不同的答案。我真的不明白为什么...

这是我的函数：

def gauss(x, mu, sigma):
    assert np.linalg.det(sigma)!=0, "determinant of sigma is 0"
    y = np.exp((-1/2)*(x-mu).T.dot(np.linalg.inv(sigma)).dot(x-mu))/np.sqrt(
      np.power(2*np.pi, len(x))*np.linalg.det(sigma)
    )
    return y

这是结果的比较：

from scipy.stats import multivariate_normal
import numpy as np

x = np.array([-0.54849176, 6.39530657])
mu = np.array([15,20])
sigma = np.array([
  [2,3],
  [4,10]
])

print(gauss(x, mu, sigma))
# output is 1.8781656851138248e-37

print(multivariate_normal.pdf(x, mu, sigma))
# output is 2.698549423643947e-61

有人注意到了吗？我的功能错了吗？任何帮助将不胜感激！

Answer 1

您用作示例的特定输入可能会有点误导，因为值太低以至于数字问题很容易导致您看到的差异。但是，即使使用密度更大的示例，您仍然会遇到问题：

In [95]: x = np.array([15.00054849176, 20.0009530657]) 
    ...: mu = np.array([15, 20]) 
    ...: sigma = np.array([ 
    ...:   [2, 3], 
    ...:   [4, 10] 
    ...: ]) 
    ...:                                                                                        

In [96]: print(gauss(x, mu, sigma)) 
    ...: print(multivariate_normal.pdf(x, mu, sigma)) 
    ...:                                                                                        
0.05626976565965294
0.07957746514880353

也许有趣的是，这种差异是 np.sqrt(2) 的一个因素，取决于数值问题，但这有点转移注意力：事实证明，这种差异只是由于您的协方差矩阵不是协方差矩阵：虽然它是半正定的，但它 不是对称的 。使用有效的输入，这两种方法确实会一致（直到数值问题）：

In [99]: x = np.array([15.00054849176, 20.0009530657]) 
    ...: mu = np.array([15, 20]) 
    ...: sigma = np.array([ 
    ...:   [2, 3], 
    ...:   [3, 10] 
    ...: ]) 
    ...:                                                                                        

In [100]: print(gauss(x, mu, sigma)) 
     ...: print(multivariate_normal.pdf(x, mu, sigma)) 
     ...:                                                                                       
0.047987017204594515
0.04798701720459451

或者，使用您的原始输入：

In [111]: x = np.array([-0.54849176, 6.39530657]) 
     ...: mu = np.array([15, 20]) 
     ...: sigma = np.array([ 
     ...:   [2, 3], 
     ...:   [3, 10] 
     ...: ]) 
     ...:                                                                                       

In [112]: print(gauss(x, mu, sigma)) 
     ...: print(multivariate_normal.pdf(x, mu, sigma)) 
     ...:                                                                                       
5.060725651214228e-32
5.060725651214157e-32

scipy.stats.multivariate_normal.pdf 与使用 numpy 编写的相同函数有何不同？

How is scipy.stats.multivariate_normal.pdf different from the same function written using numpy?

python

numpy

probability

scipy

python-3.x