时间戳

首页 > 解决方案 > scipy.stats.multivariate_normal.pdf 与使用 numpy 编写的相同函数有何不同？

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

问题描述

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

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

标签： pythonpython-3.xnumpyscipyprobability

解决方案

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

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

推荐阅读