[Numpy-discussion] strange behavior of numpy.random.multivariate_normal, ticket:1842

[Robert Kern]

On Thu, Feb 16, 2012 at 16:12, Pierre Haessig <pierre.haessig at crans.org> wrote:
> Le 16/02/2012 16:20, josef.pktd at gmail.com a écrit :
>
>> I don't see any way to fix multivariate_normal for this case, except
>> for dropping svd or for random perturbing a covariance matrix with
>> multiplicity of singular values.
>
> Hi,
> I just made a quick search in what R guys are doing. It happens there are
> several codes (http://cran.r-project.org/web/views/Multivariate.html ). For
> instance, mvtnorm
> (http://cran.r-project.org/web/packages/mvtnorm/index.html). I've attached
> the related function from the source code of this package.
>
> Interestingly enough, it seems they provide 3 different methods (svd, eigen
> values, and Cholesky).
> I don't have the time now to dive in the assessments of pros and cons of
> those three. Maybe one works for our problem, but I didn't check yet.

The main reason I used the SVD variant is because the Cholesky
decomposition failed on some covariance matrices that were nearly not
positive definite (i.e. had a nearly-0 eigenvalue). In the application
that I extracted this code from, this was a valid thing to do; the
deviates just inhabit an infinitely thin subspace of the main space,
but are otherwise multivariate-normally-distributed in that subspace.

I'm not too attached to the semantics. We should check that the
Cholesky decomposition is stable before switching, though. The
eigenvalue algorithm probably suffers from instability just as much as
the SVD one.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco