help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Multivariate pdf of a normal distribution


From: Mike Miller
Subject: Re: Multivariate pdf of a normal distribution
Date: Sat, 5 Nov 2005 23:35:45 -0600 (CST)

On Sat, 5 Nov 2005, Paul Kienzle wrote:

Thanks.  I've included mvnrnd in octave-forge.

Looking at the code it uses chol(), and if that doesn't work, it uses eig().

The R manual for the MASS package has this to say about mvrnorm:

The matrix decomposition is done via eigen; although a Choleski decomposition might be faster, the eigendecomposition is stabler.

        http://stat.ethz.ch/R-manual/R-patched/library/MASS/html/mvrnorm.html

Using the same test I did earlier on ill-conditioned positive definite matrices, eig doesn't seem to be a more accurate way to compute matrix inverses. Anyone care to comment on this?


I'll just say that I like a system where chol() is tried, and if it fails, eig() is tried. Usually chol() will work and it will be fast, but if a matrix is singular or negative definite, chol() will fail. When chol() fails, eig() should show a zero eigenvalue, but it might be estimated as a negative eigenvalue. You will then see complex eigenvectors and you will produce complex mv-normal random numbers, which is not what you want. Therefore, when the zero eigenvalue is estimated as negative, you have to take the real part of the square root of the eigenvalues and of the eigenvector matrix and use only those real parts. But when the negative eigenvalue is truly a negative number, and not just a poorly estimated zero, the function must fail, but where should we draw the line?

Here is a function my student, Gregg Lind, wrote a few days ago:


function U=eigensplit(M),

try
   U = chol(M);
catch
   [E,Lam] = eig(M);
   U = sqrt(Lam)*E';
   U = real(U);
end


I think more checking would be wise. Any opinions on how best to write such a function? What would be a good name for it?

I think the kinship package in R (also GPL) has a generalized cholesky for singular matrices, but it might only work when there is a correlation of 1.0 off the diagonal.

Thanks again, Paul, for your great work.

Mike



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------



reply via email to

[Prev in Thread] Current Thread [Next in Thread]