NovaDenizen

asked on

# Want to generate random samples consistent with covariance matrix

Say I have a fairly big nontrivial covariance matrix (50x50 or bigger) and I would like to generate a bunch of vectors (say, 10000) whose covariance is consistent with this matrix. How would I go about that?

Assume the means of all the result elements are zero, for simplicity's sake. Optimizing the incremental time for each vector generated is more important than optimizing the preparatory caclculations. Also assume I can generate normally distributed variables with no problem.

I'm guessing you could do it by generating a bunch of independent normally distributed variables, then assigning each element in the result vector as a linear combination of the independent variables, but the exact method is escaping me.

A link to a good (and free) paper or description of the method would also be acceptable.

Assume the means of all the result elements are zero, for simplicity's sake. Optimizing the incremental time for each vector generated is more important than optimizing the preparatory caclculations. Also assume I can generate normally distributed variables with no problem.

I'm guessing you could do it by generating a bunch of independent normally distributed variables, then assigning each element in the result vector as a linear combination of the independent variables, but the exact method is escaping me.

A link to a good (and free) paper or description of the method would also be acceptable.

ASKER

Bueller? Bueller? Bueller?

ASKER

Ok, I'll post my answer and ask that the question be cancelled.

Lets say we have two normally distributed variables defined as linear combinations of independent normally distributed variables.

n1 and n2 are normally distributed random variables with mean=0 and sigma=1

x = a*n1 + b*n2

y = c*n1 + d*n2

So, what is the variance of x? Look up the normal sum distribution and you'll see that's well known to be a^2*var(n1) + b^2*var(n2), which is a^2 + b^2.

What is the covariance of x and y? cov(x,y) = <x * y> - <x> * <y>

(the <...> operator is used to denote a mean)

<x> and <y> are both zero, so

cov(x,y) = <x*y>

= <(a*n1 + b*n2) * (c*n1 + d*n2)>

= <(a*c*n1^2) + (b*d*n2^2) + (ad + bc)*n1*n2>

= a*c*<n1^2> + b*d*<n2^2> + (ad + bc)*<n1*n2>

since <n1^2> is equivalent to the variance of n1, which is 1:

= a*c + b*d + (ad + bc)*<n1*n2>

<n1*n2>, since the means of n1 and n2 are zero, is equivalent to the covariance of n1 and n2. Since we know n1 and n2 are independent, their covariance must be zero.

cov(x,y) = a*c + b*d

Let's use this information to construct two random variables with the covariance matrix C=

[ c11 c21 ]

[ c21 c22 ]

where the variables x are defined as linear combinations of the independent n variables (mean=0 sigma=1):

x = A * n

[ x1 ] = [a11 a12] * [n1]

[ x2 ] [a21 a22] [n2]

We know that var(x2) = c22, so as long as we come up with a21 and a22 such that a21^2 + a22^2 = c22, we have satisfied that element in the covariance matrix. We can just assign a21 = 0 and a22 = sqrt(c22).

So, now how do we define a11 and a12? We know the covariance between x2 and x1 is c21, equal to a21*a11 + a22*a12. Since a21=0, a12 = c21 / a22.

cov(x1) = c11 = a11^2 + a12^2, so a11 = sqrt(c11 - a12^2).

This procedure can be obviously expanded to any size matrix and number of variables. It will result in an upper-triangular A matrix, which just so happens to have the property:

A * transpose(A) = C

Which means that generating the A coefficient matrix is exactly equivalent to generating a Cholesky decomposition of the covariance matrix C.

Lets say we have two normally distributed variables defined as linear combinations of independent normally distributed variables.

n1 and n2 are normally distributed random variables with mean=0 and sigma=1

x = a*n1 + b*n2

y = c*n1 + d*n2

So, what is the variance of x? Look up the normal sum distribution and you'll see that's well known to be a^2*var(n1) + b^2*var(n2), which is a^2 + b^2.

What is the covariance of x and y? cov(x,y) = <x * y> - <x> * <y>

(the <...> operator is used to denote a mean)

<x> and <y> are both zero, so

cov(x,y) = <x*y>

= <(a*n1 + b*n2) * (c*n1 + d*n2)>

= <(a*c*n1^2) + (b*d*n2^2) + (ad + bc)*n1*n2>

= a*c*<n1^2> + b*d*<n2^2> + (ad + bc)*<n1*n2>

since <n1^2> is equivalent to the variance of n1, which is 1:

= a*c + b*d + (ad + bc)*<n1*n2>

<n1*n2>, since the means of n1 and n2 are zero, is equivalent to the covariance of n1 and n2. Since we know n1 and n2 are independent, their covariance must be zero.

cov(x,y) = a*c + b*d

Let's use this information to construct two random variables with the covariance matrix C=

[ c11 c21 ]

[ c21 c22 ]

where the variables x are defined as linear combinations of the independent n variables (mean=0 sigma=1):

x = A * n

[ x1 ] = [a11 a12] * [n1]

[ x2 ] [a21 a22] [n2]

We know that var(x2) = c22, so as long as we come up with a21 and a22 such that a21^2 + a22^2 = c22, we have satisfied that element in the covariance matrix. We can just assign a21 = 0 and a22 = sqrt(c22).

So, now how do we define a11 and a12? We know the covariance between x2 and x1 is c21, equal to a21*a11 + a22*a12. Since a21=0, a12 = c21 / a22.

cov(x1) = c11 = a11^2 + a12^2, so a11 = sqrt(c11 - a12^2).

This procedure can be obviously expanded to any size matrix and number of variables. It will result in an upper-triangular A matrix, which just so happens to have the property:

A * transpose(A) = C

Which means that generating the A coefficient matrix is exactly equivalent to generating a Cholesky decomposition of the covariance matrix C.

ASKER CERTIFIED SOLUTION

membership

Create a free account to see this answer

Signing up is free and takes 30 seconds.

**No credit card required.**
ASKER

I'll leave this open for a couple days and award anybody else who figures it out.