[SciPy-user] Bias in numpy.random.multivariate_normal?

Matthieu Brucher matthieu.brucher at gmail.com
Sat Nov 10 14:00:54 EST 2007


IIRC, you computed the biaised variance, so it is normal to see biais. Try
dividing by (sample_size-1), it should be the unbiaised estimator.

Matthieu

2007/11/10, John Reid <j.reid at mail.cryst.bbk.ac.uk>:
>
> Could someone confirm I've got this right?
>
> In a multivariate normal (MVN) the covariance is  E(X.X') - E(X).E(X').
> I sample many times from a MVN with given mean and covariance. I then
> calculate the covariance of the sample. This covariance matrix should
> show no bias to be larger or smaller than the distribution's covariance
> but I see that some entries are always larger and others are smaller.
>
> See this code:
>
> from numpy.random import multivariate_normal
> from numpy.linalg import inv
> from numpy import outer
>
> mu = [10.,10.]
> precision = [[.6, .3], [.4, .3]]
> sigma = inv(precision)
> sample_size = 10000
> num_tests = 100
>
> for test in xrange(num_tests):
>    samples = multivariate_normal(mu, sigma, [sample_size])
>    sample_mean = samples.sum(axis=0)/sample_size
>    #print (sample_mean - mu) > 0.
>    sample_covariance = sum(outer(sample-sample_mean, sample-sample_mean)
> for sample in samples) / sample_size
>    print (sample_covariance - sigma) > 0
>
>
> 'True' is printed when the entry is larger and False otherwise
>
> Did I get this wrong or is this acceptable behaviour for a random
> variate generator or is it a bug?
>
> Thanks,
> John.
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>



-- 
French PhD student
Website : http://miles.developpez.com/
Blogs : http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn : http://www.linkedin.com/in/matthieubrucher
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20071110/ddfcd15f/attachment.html>


More information about the SciPy-User mailing list