# Using scipy.minimize to find the maximum likelihood estimates for multivariate gaussian

Let's say I have a 100x2 normally distributed array of data.

   In [2]: import numpy as np
...:
...: # instantiate a random number generator
...: rng = np.random.default_rng(100)
...:
...: # define mu and sigma for the dummy sample
...: mu = np.array([0.5, 0.25])
...: cov = np.array([[1, 0.5],
...:                 [0.5, 1]])
...:
...: # generate multivariate normal random sample
...: x = rng.multivariate_normal(mu, cov, size=100)
...:
...: print(x.shape)
(100, 2)


I want to find the maximum likelihood estimates of parameters $$\vec{\mu}$$ and $$\Sigma$$ using the scipy minimize function. Finding the maxima of the log-likelihood is equivalent to finding the minima of the $$-\log(\mathcal{L})$$.

$$-\log(\mathcal{L}) = -l(\vec{\mu}, \Sigma) = \frac{1}{2}[nk\ln(2\pi) + n\ln(\det(\Sigma^{-1})) + \sum_{i = 1}^{n}(\vec{x} - \vec{\mu})^{T}\Sigma^{-1}(\vec{x}-\vec{\mu})]$$

def log_likelihood(params, x):
mu, covmat = params
print(mu, covmat)
x = x[:, np.newaxis] # add a new first dimension to x
k = mu.shape[0]  # number of dimensions
n = x.shape[1]  # number of samples
inv_covmat = np.linalg.inv(covmat)  # inverse of the covariance matrix
diff = x - mu  # deviation of x from the mean
maha_dist = np.einsum('ijk, kl, ijl->ij', diff, inv_covmat, diff) # mahalabonis distance

term1 = n*k*np.log(2*np.pi)
term2 = n*np.log(np.linalg.det(inv_cov))
term3 = np.sum(maha_dist, axis=0)
return 0.5 * (term1 + term2 + term3)


My initial idea of how to minimise this likelihood function with respect to the parameters using scipy.optimize is shown below.

from scipy.optimize import minimize

# initial guesses
mu_guess = [0, 0]
covmat_guess = [[1, 0],
[0, 1]]

lik_model = minimize(log_likelihood, x0=(mu_guess, covmat_guess), args=x)

lik_model.x


But I am getting the following error. I presume it has something to do with the shapes of my initial guesses, but I don't know how else I can represent them.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: float() argument must be a string or a real number, not 'list'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Input In [20], in <cell line: 6>()
2 mu_guess = [0, 0]
3 covmat_guess = [[1, 0],
4                 [0, 1]]
----> 6 lik_model = minimize(log_likelihood, x0=(mu_guess, covmat_guess), args=x)
8 lik_model.x

File ~/miniforge3/envs/minimise/lib/python3.10/site-packages/scipy/optimize/_minimize.py:687, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
685     res = _minimize_cg(fun, x0, args, jac, callback, **options)
686 elif meth == 'bfgs':
--> 687     res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
688 elif meth == 'newton-cg':
689     res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
690                              **options)

File ~/miniforge3/envs/minimise/lib/python3.10/site-packages/scipy/optimize/_optimize.py:1296, in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, finite_diff_rel_step, **unknown_options)
1293 if maxiter is None:
1294     maxiter = len(x0) * 200
-> 1296 sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
1297                               finite_diff_rel_step=finite_diff_rel_step)
1299 f = sf.fun
1300 myfprime = sf.grad

File ~/miniforge3/envs/minimise/lib/python3.10/site-packages/scipy/optimize/_optimize.py:263, in _prepare_scalar_function(fun, x0, jac, args, bounds, epsilon, finite_diff_rel_step, hess)
259     bounds = (-np.inf, np.inf)
261 # ScalarFunction caches. Reuse of fun(x) during grad
262 # calculation reduces overall function evaluations.
--> 263 sf = ScalarFunction(fun, x0, args, grad, hess,
264                     finite_diff_rel_step, bounds, epsilon=epsilon)
266 return sf

File ~/miniforge3/envs/minimise/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:107, in ScalarFunction.__init__(self, fun, x0, args, grad, hess, finite_diff_rel_step, finite_diff_bounds, epsilon)
101     raise ValueError("Whenever the gradient is estimated via "
102                      "finite-differences, we require the Hessian "
103                      "to be estimated using one of the "
104                      "quasi-Newton strategies.")
106 # the astype call ensures that self.x is a copy of x0
--> 107 self.x = np.atleast_1d(x0).astype(float)
108 self.n = self.x.size
109 self.nfev = 0

ValueError: setting an array element with a sequence.



Thanks