Многомерная нормальная плотность в Python?

Многомерная нормаль теперь доступна на 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]) 

Я просто сделал один для своих целей, поэтому, хотя я бы поделился. Он построен с использованием «полномочий» 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) 

В общем случае диагональной ковариационной матрицы многомерный PDF можно получить простым умножением одномерных значений PDF, возвращаемых экземпляром scipy.stats.norm . Если вам нужен общий случай, вам, вероятно, придется самому закодировать его (что не должно быть сложно).

Если все еще необходимо, моя реализация будет

 import numpy as np def pdf_multivariate_gauss(x, mu, cov): ''' Caculate the multivariate normal density (pdf) Keyword arguments: x = numpy array of a "dx 1" sample vector mu = numpy array of a "dx 1" mean vector cov = "numpy array of adxd" 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() 

В случае, если я буду вносить будущие изменения, код здесь находится на GitHub

Я знаю несколько пакетов 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

По сути, если вам нужна быстрая оценка, перепишите его для вашего случая использования.

Я использую следующий код, который вычисляет значение 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 () в возвращаемое значение

Плотность можно вычислить довольно простым способом, используя функции numpy и формулу на этой странице: http://en.wikipedia.org/wiki/Multivariate_normal_distribution . Вы также можете использовать функцию правдоподобия (логарифмическая вероятность), которая менее вероятна для нижнего уровня для больших измерений и немного более простая для вычисления. Оба просто включают в себя возможность вычислить детерминант и инвертировать матрицу.

С другой стороны, CDF является совершенно другим животным …