Gęstość normalna wielowymiarowa w języku Python?

Question

Czy istnieje pakiet python, który pozwala na wydajne obliczanie wielowariantowego normalnego pliku pdf?

Wygląda na to, że nie ma go w Numpy/Scipy i nieoczekiwanie wyszukiwanie w Google nie przyniosło żadnych pożytecznych efektów.

Answer 1

wieloczynnikowej normalnym jest teraz dostępny na 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

W typowym przypadku macierzy kowariancji po przekątnej wielowymiarowy plik PDF można uzyskać, mnożąc jednokierunkowe wartości PDF zwracane przez instancję scipy.stats.norm. Jeśli potrzebujesz ogólnego przypadku, prawdopodobnie będziesz musiał sam to zakodować (co nie powinno być trudne).

Answer 3

Znam kilka pakietów Pythona, które używają go wewnętrznie, z różną ogólnością i dla różnych zastosowań, ale nie wiem, czy któryś z nich jest przeznaczony dla użytkowników.

statsmodels, na przykład, ma następującą ukrytą funkcję i klasę, ale nie jest używany przez 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

Zasadniczo, jeśli potrzebujesz szybkiej oceny, przerobić go na swój przypadek użycia.

Answer 4

Gęstość można obliczyć w prosty sposób, używając funkcji numpy i wzoru na tej stronie: http://en.wikipedia.org/wiki/Multivariate_normal_distribution. Możesz również skorzystać z funkcji prawdopodobieństwa (prawdopodobieństwo zalogowania), która jest mniej prawdopodobna w przypadku dużych rozmiarów i jest nieco prostsza do obliczenia. Oba dotyczą jedynie możliwości obliczenia wyznacznika i odwrotności macierzy.

CDF, z drugiej strony, to zupełnie inna zwierząt ...

Answer 5

Zrobiłem jeden dla moich celów, więc jakbym podzielają. Jest zbudowany przy użyciu "uprawnień" numpy, na wzorze nie-zdegenerowanego przypadku od http://en.wikipedia.org/wiki/Multivariate_normal_distribution i aso sprawdza dane wejściowe.

Oto kod wraz z przebiegu próbki

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

użyć następujące kod, który oblicza wartość logpdf, co jest korzystne w przypadku większych wymiarów. Działa również dla macierzy 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) 

Code jest od pyParticleEst, jeśli chcesz wartości pdf zamiast logpdf prostu wziąć math.exp() na zwróconej wartości

Answer 7

Jeśli nadal potrzebne, moja realizacja byłaby

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

W przypadku przyszłych zmian, kod: here on GitHub