Distribución t de estudiantes multivariada con python

Para generar muestras con distribución t multivariable utilizo esta función:

def multivariatet(mu,Sigma,N,M): ''' Output: Produce M samples of d-dimensional multivariate t distribution Input: mu = mean (d dimensional numpy array or scalar) Sigma = scale matrix (dxd numpy array) N = degrees of freedom M = # of samples to produce ''' d = len(Sigma) g = np.tile(np.random.gamma(N/2.,2./N,M),(d,1)).T Z = np.random.multivariate_normal(np.zeros(d),Sigma,M) return mu + Z/np.sqrt(g) 

pero lo que busco ahora es la auto -distribución de estudiantes multivariable para poder calcular la densidad de los elementos donde dimension > 1 .

Eso será algo así como stats.t.pdf(x, df, loc, scale) del paquete scipy pero en un espacio multidimensional.

Codifiqué la densidad por mí mismo:

 import numpy as np from math import * def multivariate_t_distribution(x,mu,Sigma,df,d): ''' Multivariate t-student density: output: the density of the given element input: x = parameter (d dimensional numpy array or scalar) mu = mean (d dimensional numpy array or scalar) Sigma = scale matrix (dxd numpy array) df = degrees of freedom d: dimension ''' Num = gamma(1. * (d+df)/2) Denom = ( gamma(1.*df/2) * pow(df*pi,1.*d/2) * pow(np.linalg.det(Sigma),1./2) * pow(1 + (1./df)*np.dot(np.dot((x - mu),np.linalg.inv(Sigma)), (x - mu)),1.* (d+df)/2)) d = 1. * Num / Denom return d 

Esto evalúa el registro pdf de la distribución multivariable de Student-T para n por d matriz de datos X:

 from scipy.special import gamma from numpy.linalg import slogdet def multivariate_student_t(X, mu, Sigma, df): #multivariate student T distribution [n,d] = X.shape Xm = X-mu V = df * Sigma V_inv = np.linalg.inv(V) (sign, logdet) = slogdet(np.pi * V) logz = -gamma(df/2.0 + d/2.0) + gamma(df/2.0) + 0.5*logdet logp = -0.5*(df+d)*np.log(1+ np.sum(np.dot(Xm,V_inv)*Xm,axis=1)) logp = logp - logz return logp 

Generalicé el código de @ farhawa para permitir múltiples entradas en x (encontré que quería consultar múltiples puntos a la vez).

 import numpy as np from math import gamma def multivariate_t_distribution(x, mu, Sigma, df): ''' Multivariate t-student density. Returns the density of the function at points specified by x. input: x = parameter (nd numpy array; will be forced to 2d) mu = mean (d dimensional numpy array) Sigma = scale matrix (dxd numpy array) df = degrees of freedom Edited from: http://stackoverflow.com/a/29804411/3521179 ''' x = np.atleast_2d(x) # requires x as 2d nD = Sigma.shape[0] # dimensionality numerator = gamma(1.0 * (nD + df) / 2.0) denominator = ( gamma(1.0 * df / 2.0) * np.power(df * np.pi, 1.0 * nD / 2.0) * np.power(np.linalg.det(Sigma), 1.0 / 2.0) * np.power( 1.0 + (1.0 / df) * np.diagonal( np.dot( np.dot(x - mu, np.linalg.inv(Sigma)), (x - mu).T) ), 1.0 * (nD + df) / 2.0 ) ) return 1.0 * numerator / denominator 

Probé las respuestas anteriores y obtuve diferentes resultados de cada una, pero no estoy seguro de por qué / qué podría estar mal. Lo siguiente, que me basé en el código scikit-learn para mezclas gaussianas, creo que funciona (para matrices de entrada de tamaño arbitrario X, y c distribuciones t con parámetros en listas de medios y covars):

 import numpy as np from scipy import linalg try: # SciPy >= 0.19 from scipy.special import gammaln as sp_gammaln except ImportError: from scipy.misc import gammaln as sp_gammaln def log_multivariate_t_density(X, means, covars, nu = 1): n_samples, n_dim = X.shape nmix = len(means) log_prob = np.empty((n_samples, nmix)) for c, (mu, cv) in enumerate(zip(means, covars)): try: cv_chol = linalg.cholesky(cv, lower=True) except linalg.LinAlgError: try: cv_chol = linalg.cholesky(cv + min_covar * np.eye(n_dim), lower=True) except linalg.LinAlgError: raise ValueError("'covars' must be symmetric, " "positive-definite") cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol))) cv_sol = linalg.solve_triangular(cv_chol, (X - mu).T, lower=True).T norm = (sp_gammaln((nu + n_dim) / 2.) - sp_gammaln(nu / 2.) - 0.5 * n_dim * np.log(nu * np.pi)) inner = - (nu + n_dim) / 2. * np.log1p(np.sum(cv_sol ** 2, axis=1) / nu) log_prob[:, c] = norm + inner - cv_log_det return log_prob