# Sampling from a $2$d normal with a given covariance matrix

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.

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.

• May I ask why in $L\mathbb{E}(xx^T)L^T$ above, $xx^T=I$? – user90593 Feb 21 '15 at 15:11
• Nice answer. But where comes the role of mean in this case ? What if mean is not zero ? – Shyamkkhadka Oct 21 '16 at 13:24
• @user90593 - we assume that x is standard normal, which means it's variance is 1. It also means that its expectancy is 0. V(x) = E(x^2) - E(x)^2 . So E(x^2) = 1 – ihadanny Jun 11 '17 at 12:08

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.

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.

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.$$