1 ## Copyright (C) 2003 Iain Murray
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
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
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/>.
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}).
24 function s = mvnrnd(mu,Sigma,K)
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
29 % May 2004 take a third arg, cases. Makes it more compatible with Matlab's.
31 % Paul Kienzle <pkienzle@users.sf.net>
33 % * Add docs for argument K
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]))
48 if (size(Sigma)~=[d,d])
49 error('Sigma must have dimensions dxd where mu is nxd.');
55 [E,Lambda]=eig(Sigma);
56 if (min(diag(Lambda))<0),error('Sigma must be positive semi-definite.'),end
60 s = randn(n,d)*U + mu;
63 % {{{ END OF CODE --- Guess I should provide an explanation:
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
71 % Our new variable y is distributed according to:
72 % y ~ B*exp(-0.5*y'*inv(U'*U)*y)
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.
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
88 % [E,Lambda]=eig(Sigma)
93 % If any Lambdas are negative then Sigma just isn't even positive semi-definite
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.