# Given a multivariate normal distribution, how can we simulate uniform random variables that hold on to the correlation structure?

Suppose that I have a multivariate normal distribution such that:

$$\mathbf{X} \sim \mathcal{N}\left(\mathbf{0}, \Sigma \right)$$

where $\mathbf{X}$ is an $n \times 1$ vector and I am wondering how I can generate uniform random variables that are correlated based on this. Would it as simple as taking $\mathbf{X}$ vector, and applying a standard normal transformation to it?

• Well, just applying (the appropriate) normal cumulativ distribution function to the components would not preserve the correlation structure. You should look into copulas (search this site). Dec 11 '17 at 18:57
• One method that uses the suggestion by @Kjetil is described at stats.stackexchange.com/questions/257779/… (for the case of two variables, where the objective is to establish a linear regression). I believe the technique would generalize to this question. Simply transforming the variables will not reproduce $\Sigma$ exactly, except in special cases.
– whuber
Dec 11 '17 at 19:46

## 1 Answer

Wikipedia introduces the method of drawing values from the distribution. I am now trying to illustrate it in more details.

If we want to draw sample from the multivariate normal distribution specified by $\mu \in \mathbb{R}^{n\times 1}$ and $\Sigma \in \mathbb{R}^{n\times n}$ and $\Sigma^\top = \Sigma$. We can do the following to generate random variable $\mathbf{x} \in \mathbb{R}^{n\times 1}$ such that $\mathbf{x} \sim \mathcal{N}(\mu, \Sigma)$.

1. Find matrix $A$ such that $A A^\top = \Sigma$. When $\Sigma$ is positive definite, the Cholesky decomposition is often used. Other way of getting such matrix $A$ includes SVD decomposition.
2. When we have such matrix $A$, we just use some random number generator, e.g. randn(n, 1) function in Matlab to generate $n$ independent standard normal variables and get a random vector $\mathbf{z} \in \mathbb{R}^{n \times 1}$, and $\mathbf{x} = \mu + A \mathbf{z}$ is what we want.

Julia code for this is

n = 2                  # 2 dimension situation
μ = [1.; 2.]           # mean value
Σ = [exp(0) exp(-1);
exp(-1) exp(0)]
Σ = Symmetric(Σ)
A = chol(Σ)'           # chol(Σ) = A^T
z = randn(2, 2000)     # generate 2000 points
x = A * z .+ μ


We now verify that data generated in this way satisfies that $$\mathbf{x} \sim \mathcal{N}(\mu, \Sigma)$$ As $\mathbf{z} \sim \mathcal{N}(0, 1)$, from affine transformation of multivariate normal distribution we know that $\mathbf{x} = \mu + A \mathbf{z}$ indeed satisfies the needed property.