--- a/main/statistics/mvnpdf.m	Mon Nov 07 21:11:53 2005 +0000
+++ b/main/statistics/mvnpdf.m	Tue Nov 08 04:42:42 2005 +0000
@@ -1,6 +1,45 @@
-% y = mvnpdf(x,mu,sigma)
-% Compute multivariate normal pdf for x given mean mu and covariance matrix sigma.
-% x is dimension d x p, is 1 x p and sigma is p x p.
+% y = mvnpdf(x,mu,Sigma)
+% Compute multivariate normal pdf for x given mean mu and covariance matrix 
+% sigma.  The dimension of x is d x p, mu is 1 x p and sigma is p x p.
+%
+% 1/y^2 = (2 pi)^p |\Sigma| exp { (x-\mu)' inv(\Sigma) (x-\mu) }
+%
+% Ref: NIST Engineering Statistics Handbook 6.5.4.2
+% http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc542.htm
+
+% Algorithm
+%
+% 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 A = x' inv(sigma) x
+%	
+%    A = x' inv(sigma) x = x' inv(r' * r) x = x' inv(r) inv(r') x
+%
+% Given that inv(r') == inv(r)', at least in theory if not numerically,
+%
+%    A  = (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:
+%
+%    A = sumsq(x/r)
+%
+% and the final result is:
+%
+%    r = chol(Sigma);
+%    y = (2*pi)^(-p/2) * exp(-sumsq((x-mu)/r,2)/2) / prod(diag(r));
+
 
 % This program is public domain
 % Author: Paul Kienzle
@@ -10,10 +49,10 @@
   if nargin == 1, mu = 0; end
   if all(size(mu) == [1,p]), mu = repmat(mu,[d,1]); end
   if nargin < 3
-    pdf = (2*pi)^(-p/2) * exp(-0.5 * sumsq(x-mu,2));
+    pdf = (2*pi)^(-p/2) * exp(-sumsq(x-mu,2)/2);
   else
     r = chol(sigma);
-    pdf = (2*pi)^(-p/2) * exp(-0.5 * sum(((x-mu)/r).^2,2)) / prod(diag(r));
+    pdf = (2*pi)^(-p/2) * exp(-sumsq((x-mu)/r,2)/2) / prod(diag(r));
   end
 
 %!test