We typically create a p-variate multivariate normal vector by Cholesky
factorizing the var-covar matrix and multiplying the resultant by an
iid normal(0,1) vector:
mvnorm = chol(cov)'*rand(p,1);
or for N mv-normal row vectors:
mvnorm = rand(N,p)*chol(cov);
That works great unless the variance covariance matrix is singular.
For that case, I have the answer given in comments at the end of the
page at this URL:
http://www.gatsby.ucl.ac.uk/~iam23/code/mvnrnd.m
% But the Cholesky decomposition function doesn't always work...
% Consider Sigma=[1 1;1 1]. Now inv(Sigma) doesn't actually exist, but
Matlab's
% mvnrnd provides samples with this covariance st x(1)~N(0,1)
x(2)=x(1). The
% fast way to deal with this would do something similar to chol but be
clever
% when the rows aren't linearly independent. However, I can't be
bothered, so
% another way of doing the decomposition is by diagonalising Sigma
(which is
% slower but works).
% if
% [E,Lambda]=eig(Sigma)
% then
% Sigma = E*Lambda*E'
% so
% U = sqrt(Lambda)*E'
% If any Lambdas are negative then Sigma just isn't even positive
semi-definite
% so we can give up.
That seems like the perfect solution until we realize that numerical
imprecision causes the lambdas to be small negative numbers when they
really should be zeros:
[e,lambda]=eig(ones(3)); sqrt(lambda(1))
ans = 0.00000000000000e+00 + 3.25207016166203e-09i
Thus, checking for negative lambdas will sometimes fail us but not
checking for them will cause some complex-valued mv-normal data, which
is very undesirable. What do you guys think is the best way of
dealing with this problem? I would like to return an error if the
matrix has a clearly negative eigenvalue, but not if the eigenvalue
was an imperfectly estimated zero. When the zero eigenvalue is
estimated as a small negative number, I could either use
sqrt(abs(lambda)) instead of sqrt(lambda) or maybe
real(sqrt(lambda)). Or maybe I could replace the near-zero values
with zeros. I'm not concerned about speed for this step.
Any suggestions?
I would also love to know if any of you can suggest a way to produce a
generalized Cholesky factorization for non-negative definite matrices
that work when the matrix is singular. For singular non-ND matrices,
I think there must be infinitely many possible upper triangular
matrices, U, such that U'U returns the original non-ND matrix. I'd be
happy with any one of those solutions so long as it contains only real
values!
Mike
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------