I am trying to fit the parameters of a multivariate normal distribution using MLE.
import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize
from scipy.stats import multivariate_normal as mnorm
def estimation(obs,fun,init,method='Nelder-Mead'):
mle = lambda param: -np.sum(fun(*[obs,param])) ## negate since we will minimize
result = minimize(mle,init,method=method)
return result.x
Fitting univariate normal distribution is fine:
obs = np.random.normal(1,4,50000)
ini = [0,1]
print(estimation(obs,lambda ob,p:norm.logpdf(ob,p[0],p[1]),ini))
But ran into some issues for multivariate (error assigning array to variable):
obs_m = np.random.multivariate_normal([0,0],[[1,0],[0,100]],50000)
ini_m = [[0,0],[[1,0],[0,100]]]
print(estimation(obs_m,lambda ob,p:mnorm.logpdf(ob,p[0],p[1],ini_m))
It seems the optimization algorithm doesn't work for arbitrary arrays/matrices. I have to unpack the mean array and covariance matrix into a flat array to feed to minimize.
ini_m = [0,0,1,0,0,100]
print(estimation(obs_m,lambda ob,p:mnorm.logpdf(ob,[p[0],p[1]],[[p[2],p[3]],[p[4],p[5]]]),ini_m))
Obviously this will quickly go out of hand when the dimension increases or for some more complicated distribution without closed form solution. What is the best to do here? Thank you.
