scratch - python machine learning packt pdf
¿Densidad normal multivariada en Python? (7)
¿Hay algún paquete de Python que permita el cálculo eficiente del pdf normal multivariante?
No parece estar incluido en Numpy / Scipy, y sorprendentemente una búsqueda en Google no apareció nada útil.
Acabo de hacer uno para mis propósitos, así que pensé que lo compartiría. Está construido usando "los poderes" de numpy, en la fórmula del caso no degenerado de http://en.wikipedia.org/wiki/Multivariate_normal_distribution y así valida la entrada.
Aquí está el código junto con una ejecución de muestra
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)
Conozco varios paquetes de Python que lo usan internamente, con diferente generalidad y para diferentes usos, pero no sé si alguno de ellos está destinado a usuarios.
Los modelos de estadísticas, por ejemplo, tienen la siguiente función y clase ocultas, pero no es utilizada por los modelos de estadísticas:
https://github.com/statsmodels/statsmodels/blob/master/statsmodels/miscmodels/try_mlecov.py#L36
Básicamente, si necesita una evaluación rápida, vuelva a escribirla para su caso de uso.
En el caso común de una matriz de covarianza diagonal, el PDF multivariable se puede obtener simplemente multiplicando los valores PDF univariados devueltos por una instancia de scipy.stats.norm . Si necesita el caso general, probablemente deba codificarlo usted mismo (lo que no debería ser difícil).
La densidad se puede calcular de una manera bastante directa usando las funciones numpy y la fórmula en esta página: http://en.wikipedia.org/wiki/Multivariate_normal_distribution . También es posible que desee utilizar la función de verosimilitud (probabilidad de registro), que es menos probable que desborde para grandes dimensiones y es un poco más fácil de calcular. Ambos solo implican poder calcular el determinante y el inverso de una matriz.
El CDF, por otro lado, es un animal completamente diferente ...
La normal multivariada ahora está disponible en 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])
Si todavía es necesario, mi implementación sería
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()
En caso de que haga cambios futuros, el código está aquí en GitHub
Utilizo el siguiente código que calcula el valor logpdf, que es preferible para dimensiones más grandes. También funciona para matrices 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)
El código es de pyParticleEst , si desea el valor pdf en lugar del logpdf simplemente tome math.exp () en el valor devuelto