## [Octave-cvsupdate] octave-forge/main/statistics mvnpdf.m,1.1,1.2

 [Octave-cvsupdate] octave-forge/main/statistics mvnpdf.m,1.1,1.2 From: Paul Kienzle - 2005-11-08 04:42:53 Update of /cvsroot/octave/octave-forge/main/statistics In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29349 Modified Files: mvnpdf.m Log Message: Improve docs; cite ref. Index: mvnpdf.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/statistics/mvnpdf.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- mvnpdf.m 5 Nov 2005 16:57:46 -0000 1.1 +++ mvnpdf.m 8 Nov 2005 04:42:42 -0000 1.2 @@ -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 

 [Octave-cvsupdate] octave-forge/main/statistics mvnpdf.m,1.1,1.2 From: Paul Kienzle - 2005-11-08 04:42:53 Update of /cvsroot/octave/octave-forge/main/statistics In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29349 Modified Files: mvnpdf.m Log Message: Improve docs; cite ref. Index: mvnpdf.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/statistics/mvnpdf.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- mvnpdf.m 5 Nov 2005 16:57:46 -0000 1.1 +++ mvnpdf.m 8 Nov 2005 04:42:42 -0000 1.2 @@ -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