I am working in a implementation of the Expectation-Maximization algorithm with Missing Data for Mixture of MVNs. You don't have to know anything about this algorithm to help me with my issue.
Let say that my dataset has shape
D x N with
D = 6.
I compute the estimation of sigma (the covariance matrix). The result is a
6 x 6 matrix (everything ok).
Then I have to compute the PDF of some samples for the distribution with the mean and covariance matrix that I just computed (estimated parameters).
In order to compute the PDF, I use
scipy.stats.multivariate_normal because it has a pdf() method.
But I always get the following error:
ValueError: the input matrix must be positive semidefinite
This error is due to the covariance matrix. I have checked its eigenvalues, and the first eigenvalue is always a negative value. Some real examples that I got:
import numpy as np eigvals = np.linalg.eigvals(sigma) # shape (6,) ''' Three different results that I got in different executions (the parameters are randomly generated) [-406.73080893 94.43623712 57.06170498 73.75668673 69.21443981 70.60878445] [-324.74509361 104.75315794 50.43203113 67.92014983 63.35285505 65.41071698] [-277.14957619 98.17501755 69.8623394 54.49958827 59.4808295 57.28174734] '''
Does anyone know why this is happening?
I am not very familiar with the eigenvalues, so I would greatly appreciate your help.
Here I include some code of how I compute sigma:
# This is how I compute the covariance matrix def estimate_sigma_with_nan(xx_est, gamma_k, mu): sigma = (xx_est / gamma_k) - np.dot(mu, np.transpose(mu)) return sigma # = (d, d) # This is how I compute xx_est. I will include a image of the maths behind this exp_prod = np.zeros_like(xx_est_k) x_i_norm[m] = estimate_x_norm(x_i_norm, mu_k, sigma_k, m, o) exp_prod[np.ix_(m, m)] = estimate_xx_norm(reshape_(x[:d, i], keep_axes=range(2)), mu_k, sigma_k, m, o) exp_prod[np.ix_(o, o)] = np.dot(x_i_norm[o], np.transpose(x_i_norm[o])) exp_prod[np.ix_(o, m)] = np.dot(x_i_norm[o], np.transpose(x_i_norm[m])) exp_prod[np.ix_(m, o)] = np.dot(x_i_norm[m], np.transpose(x_i_norm[o])) xx_est = xx_est + (exp_prod * gamma[k, i]) def estimate_xx_norm(x_i, mu_k, sigma_k, m, o): assert x_i.ndim == 2 and mu_k.ndim == 2 and sigma_k.ndim == 2 est_xx = sigma_k[np.ix_(m, m)] - np.dot(np.dot(sigma_k[np.ix_(m, o)], np.linalg.inv(sigma_k[np.ix_(o, o)])), np.transpose(sigma_k[np.ix_(m, o)])) est_xx = est_xx + np.dot(estimate_x_norm(x_i, mu_k, sigma_k, m, o), np.transpose(estimate_x_norm(x_i, mu_k, sigma_k, m, o))) assert est_xx.ndim == 2 return est_xx
This work is based on this paper: