/usr/share/octave/packages/statistics-1.3.0/mvnpdf.m is in octave-statistics 1.3.0-4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 | ## Author: Paul Kienzle <pkienzle@users.sf.net>
## This program is granted to the public domain.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{y} =} mvnpdf (@var{x})
## @deftypefnx{Function File} {@var{y} =} mvnpdf (@var{x}, @var{mu})
## @deftypefnx{Function File} {@var{y} =} mvnpdf (@var{x}, @var{mu}, @var{sigma})
## Compute multivariate normal pdf for @var{x} given mean @var{mu} and covariance matrix
## @var{sigma}. The dimension of @var{x} is @var{d} x @var{p}, @var{mu} is
## @var{1} x @var{p} and @var{sigma} is @var{p} x @var{p}. The normal pdf is
## defined as
##
## @example
## @iftex
## @tex
## $$ 1/y^2 = (2 pi)^p |\Sigma| \exp \{ (x-\mu)^T \Sigma^{-1} (x-\mu) \} $$
## @end tex
## @end iftex
## @ifnottex
## 1/@var{y}^2 = (2 pi)^@var{p} |@var{Sigma}| exp @{ (@var{x}-@var{mu})' inv(@var{Sigma})@
## (@var{x}-@var{mu}) @}
## @end ifnottex
## @end example
##
## @strong{References}
##
## NIST Engineering Statistics Handbook 6.5.4.2
## http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc542.htm
##
## @strong{Algorithm}
##
## Using Cholesky factorization on the positive definite covariance matrix:
##
## @example
## @var{r} = chol (@var{sigma});
## @end example
##
## where @var{r}'*@var{r} = @var{sigma}. Being upper triangular, the determinant
## of @var{r} is trivially the product of the diagonal, and the determinant of
## @var{sigma} is the square of this:
##
## @example
## @var{det} = prod (diag (@var{r}))^2;
## @end example
##
## The formula asks for the square root of the determinant, so no need to
## square it.
##
## The exponential argument @var{A} = @var{x}' * inv (@var{sigma}) * @var{x}
##
## @example
## @var{A} = @var{x}' * inv (@var{sigma}) * @var{x}
## = @var{x}' * inv (@var{r}' * @var{r}) * @var{x}
## = @var{x}' * inv (@var{r}) * inv(@var{r}') * @var{x}
## @end example
##
## Given that inv (@var{r}') == inv(@var{r})', at least in theory if not numerically,
##
## @example
## @var{A} = (@var{x}' / @var{r}) * (@var{x}'/@var{r})' = sumsq (@var{x}'/@var{r})
## @end example
##
## The interface takes the parameters to the multivariate normal in columns rather than
## rows, so we are actually dealing with the transpose:
##
## @example
## @var{A} = sumsq (@var{x}/r)
## @end example
##
## and the final result is:
##
## @example
## @var{r} = chol (@var{sigma})
## @var{y} = (2*pi)^(-@var{p}/2) * exp (-sumsq ((@var{x}-@var{mu})/@var{r}, 2)/2) / prod (diag (@var{r}))
## @end example
##
## @seealso{mvncdf, mvnrnd}
## @end deftypefn
function pdf = mvnpdf (x, mu = 0, sigma = 1)
## Check input
if (!ismatrix (x))
error ("mvnpdf: first input must be a matrix");
endif
if (!isvector (mu) && !isscalar (mu))
error ("mvnpdf: second input must be a real scalar or vector");
endif
if (!ismatrix (sigma) || !issquare (sigma))
error ("mvnpdf: third input must be a square matrix");
endif
[ps, ps] = size (sigma);
[d, p] = size (x);
if (p != ps)
error ("mvnpdf: dimensions of data and covariance matrix does not match");
endif
if (numel (mu) != p && numel (mu) != 1)
error ("mvnpdf: dimensions of data does not match dimensions of mean value");
endif
mu = mu (:).';
if (all (size (mu) == [1, p]))
mu = repmat (mu, [d, 1]);
endif
if (nargin < 3)
pdf = (2*pi)^(-p/2) * exp (-sumsq (x-mu, 2)/2);
else
r = chol (sigma);
pdf = (2*pi)^(-p/2) * exp (-sumsq ((x-mu)/r, 2)/2) / prod (diag (r));
endif
endfunction
%!demo
%! mu = [0, 0];
%! sigma = [1, 0.1; 0.1, 0.5];
%! [X, Y] = meshgrid (linspace (-3, 3, 25));
%! XY = [X(:), Y(:)];
%! Z = mvnpdf (XY, mu, sigma);
%! mesh (X, Y, reshape (Z, size (X)));
%! colormap jet
%!test
%! mu = [1,-1];
%! sigma = [.9 .4; .4 .3];
%! x = [ 0.5 -1.2; -0.5 -1.4; 0 -1.5];
%! p = [ 0.41680003660313; 0.10278162359708; 0.27187267524566 ];
%! q = mvnpdf (x, mu, sigma);
%! assert (p, q, 10*eps);
|