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' )