Re: Multivariate pdf of a normal distribution

[Paul Kienzle]

On Nov 5, 2005, at 5:10 PM, Mike Miller wrote:

> On Sat, 5 Nov 2005, Paul Kienzle wrote:
>
> > The math looks pretty easy to implement:
> >
> >         http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc542.htm
> >
> > Using Cholesky factorization on the positive definite covariance 
> > matrix:
>
> why not this?:
>
> density = mvnorm_pdf(x, mu, Sigma)
>
> # first check array sizes, but I'm skipping that
>
> p = length(x);
>
> density = (2*pi)^(-p/2) * (1/sqrt(det(Sigma))) * 
> exp(-.5*(x-mu)'*inv(Sigma)*(x-mu));
>
>
> 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.

- Paul



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