From help-request at octave dot org Sat Nov  5 17:46:18 2005
Subject: Re: Multivariate pdf of a normal distribution
From: Mike Miller <mbmiller at taxa dot epi dot umn dot edu>
To: Paul Kienzle <pkienzle at users dot sourceforge dot net>
cc: Octave Help <help at octave dot org>, brumen@hispeed.ch
Date: Sat, 5 Nov 2005 17:38:37 -0600 (CST)

On Sat, 5 Nov 2005, Paul Kienzle wrote:

>> Of course, the density goes to infinity when Sigma is singular.  Is your 
>> use of chol() just meant to check that the matrix is PD?
>
> The wikipedia entry on cholesky decomposition 
> (http://en.wikipedia.org/wiki/Choleskey_decomposition) suggests it is faster 
> and more stable than the lu decomposition which would be used to compute the 
> inverse.  The speed doesn't matter in this case, but accuracy is always a 
> concern.  The side effect of checking positive definiteness of sigma is a 
> bonus.
>
> Some numerical tests with e.g.,
>
> 	n=11; x = prolate(n); cn=cond(x), d=norm(x*inv(x)-eye(n)), r=chol(x); 
> c=norm(x*inv(r)*inv(r)'-eye(n)),
>
> shows that it chol() is indeed a little better than inv() for ill-conditioned 
> positive definite matrices. The function prolate() is from higham's test 
> matrix toolbox.
>
> Similarly for hilb(), though norm(inv(x) - invhilb(n)) and 
> norm(inv(r)*inv(r)' - invhilb(n)) are both pretty bad.


Interesting.  Thanks for sharing this information.

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