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

【问题标题】：How is scipy.stats.multivariate_normal.pdf different from the same function written using numpy?scipy.stats.multivariate_normal.pdf 与使用 numpy 编写的相同函数有何不同？
【发布时间】：2020-03-16 20:54:03
【问题描述】：

我需要在脚本中使用多元正态分布。我注意到我的版本给出了与 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

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

【问题讨论】：

当x 和mu 是一维数组时，(x-mu).T 不会像您认为的那样做。如果您想安全，请使用reshape(1, -1)。
另外，不是断言数组，而是转换为数组而不需要额外的副本。这实际上是相当标准的。
@MadPhysicist 感谢 cmets。我试过reshape(1, -1)，但似乎没有帮助
我会玩它，当我使用桌面时会通知你
密度非常低，因此先验地，人们可能会猜测您所看到的只是数字问题：以不同的方式设置一对括号就足以导致差异。但是，这里发生了其他事情，因为您也有 x = np.array([15.054849176, 20.39530657]) 的差异@

标签： python python-3.x numpy scipy probability

【解决方案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)达到数值问题的因素，但这是一个红鲱鱼的一点：事实证明，差异是简单地由你的协方差矩阵而不是协方差矩阵：虽然它是正半定的，但它是不对称 em>。使用有效输入，两种方法确实同意（最多为数字问题）：

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

【讨论】：

正确！我应该在我的问题中提到我尝试使用X作为获得最高概率，我仍然看到了差异。而且我确实尝试了身份矩阵作为协方差，这使得两个输出相同。但是一旦我改变了它，差异就回来了！ span>
伟大;如果答案很有帮助，您可以mark it as accepted;除了为我提供人工互联网积分外，这有助于用户确定哪些问题仍需要答案。 span>