Densità normale multivariata in Python?

Question

Esiste un pacchetto Python che consente il calcolo efficiente del pdf normale multivariato?

Non sembra essere incluso in Numpy/Scipy, e sorprendentemente una ricerca su Google non ha rivelato nulla di utile.

Answer 1

Il multivariata normale è ora disponibile su 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

Nel caso comune di una matrice di covarianza diagonale, il PDF multivariato può essere ottenuto semplicemente moltiplicando i valori PDF univariati restituiti da un'istanza scipy.stats.norm. Se hai bisogno del caso generale, probabilmente dovrai codificarlo da solo (il che non dovrebbe essere difficile).

Answer 3

Conosco diversi pacchetti Python che lo utilizzano internamente, con diverse generalità e per diversi usi, ma non so se qualcuno di essi è destinato agli utenti.

statsmodels, ad esempio, ha la seguente funzione nascosta e di classe, ma non è utilizzato da 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

In sostanza, se avete bisogno di valutazione veloce, riscriverlo per la vostra caso d'uso.

Answer 4

La densità può essere calcolata in modo piuttosto semplice utilizzando le funzioni numpy e la formula in questa pagina: http://en.wikipedia.org/wiki/Multivariate_normal_distribution. Potresti anche voler utilizzare la funzione di verosimiglianza (probabilità di log), che ha meno probabilità di underflow per le grandi dimensioni ed è un po 'più semplice da calcolare. Entrambi implicano solo la capacità di calcolare il determinante e l'inverso di una matrice.

Il CDF, d'altra parte, è un animale completamente diverso ...

Answer 5

ho appena fatto uno per i miei scopi così ho pensato che vorrei condividere. È costruito usando "i poteri" di numpy, sulla formula del caso non degenerato da http://en.wikipedia.org/wiki/Multivariate_normal_distribution e convalida anche l'input.

Ecco il codice con un campione di corsa

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

uso il seguente codice che calcola il valore logpdf, che è preferibile per dimensioni maggiori. Funziona anche con matrici 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) 

Codice è da pyParticleEst, se si desidera che il valore pdf invece del logpdf basta prendere math.exp() sul valore restituito

Answer 7

Se ancora bisogno, la mia implementazione sarebbe

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() 

In caso di modifiche future, il codice è here on GitHub