 Solution

Let $Z=(Z_1,\dots ,Z_p{)}^T$. Independence of the standard normal distributions gives $\mathbb{E}Z=0$, $\mathrm{Cov}(Z)=I_p$, and joint densityChoose a matrix $C$ with $C{C}^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 isIt has mean $\mu$ and covariance matrix $C{C}^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}$, andHenceThe nonsingularity assumption is needed for a density with respect to $p$-dimensional Lebesgue measure.