Repository navigation
numpy.random.multivariate_normal accepts indefinite covariance matrices (Trac #1223) #1821
Description
Original ticket http://projects.scipy.org/numpy/ticket/1223 on 2009-09-15 by trac user zero79, assigned to unknown.
numpy.random.multivariate_normal will happily accept non-physical covariance matrices (i.e. non-symetric and those with det(C)<0), which results in erroneous distributions. attached is code that demonstrates the problem for a covariance matrix with det(C)<0 (except when rhoxz=1, which produces a valid covariance matrix).
!/usr/bin/python
import numpy , pylab
nmcsamples = 2_10_*6
xavg = 1.0
yavg = 2.0
zavg = 3.0
xdev = 0.1
ydev = 0.1
zdev = 0.1
rhoxy = 1.0
rhoxz = numpy.arange( 0.0 , 1.05 , 0.1 )
rhoyz = 1.0
nsamples = len( rhoxz )
xdevmc = numpy.zeros( nsamples )
ydevmc = numpy.zeros( nsamples )
zdevmc = numpy.zeros( nsamples )
rhoxymc = numpy.zeros( nsamples )
rhoxzmc = numpy.zeros( nsamples )
rhoyzmc = numpy.zeros( nsamples )
for n in range( 0 , nsamples ):
cxx = xdev**2
cyy = ydev**2
czz = zdev**2
cxy = rhoxy*xdev*ydev
cxz = rhoxz[n]*xdev*zdev
cyz = rhoyz*ydev*zdev
c = numpy.array( [ [ cxx , -cxy , cxz ] , [ cxy , cyy , cyz ] , [ cxz , cyz , czz ] ] )
x,y,z = numpy.transpose( numpy.random.multivariate_normal( [ xavg , yavg , zavg ] , c , nmcsamples ) )
print rhoxz[n],numpy.linalg.det( c )
xavgmc = numpy.mean( x )
yavgmc = numpy.mean( y )
zavgmc = numpy.mean( z )
xdevmc[n] = numpy.sqrt( numpy.mean( ( x - xavgmc )**2 ) )
ydevmc[n] = numpy.sqrt( numpy.mean( ( y - yavgmc )**2 ) )
zdevmc[n] = numpy.sqrt( numpy.mean( ( z - zavgmc )**2 ) )
rhoxymc[n] = ( numpy.mean( x*y ) - xavgmc*yavgmc )/xdevmc[n]/ydevmc[n]
rhoxzmc[n] = ( numpy.mean( x*z ) - xavgmc*zavgmc )/xdevmc[n]/zdevmc[n]
rhoyzmc[n] = ( numpy.mean( y*z ) - yavgmc*zavgmc )/ydevmc[n]/zdevmc[n]
pylab.subplots_adjust( wspace = 0.4 )
pylab.rc( 'font' , size = 9 )
pylab.rc( 'legend' , fontsize = 9 )
pylab.subplot( 2 , 3 , 1 )
pylab.plot( rhoxz , numpy.ones( nsamples )xdev , rhoxz , xdevmc )
pylab.xlabel( '$\rho{xz}$' )
pylab.ylabel( '$\Delta x$' )
pylab.axis( [ 0.0 , 1.0 , 0.09 , 0.12 ] )
pylab.subplot( 2 , 3 , 2 )
pylab.plot( rhoxz , numpy.ones( nsamples )ydev , rhoxz , ydevmc )
pylab.xlabel( '$\rho{xz}$' )
pylab.ylabel( '$\Delta y$' )
pylab.axis( [ 0.0 , 1.0 , 0.09 , 0.12 ] )
pylab.subplot( 2 , 3 , 3 )
pylab.plot( rhoxz , numpy.ones( nsamples )zdev , rhoxz , zdevmc )
pylab.xlabel( '$\rho{xz}$' )
pylab.ylabel( '$\Delta z$' )
pylab.legend( ( 'input' , 'monte carlo' ) )
pylab.axis( [ 0.0 , 1.0 , 0.09 , 0.12 ] )
pylab.subplot( 2 , 3 , 4 )
pylab.plot( rhoxz , numpy.ones( nsamples )rhoxy , rhoxz , rhoxymc )
pylab.xlabel( '$\rho{xz}$' )
pylab.ylabel( '$\rho_{xy}$' )
pylab.subplot( 2 , 3 , 5 )
pylab.plot( rhoxz , rhoxz , rhoxz , rhoxzmc )
pylab.xlabel( '$\rho_{xz}$' )
pylab.ylabel( '$\rho_{xz}$' )
pylab.subplot( 2 , 3 , 6 )
pylab.plot( rhoxz , numpy.ones( nsamples )*rhoyz , rhoxz , rhoyzmc )
pylab.xlabel( '$\rho_{xz}$' )
pylab.ylabel( '$\rho_{yz}$' )
pylab.savefig( 'multivariate_normal-problem.pdf' )
Activity
thouis commentedon Oct 19, 2012
Attachment added by trac user zero79 on 2009-09-15: multivariate_normal-problem
thouis commentedon Oct 19, 2012
@charris wrote on 2009-09-15
copy and paste formatting didn't work.
You need to highlight the code and hit the code block button.
Your example code works fine here. If there is a negative singular value the svd code is likely the culprit and that indicates a problem with lapack/atlas. Where did you get those packages? What os/distro are you running? What numpy version?
thouis commentedon Oct 19, 2012
trac user zero79 wrote on 2009-09-16
this code should execute just fine. the real problem is actually demonstrated in the plot file saved to multivariate_normal-problem.pdf.
this is not a problem with the svd code itself. its that non-physical covariances are permissable as inputs to multivariate_normal, where that should not be allowed. i've intentionally created bad covariance matrices (with negative determinants, and one with valid determinant, to demonstrate the problem). i would like to see numpy check for these cases and throw an error. otherwise, the user gets something that looks feasible, but is completely wrong.
fyi, this was tested on numpy 1.3.0 on debian lenny and sid.
thouis commentedon Oct 19, 2012
trac user dgoldsmith wrote on 2010-06-29
Let me be more succinct: the problem is not that the code doesn't run, the problem is that the code does run! Here's a much, much simpler example:
This shouldn't run because the covariance matrix is not "non-negative definite" (aka, positive semidefinite).
I discovered this bug independently while working on this function's docstring which states at the end of the Notes (i.e., almost as an afterthought):
"Note that the covariance matrix must be non-negative definite."
I thought it odd that such a requirement should be relegated to the Notes and not be included in the discussion of the cov function parameter, so I checked the code and found no "enforcement," so I checked the behavior, and found that, sure enough, the function is all too "happy" to accept input which the Notes state it "must" not be fed.
Now, I can see the reason - speed: my first call of this function w/ a two-component integer mean and a very simple integer-valued cov resulted in a brief but noticeable delay in being provided the answer; is this due to the use of svd, which a note in the code says should be replaced by Cholesky decomp.? - why we may want to leave it up to the user to enforce positive semidefiniteness (PSDness), but IMHO that's a dangerous road to go down, and if we do go down it, then the warning sign needs to be featured much more prominently, i.e., in the Extended summary and again in the cov Parameter description.
However, a better alternative, IMO, is to add code that checks for PSDness and raises a warning if it's detected - that should be the default behavior, but we can add a keyword parameter which allows the user to by-pass that check if they want to.
thouis commentedon Oct 19, 2012
Title changed from
need errors for non-physical numpy.random.multivariate_normal covariance matricestonumpy.random.multivariate_normal accepts indefinite covariance matricesby trac user dgoldsmith on 2010-06-29thouis commentedon Oct 19, 2012
@vincentdavis wrote on 2010-06-29
I agree with this. (including adding the keyword)
"However, a better alternative, IMO, is to add code that checks for PSDness and raises a warning if it's detected - that should be the default behavior, but we can add a keyword parameter which allows the user to by-pass that check if they want to. "
thouis commentedon Oct 19, 2012
@josef-pkt wrote on 2010-06-29
some discussion is in the thread at
http://mail.scipy.org/pipermail/numpy-discussion/2009-September/045255.html
thouis commentedon Oct 19, 2012
Milestone changed to
Unscheduledby @mwiebe on 2011-03-24thouis commentedon Oct 19, 2012
trac user dgoldsmith wrote on 2010-06-30
Thanks, josef. Based on that discussion, might I suggest a hierarchical approach: first, Cholesky is tried to take the matrix square-root - according to CH in the above thread, IIUC, this will automatically provide the function w/ a check on the "physicality" of
cov. If Cholesky works, great, we're done; if it throws an exception, to satisfy RK, an independent check oncovis performed (what would be the speediest way to independently assesscovfor both symmetry and PSDness?). Ifcovpasses that test, eigh (as suggested by CH) is tried; if it succeeds, the result is returned, along w/ a message that Cholesky failed,covmight be asymmetric and/or not PSD; if eigh also fails, give up: return np.empty() and a message thatcovis likely asymmetric and/or not PSD.rgommers commentedon Jan 4, 2017
Fixed by gh-5726, so closing.