= Solution
Let $Z=(Z_1,\ldots,Z_p)^T$. Independence of the standard <normal distributions> gives $\mathbb EZ=0$, $\operatorname{Cov}(Z)=I_p$, and joint density
$$
f_Z(z)=(2\pi)^{-p/2}\exp(-z^Tz/2).
$$
Choose a matrix $C$ with $CC^T=\Sigma$, for example the lower-triangular factor from the <Cholesky decomposition>. It exists and is nonsingular because $\Sigma$ is a <symmetric positive-definite matrix>. Then the required simulation is
$$
\boxed{X=\mu+CZ.}
$$
It has mean $\mu$ and <covariance matrix> $CC^T=\Sigma$, and every linear combination of its coordinates is Gaussian, so it has the required <multivariate normal distribution>.
For the <multivariate normal density>, apply the <change of variables formula> to $z=C^{-1}(x-\mu)$. Its absolute Jacobian determinant is $|\det C|^{-1}=|\Sigma|^{-1/2}$, and
$$
z^Tz=(x-\mu)^TC^{-T}C^{-1}(x-\mu)=(x-\mu)^T\Sigma^{-1}(x-\mu).
$$
Hence
$$
\boxed{f_X(x)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}\exp\left[-\frac12(x-\mu)^T\Sigma^{-1}(x-\mu)\right],\qquad x\in\mathbb R^p.}
$$
The nonsingularity assumption is needed for a density with respect to $p$-dimensional Lebesgue measure.
Back to article page
