On Thu, Feb 16, 2012 at 16:12, Pierre Haessig <pierre.haes...@crans.org> wrote: > Le 16/02/2012 16:20, josef.p...@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 _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
