CODEWITHSUNDEEP
pythonmatplotlibscipystatistics
alberto
alberto

Reputation: 2645

Why bivariate_normal returns NaNs even if covariance is semi-positive definite?

I have the following normal distributed points:

import numpy as np
from matplotlib import pylab as plt
from matplotlib import mlab

mean_test = np.array([0,0])
cov_test = array([[ 0.6744121 , -0.16938146],
                  [-0.16938146,  0.21243464]])

The covariance matrix is definite semi-positive so it can be used as a covariance

# Semi-positive definite if all eigenvalues are 0 or 
# if there exists a Cholesky decomposition 
print np.linalg.eigvals(cov_test)
print np.linalg.cholesky(cov_test)

[ 0.72985988 0.15698686]

[[ 0.82122597 0. ] [-0.20625439 0.41218172]]

If I generate some points I get:

 data_test = np.random.multivariate_normal(mean_test, cov_test, 1000)
 plt.scatter(data_test[:,0],data_test[:,1])

data

Question:

Why does bivariate_normal method fail (return NaNs) when I try to plot the covariance contour? 

x = np.arange(-3.0, 3.0, 0.1)
y = np.arange(-3.0, 3.0, 0.1)
X, Y = np.meshgrid(x, y)
Z = mlab.bivariate_normal(X, Y, 
                      cov_test[0,0], cov_test[1,1],
                      0, 0, cov_test[0,1])
print Z
plt.contour(X, Y, Z)

Output:

 [[ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]
 ..., 
 [ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]]

 ValueError: zero-size array to reduction operation minimum which has no identity

Upvotes: 3

Views: 1228

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114911

The diagonals of the covariance matrix are the variances, but the arguments sigmax and sigmay of mlab.bivariate_normal are the square roots of the variances. Change this:

Z = mlab.bivariate_normal(X, Y, 
                      cov_test[0,0], cov_test[1,1],
                      0, 0, cov_test[0,1])

to this:

Z = mlab.bivariate_normal(X, Y, 
                      np.sqrt(cov_test[0,0]), np.sqrt(cov_test[1,1]),
                      0, 0, cov_test[0,1])

Upvotes: 5

Related Questions