by Iago González Rodríguez
Last Updated August 14, 2019 00:19 AM - source

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:

- Serverfault Help
- Superuser Help
- Ubuntu Help
- Webapps Help
- Webmasters Help
- Programmers Help
- Dba Help
- Drupal Help
- Wordpress Help
- Magento Help
- Joomla Help
- Android Help
- Apple Help
- Game Help
- Gaming Help
- Blender Help
- Ux Help
- Cooking Help
- Photo Help
- Stats Help
- Math Help
- Diy Help
- Gis Help
- Tex Help
- Meta Help
- Electronics Help
- Stackoverflow Help
- Bitcoin Help
- Ethereum Help