generation of normally distributed random vector with covariance matrix

Question

In matlab it is easy to generate a normally distributed random vector with a mean and a standard deviation. From the help randn:

Generate values from a normal distribution with mean 1 and standard deviation 2. r = 1 + 2.*randn(100,1);

Now I have a covariance matrix C and I want to generate N(0,C).

But how could I do this?

From the randn help: Generate values from a bivariate normal distribution with specified mean vector and covariance matrix. mu = [1 2]; Sigma = [1 .5; .5 2]; R = chol(Sigma); z = repmat(mu,100,1) + randn(100,2)*R;

But I don't know exactly what they are doing here.

Which part do you not understand? mu is the mean vector (in your case 0, so leave it off), Sigma is the covariance matrix, and they're generating 100 pairs of random numbers. — Donnie
I don't understand the repmat part. Is it also possible to do it this way: chol(C, 'lower') + randn(N,1); with C the covariance matrix — Derk
Sorry, I think I do understand now. The repmat is used to build a mean matrix for 100 pairs of random numbers. — Derk
you can use the MVNRND function from the Statistics Toolbox, see this related question: stackoverflow.com/questions/4041866/gaussian-basis-function — Amro
@Donnie or Amro : you should post an answer so Derk can accept it and this question will be archived for anyone to consult :) — Bruno

Andrew Mao Andrew Mao · Accepted Answer · 2013-01-25T07:56:29

This is somewhat a math question, not a programming question. But I'm a big fan of writing great code that requires both solid math and programming knowledge, so I'll write this for posterity.

You need to take the Cholesky decomposition (or any decomposition/square root of a matrix) to generate correlated random variables from independent ones. This is because if X is a multivariate normal with mean m and covariance D, then Y = AX is a multivariate normal with mean Am and covariance matrix ADA' where A' is the transpose. If D is the identity matrix, then the covariance matrix is just AA' which you want to be equal to the covariance matrix C you are trying to generate.

The Cholesky decomposition computes such a matrix A and is the most efficient way to do it.

For more information, see: http://web.as.uky.edu/statistics/users/viele/sta601s03/multnorm.pdf

generation of normally distributed random vector with covariance matrix

2 Answers