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?

  • $\begingroup$ 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). $\endgroup$ Dec 11 '17 at 18:57
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Dec 11 '17 at 19:46

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.