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: Paul Kienzle
Subject: Re: Multivariate pdf of a normal distribution
Date: Sat, 5 Nov 2005 16:24:47 -0500

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:

        r = chol(sigma);

where r'*r = sigma. Being upper triangular, the determinant of r is trivially the product of the diagonal, and the determinant of sigma is the square of this:

        det = prod(diag(r))^2;

The formula asks for the square root of the determinant, so no need to square it.

The exponential argument y = x' inv(sigma) x
        
        y = x' inv(sigma) x = x' inv(r' * r) x = x' inv(r) inv(r') x

I don't know for sure that inv(r') == inv(r)' for r upper triangular. Numerically it is not the case in octave:

        octave:34> x = triu(rand(10)); norm(inv(x') - inv(x)')
        ans =  3.3466e-14

Assuming that it is, then

        y = (x' / r) * (x'/r)' = sumsq(x'/r)

The interface takes the parameters to the mvn in columns rather than rows, so we are actually dealing with the transpose:

        y = sumsq(x/r)

and the final result is:

    r = chol(sigma);
pdf = (2*pi)^(-p/2) * exp(-0.5 * sum(((x-mu)/r).^2,2)) / prod(diag(r));

I've added mvnpdf to octave-forge (see attached).

Note: Google found a cdf here:

http://alex.strashny.org/a/Multivariate-normal-cumulative- distribution-function-(cdf)-in-MATLAB.html


- Paul


On Nov 5, 2005, at 8:18 AM, Gorazd Brumen wrote:


Hi!

Is there a function in octave or octave forge that gives the density of
a *multivariate* normal distribution. normal_pdf only gives a one-dimensional.
I know it is a joke to do get it from there, but nevertheless.

Attachment: mvnpdf.m
Description: Binary data



reply via email to

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