Python中的多元正态密度？

Question

是否有任何允许高效计算多元正态pdf的python包？

它似乎并未包含在Numpy/Scipy中，令人惊讶的是谷歌搜索没有发现任何有用的东西。

Answer 1

多元正态分布现已有SciPy 0.14.0.dev-16fc0af：

from scipy.stats import multivariate_normal 
var = multivariate_normal(mean=[0,0], cov=[[1,0],[0,1]]) 
var.pdf([1,0]) 

Answer 2

在对角协方差矩阵的常见情况下，多元PDF可以通过简单乘以scipy.stats.norm实例返回的单变量PDF值来获得。如果你需要一般情况下，你可能必须自己编码（这不应该很难）。

Answer 3

我知道几个python包内部使用它，具有不同的通用性和不同的用途，但我不知道它们中的任何一个是否打算供用户使用。

statsmodels，例如，具有以下隐藏函数和类，但它不使用statsmodels：

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/miscmodels/try_mlecov.py#L36

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/sandbox/distributions/mv_normal.py#L777

从本质上讲，如果你需要快速评估，把它改写为你用例。

Answer 4

使用numpy函数和本页上的公式可以非常简单地计算密度：http://en.wikipedia.org/wiki/Multivariate_normal_distribution。您可能还需要使用似然函数（对数概率），对于较大的维度来说下溢的可能性较小，并且计算起来更简单一些。两者都只涉及能够计算矩阵的行列式和逆矩阵。

的CDF，在另一方面，是一个完全不同的动物......

Answer 5

我只是做了一个我的目的，所以我虽然我会分享。它使用numpy的“权力”构建，在http://en.wikipedia.org/wiki/Multivariate_normal_distribution的非简并案例的公式中，它也验证输入。

这里是一个样品运行

from numpy import * 
import math 
# covariance matrix 
sigma = matrix([[2.3, 0, 0, 0], 
      [0, 1.5, 0, 0], 
      [0, 0, 1.7, 0], 
      [0, 0, 0, 2] 
      ]) 
# mean vector 
mu = array([2,3,8,10]) 

# input 
x = array([2.1,3.5,8, 9.5]) 

def norm_pdf_multivariate(x, mu, sigma): 
    size = len(x) 
    if size == len(mu) and (size, size) == sigma.shape: 
     det = linalg.det(sigma) 
     if det == 0: 
      raise NameError("The covariance matrix can't be singular") 

     norm_const = 1.0/ (math.pow((2*pi),float(size)/2) * math.pow(det,1.0/2)) 
     x_mu = matrix(x - mu) 
     inv = sigma.I   
     result = math.pow(math.e, -0.5 * (x_mu * inv * x_mu.T)) 
     return norm_const * result 
    else: 
     raise NameError("The dimensions of the input don't match") 

print norm_pdf_multivariate(x, mu, sigma) 

Answer 6

我使用下面的代码，其计算logpdf值，这是优选为更大的尺寸沿着代码。它也适用于scipy.sparse矩阵。

import numpy as np 
import math 
import scipy.sparse as sp 
import scipy.sparse.linalg as spln 

def lognormpdf(x,mu,S): 
    """ Calculate gaussian probability density of x, when x ~ N(mu,sigma) """ 
    nx = len(S) 
    norm_coeff = nx*math.log(2*math.pi)+np.linalg.slogdet(S)[1] 

    err = x-mu 
    if (sp.issparse(S)): 
     numerator = spln.spsolve(S, err).T.dot(err) 
    else: 
     numerator = np.linalg.solve(S, err).T.dot(err) 

    return -0.5*(norm_coeff+numerator) 

代码为pyParticleEst，如果你想在PDF值而不是logpdf只取math.exp（）的返回值

Answer 7

如果还需要，我的实现将是

import numpy as np 

def pdf_multivariate_gauss(x, mu, cov): 
    ''' 
    Caculate the multivariate normal density (pdf) 

    Keyword arguments: 
     x = numpy array of a "d x 1" sample vector 
     mu = numpy array of a "d x 1" mean vector 
     cov = "numpy array of a d x d" covariance matrix 
    ''' 
    assert(mu.shape[0] > mu.shape[1]), 'mu must be a row vector' 
    assert(x.shape[0] > x.shape[1]), 'x must be a row vector' 
    assert(cov.shape[0] == cov.shape[1]), 'covariance matrix must be square' 
    assert(mu.shape[0] == cov.shape[0]), 'cov_mat and mu_vec must have the same dimensions' 
    assert(mu.shape[0] == x.shape[0]), 'mu and x must have the same dimensions' 
    part1 = 1/(((2* np.pi)**(len(mu)/2)) * (np.linalg.det(cov)**(1/2))) 
    part2 = (-1/2) * ((x-mu).T.dot(np.linalg.inv(cov))).dot((x-mu)) 
    return float(part1 * np.exp(part2)) 

def test_gauss_pdf(): 
    x = np.array([[0],[0]]) 
    mu = np.array([[0],[0]]) 
    cov = np.eye(2) 

    print(pdf_multivariate_gauss(x, mu, cov)) 

    # prints 0.15915494309189535 

if __name__ == '__main__': 
    test_gauss_pdf() 

如果我做将来的更改，代码是here on GitHub