octave_packages/statistics-1.1.3/mvnpdf.m

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