trying to plot contours of bivariate normal, won't work with a correlation term

Question

refer to this tutorial: http://matplotlib.org/1.4.0/examples/pylab_examples/contour_demo.html

Here is the prototype for the bivariate_normal function from mplotlib.mlab:

bivariate_normal(X, Y, sigmax=1.0, sigmay=1.0, mux=0.0, muy=0.0, sigmaxy=0.0)

X and Y define the grid, and we have arguments for the 2 dimensional means and covariance terms. As you can see, there is an argument at the end for a the covariance between x and y. Here's the thing: plt.contour() will plot bivariate normal contours if sigmaxy = 0. However, if sigmaxy has any other value, I get a

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

For example,

Z =  bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0, 0.0)
plt.contour(X,Y,Z)

works

But, the following does not work:

Z = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0, 1.0)
plt.contour(X,Y,Z)

Anyone familiar with matplotlib have any ideas? Thanks.

Answer 1

To take the notation from the accepted answer:

cov_test = np.array([[..., ...],
                     [..., ...]])

You need to pass the standard deviation (in the arguments sigmax and sigmay), not the variance:

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

For whatever reason, the sigmaxy argument is the actual entry from the covariance matrix (actually rho*sigmax*sigmay).

For any nontrivial example, the code will blow up when it tries to calculate rho if you try to pass in either all covariance or all sqrt covariance.

Answer 2

It does not work because your covariance matrix is not positive definite. To see if your matrix is positive definite you can check whether all its eigenvalues are greater than zero.

Extreme case

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


cov_test = np.array([[1,0.999],
                     [0.999,1]])

print np.linalg.eigvals(cov_test)

[ 1.99900000e+00 1.00000000e-03]

You can see the second eigenvalue is super close to zero. Actually, if you plot it, you see this is really an extreme case of covariance:

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 = bivariate_normal(X, Y, cov_test[0,0], cov_test[1,1], 0.0, 0.0, cov_test[0,1])
plt.contour(X,Y,Z)

Non positive definite case:

And if you go a bit further...

import numpy as np
from matplotlib import pylab as plt

cov_test = np.array([[1,1],
                 [1,1]])

print np.linalg.eigvals(cov_test)

[ 2. 0.]

Then the second eigenvalue reaches 0, this is not positive definite and if you try to plot it then:

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 = bivariate_normal(X, Y, cov_test[0,0], cov_test[1,1], 0.0, 0.0, cov_test[0,1])
plt.contour(X,Y,Z)

You get the error:

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

Actually, my Z is now full of NaN