Martin Jansche wrote:
Does this mean you want to sample from the multivariate normal
distribution? If so, you're right that there is no code in the GSL
(except for the bivariate case), unless I'm missing something. To
simulate a random draw from the multivariate normal, consult any
statistics book that discusses computational techniques. For example,
Gelman, Carlin, Stern & Rubin ("Bayesian Data Analysis", 1995, p.478)
give the following algorithm. The basic idea is to use the Cholesky
decomposition of the covariance matrix (GSL can do that) in
conjunction with several draws from the univariate normal (GSL can do
that too).
More specifically, gsl_linalg_cholesky_decomp(Sigma) computes the
Cholesky decomposition of the covariance matrix Sigma and stores it in
the lower triangular part of Sigma. Call that matrix (the lower
triangular one) A. Perform d independent draws from a standard normal
(using gsl_ran_ugaussian()) and store them in a d-dimensional vector
x. Then transform x into a new vector y as by letting y = mu + A x,
where mu is the d-dimensional mean vector. The d-dimensional vector y
is then a draw from the mulitvariate normal with covariance matrix
Sigma and mean vector mu.
Does that answer your question?
- martin
Dear Martin,
Thank you very much! You focused exactly the problem, and now I'm going
to try the code you have outlined. Using the GSL gives to me much
control over the generated data, i.e. the possibility to select a random
seed for the generation of the x vector, and better error control.
Unfortunately sometimes the covariance matrix is not positive-definite,
i.e. because two variables are too-much correlated, and the
decomposition fails. I've read that the correct approach to this problem
should be to reduce the number of variables by 1, i.e. from 6 to 5,
removing one of the two too-much correlated variables, and then use the
reduced covariance matrix to generate 5-vectors instead of 6-vectors.
The only thing I don't understand is how to correctly recover, after the
data generation, a 6-vector from a 5-vector; I suppose I should use a
formula relating the mu, the standard deviation of the removed variable
and the generated 5-vector to the 6th element. Maybe there is a better
solution to this. The solution to this issue would definitely solve the
problem.
Best regards
--Pasquale