Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

How would one sample from the $2$-dimensional normal distribution with mean $0$ and covariance matrix $$\begin{bmatrix} a & b\\b & c \end{bmatrix}$$ given the ability to sample from the standard ($1$-dimensional) normal distribution?

This seems like it should be quite simple, but I can't actually find an answer anywhere.

share|cite|improve this question
(This is not homework.) $\:$ – Ricky Demer Dec 31 '12 at 18:09
up vote 5 down vote accepted

Typically, given a covariance matrix $Q$, one perform the Cholesky decomposition on $Q$ i.e. get $Q = LL^T$ where $L$ is a lower triangular matrix.

Once we have $L$, generate a standard normal vector say $x$. Then $y=Lx$ is the desired random vector with covariance matrix $Q$.

This is because we have $$\mathbb{E}(y) = \mathbb{E}(Lx) = L\mathbb{E}(x) = L \times 0 = 0$$ and $$\mathbb{E}(yy^T)=\mathbb{E}(Lx(Lx)^T)=\mathbb{E}(Lxx^TL^T) = L\mathbb{E}(xx^T)L^T = L \times I \times L^T = LL^T$$

Note that any other symmetric square root decomposition of $Q$ also gives you the desired result i.e. if we can decompose $Q$ into $GG^T$, then $Gx$ will also gives us the random vector with zero mean and covariance $Q$.

In your case, if we perform Cholesky decomposition, we get $$L = \begin{bmatrix}\sqrt{a} & 0\\ \dfrac{b}{\sqrt{a}} & \dfrac{\sqrt{ac-b^2}}{\sqrt{a}} \end{bmatrix}$$ Now generate $x = \begin{bmatrix} x_1 \\ x_2\end{bmatrix}$ where $x_,x_2 \sim \mathcal{N}(0,1)$, then $$y = Lx = \begin{bmatrix}\sqrt{a} & 0\\ \dfrac{b}{\sqrt{a}} & \dfrac{\sqrt{ac-b^2}}{\sqrt{a}} \end{bmatrix} \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} \sqrt{a} x_1 \\ \dfrac{bx_1 + \sqrt{ac-b^2}x_2}{\sqrt{a}}\end{bmatrix}$$is the desired vector you are looking for.

share|cite|improve this answer
May I ask why in $L\mathbb{E}(xx^T)L^T$ above, $xx^T=I$? – user90593 Feb 21 '15 at 15:11

Wikipedia has an outline. You need a square root (of sorts) $A$ of your covariance matrix, which you can get using the Cholesky decomposition. Then sample two iid $N(0,1)$ random variables, multiply them (as a vector) by $A$, and you get a bivariate draw from your bivariate normal distribution.

If you work in R, you can use the mvtnorm package.

share|cite|improve this answer

Say you have a random variable $X\sim N(0,E)$ where $E$ is the identity matrix. Let $A$ be a matrix. Then $Y:=AX\sim N(0,AA^T)$. Hence you need to find a matrix $A$ with $AA^T = \left[\matrix{a & b \\ b & c}\right]$. There is no unique solution to this problem. One popular method is the Cholesky decomposition, where you find a triangular matrix $L$ with a given covariance matrix. Another method is to perform a principle axis transform $$ \left[\matrix{a & b \\ b & c}\right] = U^T\left[\matrix{\lambda_1 & 0 \\ 0 & \lambda_2}\right]U $$ with $UU^T=E$ and then take $$ A = U^T\left[\matrix{\sqrt{\lambda_1} & 0 \\ 0 & \sqrt{\lambda_2}}\right]U $$ as a solution. This is the only symmetric positive definite solution then. It is also called the square root of the positive-definite symmetric matrix $AA^T$.

More generally for distribution $X\sim N(\mu,\Sigma)$ it holds $AX+b\sim N(A\mu+b,B\Sigma B^T)$, see Wikipedia.

share|cite|improve this answer

Several methods have been posted; here's another.

If you can sample from the $2$-dimensional normal distribution whose variance is $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix}$, you can then rescale the $x$-coordinate and the $y$-coordinate to get the standard deviations you want, and that won't alter the correlation $\rho$. If $A$ is a square root of the given correlation matrix and $\begin{bmatrix} Z_1 \\ Z_1 \end{bmatrix}$ has a $2$-dimensional normal distribution whose variance is the $2\times 2$ identity matrix, then $A\begin{bmatrix} Z_1 \\ Z_1 \end{bmatrix}$ has a $2$-dimensional normal distribution whose variance is $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix}$.

So the remaining problem is finding a square root of $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix}$.

Let $P=\begin{bmatrix} 1/2 & 1/2 \\ 1/2 & 1/2\end{bmatrix}$ and $Q=I-P$. Then $P^2=P$, $Q^2=Q$, and $PQ=QP=0$.

Then we have $$\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix} = (1+\rho)P + (1-\rho)Q.$$ Since $P^2=P$, $Q^2=Q$, and $PQ=QP=0$, we can find a square root just by taking square roots of the coefficients: $$ \sqrt{1+\rho}\,P + \sqrt{1-\rho}\,Q. $$

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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