Drawing a normal vector

Section of Estimation and inference for the multivariate normal, part of Statistical learning

Type	Section
Taught in	Estimation and inference for the multivariate normal
Part	Statistical learning

The distributions just stated are also the distributions a simulation should reproduce. A generator usually supplies independent uniforms, or independent standard normals, and the step from those draws to \operatorname{N}_p(\boldsymbol{\mu},\Sigma) is the linear-image proposition. Start from uniforms by the inverse transform. Let U be uniform on [0,1] and set X=\Phi^{-1}(U), with \Phi the standard normal distribution function. Because \Phi is increasing,

\mathbb{P}(X\le x)=\mathbb{P}\bigl(\Phi^{-1}(U)\le x\bigr)=\mathbb{P}\bigl(U\le\Phi(x)\bigr)=\Phi(x),

so X is standard normal. This is the inverse-transform method. Draw p independent copies and stack them as \mathbf{Z}\sim\operatorname{N}_p(\mathbf{0},I_p). Any matrix B with BB^\top=\Sigma then has the covariance required by the linear-image formula. The spectral factor B=\Gamma\Lambda^{1/2}, whose columns are eigenvectors scaled by square-root eigenvalues, is one such factor, and the Cholesky factor is another. Set

\mathbf{X}=\boldsymbol{\mu}+B\mathbf{Z}.

The mean of B\mathbf{Z}+\boldsymbol{\mu} is \boldsymbol{\mu} and its covariance is BB^\top=\Sigma. Normality follows from the linear-image proposition. In software, rnorm supplies pseudorandom draws. A seed reproduces the computation. A sample covariance close to its target checks this implementation. It is not a proof that every property of the generator is correct.

The spectral and Cholesky factors both satisfy BB^\top=\Sigma. Such a factor need not satisfy B^2=\Sigma. The name symmetric square root is reserved for \Gamma\Lambda^{1/2}\Gamma^\top. A third construction never forms B. Draw the first coordinate from its univariate marginal normal, then each later coordinate from the univariate conditional normal given the coordinates already drawn, using the conditional formula on Covariance and the multivariate normal. The constructions target the same \operatorname{N}_p(\boldsymbol{\mu},\Sigma). They factor it differently.

A goodness-of-fit test for normality does not decide whether the exact F calibration below applies. With small n the test has little power, and a rejection does not name a replacement model. The exact results on this page assume normality in the statements that use it. The next section records what remains for the sample mean when that assumption is dropped and only a finite covariance is kept.

The fragment below is the Cholesky route. The lower-triangular factor satisfies BB^\top=\Sigma, with \boldsymbol{\mu}=(1,2)^\top and \Sigma equal to the matrix with entries 4,2 in the first row and 2,3 in the second. The column means and the sample covariance come out close to that target. Setting a seed before rnorm makes a later run reproduce the same draws. The laboratory at the end of the page evaluates the one-sample test on data, rather than checking a simulated covariance.

Sigma <- matrix(c(4, 2, 2, 3), 2, 2); mu <- c(1, 2)
B <- t(chol(Sigma))                 # lower-triangular factor, B %*% t(B) = Sigma
Z <- matrix(rnorm(2 * 1000), 2)     # columns are independent N(0, I) vectors
X <- t(mu + B %*% Z)                # rows are N(mu, Sigma) draws
colMeans(X); cov(X)                 # close to mu and Sigma

Drawing a normal vector

Used in

See also