python - 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