title: Maximum Likelihood Estimation for Multivariate Normal Distribution in Python
body: <p>I am trying to fit the parameters of a multivariate normal distribution using MLE. </p>
  
  <pre><code>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
  </code></pre>
  
  <p>Fitting univariate normal distribution is fine:</p>
  
  <pre><code>obs = np.random.normal(1,4,50000)
  ini = [0,1]
  print(estimation(obs,lambda ob,p:norm.logpdf(ob,p[0],p[1]),ini))
  </code></pre>
  
  <p>But ran into some issues for multivariate (error assigning array to variable):</p>
  
  <pre><code>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))
  </code></pre>
  
  <p>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. </p>
  
  <pre><code>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))
  </code></pre>
  
  <p>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.</p>