octave_packages/statistics-1.1.3/mvnrnd.m

   1 ## Copyright (C) 2003 Iain Murray
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} @var{s} = mvnrnd (@var{mu}, @var{Sigma})
  18 ## @deftypefnx{Function File} @var{s} = mvnrnd (@var{mu}, @var{Sigma}, @var{n})
  19 ## Draw @var{n} random @var{d}-dimensional vectors from a multivariate Gaussian
  20 ## distribution with mean @var{mu}(@var{n}x@var{d}) and covariance matrix
  21 ## @var{Sigma}(@var{d}x@var{d}).
  22 ## @end deftypefn
  23
  24 function s = mvnrnd(mu,Sigma,K)
  25
  26   % Iain Murray 2003 -- I got sick of this simple thing not being in Octave and
  27   %                     locking up a stats-toolbox license in Matlab for no good
  28   %                     reason.
  29   % May 2004 take a third arg, cases. Makes it more compatible with Matlab's.
  30
  31   % Paul Kienzle <pkienzle@users.sf.net>
  32   % * Add GPL notice.
  33   % * Add docs for argument K
  34
  35   % If mu is column vector and Sigma not a scalar then assume user didn't read
  36   % help but let them off and flip mu. Don't be more liberal than this or it will
  37   % encourage errors (eg what should you do if mu is square?).
  38   if ((size(mu,2)==1)&&(size(Sigma)~=[1,1]))
  39     mu=mu';
  40   end
  41
  42   if nargin==3
  43     mu=repmat(mu,K,1);
  44   end
  45
  46   [n,d]=size(mu);
  47
  48   if (size(Sigma)~=[d,d])
  49     error('Sigma must have dimensions dxd where mu is nxd.');
  50   end
  51
  52   try
  53     U=chol(Sigma);
  54   catch
  55     [E,Lambda]=eig(Sigma);
  56     if (min(diag(Lambda))<0),error('Sigma must be positive semi-definite.'),end
  57     U = sqrt(Lambda)*E';
  58   end
  59
  60   s = randn(n,d)*U + mu;
  61 endfunction
  62
  63 % {{{ END OF CODE --- Guess I should provide an explanation:
  64 %
  65 % We can draw from axis aligned unit Gaussians with randn(d)
  66 %   x ~ A*exp(-0.5*x'*x)
  67 % We can then rotate this distribution using
  68 %   y = U'*x
  69 % Note that
  70 %   x = inv(U')*y
  71 % Our new variable y is distributed according to:
  72 %   y ~ B*exp(-0.5*y'*inv(U'*U)*y)
  73 % or
  74 %   y ~ N(0,Sigma)
  75 % where
  76 %   Sigma = U'*U
  77 % For a given Sigma we can use the chol function to find the corresponding U,
  78 % draw x and find y. We can adjust for a non-zero mean by just adding it on.
  79 %
  80 % But the Cholsky decomposition function doesn't always work...
  81 % Consider Sigma=[1 1;1 1]. Now inv(Sigma) doesn't actually exist, but Matlab's
  82 % mvnrnd provides samples with this covariance st x(1)~N(0,1) x(2)=x(1). The
  83 % fast way to deal with this would do something similar to chol but be clever
  84 % when the rows aren't linearly independent. However, I can't be bothered, so
  85 % another way of doing the decomposition is by diagonalising Sigma (which is
  86 % slower but works).
  87 % if
  88 %   [E,Lambda]=eig(Sigma)
  89 % then
  90 %   Sigma = E*Lambda*E'
  91 % so
  92 %   U = sqrt(Lambda)*E'
  93 % If any Lambdas are negative then Sigma just isn't even positive semi-definite
  94 % so we can give up.
  95 %
  96 % Paul Kienzle adds:
  97 %   Where it exists, chol(Sigma) is numerically well behaved.  chol(hilb(12))
  98 %   for doubles and for 100 digit floating point differ in the last digit.
  99 %   Where chol(Sigma) doesn't exist, X*sqrt(Lambda)*E' will be somewhat
 100 %   accurate.  For example, the elements of sqrt(Lambda)*E' for hilb(12),
 101 %   hilb(55) and hilb(120) are accurate to around 1e-8 or better.  This was
 102 %   tested using the TNT+JAMA for eig and chol templates, and qlib for
 103 %   100 digit precision.
 104 % }}}
 105