TAOCP 3.4.1 Exercise 31

Exercise 31. [HM39] (Direct generation of normal deviates.)

a) Prove that if $a_1^2 + \cdots + a_n^2 = 1$ and if $X_1, \ldots, X_n$ are independent normal deviates with mean 0 and variance 1, then $a_1 X_1 + \cdots + a_n X_n$ is a normal deviate with mean 0 and variance 1.

b) The result of (a) suggests that we can generate new normal deviates from old ones, just as we obtain new uniform deviates from old ones. For example, we might use the idea of 3.2.2–(7), but with a recurrence like

$$X_n = (X_{n-24} + X_{n-55})/\sqrt{2} \qquad \text{or} \qquad X_n = \tfrac{1}{2}X_{n-21} + \tfrac{1}{2}X_{n-55},$$

after a set of normal deviates $X_1, \ldots, X_{54}$ has been initialized. Explain why this is not a good idea.

c) Show, however, that there is a suitable way to generate normal deviates quickly from other normal deviates, by using a refinement of the idea in (a) and (b). [Hint: If $X$ and $Y$ are independent normal deviates, so are $X' = X\cos\theta + Y\sin\theta$ and $Y' = -X\sin\theta + Y\cos\theta$, for any angle $\theta$.]

Verified: yes
Solve time: 4m21s

(a)

Let

$$ Y=a_1X_1+\cdots+a_nX_n, $$

where $X_1,\ldots,X_n$ are independent $N(0,1)$ random variables and

$$ a_1^2+\cdots+a_n^2=1. $$

We compute the moment generating function of $Y$. Since the $X_k$ are independent,

$$ M_Y(t)=E(e^{tY}) =E!\left(e^{t(a_1X_1+\cdots+a_nX_n)}\right) =\prod_{k=1}^n E(e^{ta_kX_k}). $$

For a standard normal variable $X$,

$$ E(e^{sX})=e^{s^2/2}. $$

Hence

$$ M_Y(t) =\prod_{k=1}^n e^{(a_k t)^2/2} =e^{\frac{t^2}{2}\sum_{k=1}^n a_k^2} =e^{t^2/2}. $$

But $e^{t^2/2}$ is exactly the moment generating function of a standard normal variable $N(0,1)$. Therefore $Y$ is normally distributed with mean $0$ and variance $1$.

The same conclusion can be verified directly from means and variances:

$$ E(Y)=\sum_{k=1}^n a_kE(X_k)=0, $$

and, using independence,

$$ \operatorname{Var}(Y) =\sum_{k=1}^n a_k^2\operatorname{Var}(X_k) =\sum_{k=1}^n a_k^2 =1. $$

Thus $Y\sim N(0,1)$.

(b)

Each recurrence preserves the one-dimensional distribution.

For example,

$$ X_n=\frac{X_{n-24}+X_{n-55}}{\sqrt2} $$

is a linear combination of two independent $N(0,1)$ variables with coefficients whose squares sum to $1$. By part (a), $X_n$ is again $N(0,1)$. The same is true for

$$ X_n=\frac12X_{n-21}+\frac12X_{n-55}, $$

after normalization of the variance.

The difficulty is not the marginal distribution but the dependence structure.

Every new value is a deterministic linear combination of earlier values. Consequently all later variables lie in the linear span of the initial values

$$ X_1,\ldots,X_{54}. $$

No new randomness is ever introduced. After many steps, strong correlations develop among the generated values. In fact, if

$$ \mathbf X=(X_1,\ldots,X_{54})^T, $$

every subsequent deviate has the form

$$ X_n=\sum_{j=1}^{54} c_{n,j}X_j. $$

Thus the entire infinite sequence is determined by only $54$ independent normal variables.

This is analogous to the defect of additive generators for uniform random numbers. Although each output individually has the desired distribution, the vectors

$$ (X_n,X_{n+1},\ldots) $$

occupy a low-dimensional linear structure and exhibit substantial correlations. Such a sequence is unsuitable when independent normal deviates are required.

(c)

The hint suggests using orthogonal transformations rather than one-dimensional recurrences.

Suppose $X$ and $Y$ are independent $N(0,1)$ variables. Define

$$ X' = X\cos\theta + Y\sin\theta, $$

$$ Y' = -X\sin\theta + Y\cos\theta. $$

In matrix form,

$$ \begin{pmatrix} X'\ Y' \end{pmatrix} = R \begin{pmatrix} X\ Y \end{pmatrix}, \qquad R= \begin{pmatrix} \cos\theta & \sin\theta\ -\sin\theta & \cos\theta \end{pmatrix}. $$

Since $R$ is orthogonal,

$$ RR^T=I. $$

The vector $(X,Y)$ has the two-dimensional standard normal distribution with density

$$ \frac1{2\pi}e^{-(x^2+y^2)/2}, $$

which depends only on $x^2+y^2$. Orthogonal transformations preserve $x^2+y^2$, hence preserve this distribution. Therefore $(X',Y')$ has the same joint distribution as $(X,Y)$. In particular, $X'$ and $Y'$ are again independent $N(0,1)$ variables.

This leads to a practical method. Maintain a pool of independent normal deviates. Repeatedly select pairs $(X,Y)$ and replace them by

$$ (X',Y') = (X\cos\theta+Y\sin\theta,, -X\sin\theta+Y\cos\theta), $$

using angles $\theta$ chosen independently. Each transformation is orthogonal, so it preserves the joint standard normal distribution and does not reduce the dimension of the underlying randomness. The resulting variables remain independent standard normals.

More generally, if

$$ \mathbf X=(X_1,\ldots,X_n)^T $$

has independent $N(0,1)$ components and $Q$ is any orthogonal matrix, then

$$ \mathbf X' = Q\mathbf X $$

also has independent $N(0,1)$ components. Applying sparse orthogonal transformations, such as successive random plane rotations, provides a fast way to generate large numbers of normal deviates from an existing pool while preserving the correct joint distribution.

$\square$