1
$\begingroup$

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

$\endgroup$

1 Answer 1

1
$\begingroup$

The problem actually was the problem.

I have since learned that the covariance matrix was to be assumed as a constant and it is in fact only the mean which is to be determined.

Removing this requirement also removes the exception detailed above.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.