3. Random Number Generation

Every random variate in a simulation is built from i.i.d. U(0,1) numbers; these notes cover where those numbers come from and what makes a generator good.

Motivation. Consider simulating an $M/M/1$ queue (Markovian, i.e. exponential, interarrival times; exponential service times; one server). We must repeatedly generate $X \sim \text{Exp}(\lambda)$, whose c.d.f. is $F(x) = 1 - e^{-\lambda x}$, $x \ge 0$. This can be done from a single uniform random number:

1. Generate $U \sim U(0,1)$.
2. Return $X = -\dfrac{1}{\lambda}\ln U$.

Why it works: setting $U = F(X)$ and solving for $X$ gives $X = F^{-1}(U) = -\frac{1}{\lambda}\ln(1-U)$; since $1-U \sim U(0,1)$ as well, we may use $\ln U$ instead. Then

$$P(X \le x) = P\left(-\tfrac{1}{\lambda}\ln U \le x\right) = P\left(U \ge e^{-\lambda x}\right) = 1 - e^{-\lambda x} = F(x).$$

The same idea applies to other distributions: every random variate in a simulation is built from i.i.d. $U(0,1)$ numbers.

Goal of this chapter. Produce a sequence of i.i.d. $U(0,1)$ random numbers, i.e. numbers with p.d.f.

$$f(x) = \begin{cases} 1, & 0 \le x \le 1, \\ 0, & \text{otherwise.} \end{cases}$$

Three alternative sources of “random” numbers

1. True random numbers (from physical experiments)

• Draw numbered balls from an urn.
• Electrical circuits: a hardware random number machine, e.g. ERNIE (Electronic Random Number Indicator Equipment).
• Disadvantages:
  • Not reproducible. Reproducibility matters in simulation: it makes debugging easier, and it gives sharper comparisons when two alternative systems are run on the same random inputs.
  • Difficult to implement.
  • Slow.

2. Pseudo-random numbers — a completely deterministic sequence that is statistically indistinguishable from a true random sequence. This is what simulation uses in practice, and what the rest of the chapter is about.

3. Quasi-random numbers (used in Monte Carlo integration rather than in simulation)

An integral over $[0,1]$ is an expectation under $U(0,1)$, because the $U(0,1)$ p.d.f. equals $1$ there:

$$\int_0^1 g(x)\,dx = \int_0^1 g(x)\cdot 1\,dx = \mathbf{E}[g(X)], \qquad X \sim U(0,1).$$

So generate $X_1, \dots, X_n \overset{\text{i.i.d.}}{\sim} U(0,1)$ and estimate the integral by

$$\bar{g}(n) = \frac{1}{n}\sum_{i=1}^n g(X_i) \;\longrightarrow\; \mathbf{E}[g(X)].$$

By the CLT the error of this estimator is of order $1/\sqrt{n}$. The idea of quasi-random numbers is to replace the random sequence by a deterministic sequence that covers $[0,1]$ more evenly; this improves the error rate to order $(\ln n)/n$. Such sequences are not i.i.d., so they are useless as simulation inputs, but for integration the even coverage is exactly what helps.

3.1. Random Number Generators (RNGs)

Desired properties of an RNG

1. Uniformity: the sequence should appear to be uniformly distributed on $(0,1)$.
2. Independence: terms in the sequence should not be correlated.
3. Reproducibility: one must be able to reproduce a particular stream of random numbers (e.g. from a seed).
4. Fast, with low memory usage.
5. Long period. A deterministic generator eventually loops. The segment of non-repeating numbers is called a cycle, and its length is the period of the RNG.

flowchart LR
    u1(("u₁")) --> u2(("u₂")) --> u3(("u₃")) --> d[". . ."] --> uP(("u_P")) --> u1

The sequence repeats after $P$ terms, so the period must be far longer than the number of random numbers a simulation will consume.

3.2. Linear Congruential Generators (LCGs)

An LCG produces integers by the recursion

$$z_i = (a\, z_{i-1} + c) \bmod m, \qquad i = 1, 2, \dots$$

where $z_0$ is the initial seed, $a$ is the multiplier, $c$ is the increment, and $m$ is the modulus. All four are non-negative integers, and $x \bmod y$ is the remainder when $x$ is divided by $y$.

3.2.1. From integers to uniforms.

Since $z_i$ is a remainder after division by $m$,

$$z_0, z_1, z_2, \dots \in \{0, 1, \dots, m-1\}.$$

Dividing by $m$ scales these into the unit interval:

$$u_i = \frac{z_i}{m} \in [0, 1), \qquad i = 0, 1, 2, \dots$$

The $u_i$ are the “random numbers” delivered to the simulation; the $z_i$ are the generator’s internal state.

Example. Take $m = 16$, $a = 5$, $c = 3$, $z_0 = 7$:

$i$	0	1	2	3	4	$\cdots$	16
$z_i$	7	6	1	8	11	$\cdots$	7
$u_i$	0.4375	0.375	0.0625	0.5	0.6875	$\cdots$	0.4375

e.g. $z_1 = (5\cdot 7 + 3) \bmod 16 = 38 \bmod 16 = 6$. Here $z_{16} = z_0$, so the sequence has period $16 = m$.

Observations.

1. Values on a grid. The $u_i$ can only take the $m$ values $\left\{0, \tfrac{1}{m}, \tfrac{2}{m}, \dots, \tfrac{m-1}{m}\right\}$: rational numbers on an equally spaced grid, not a genuinely continuous $U(0,1)$ variable. In practice one takes $m$ large, typically $m \ge 10^9$; then the grid is so fine that $u_i \approx U(0,1)$ for all practical purposes.
2. Looping behaviour. Since $z_i$ depends only on $z_{i-1}$, as soon as some earlier value reappears the whole sequence after it repeats. There are only $m$ possible values, so a value must recur within $m$ steps. The length of the repeating segment is the period of the LCG, and by this argument the period is at most $m$. An LCG whose period equals $m$ is called a full-period LCG: it visits every value $0, 1, \dots, m-1$ exactly once per cycle.
3. Determinism. The whole sequence is fixed by $(a, c, m, z_0)$. Unrolling the recursion gives the closed form (for $a \ne 1$)

   $$z_i = \left[a^i z_0 + c\,\frac{a^i - 1}{a - 1}\right] \bmod m,$$

   so $z_i$ can be computed directly from the seed without generating $z_1, \dots, z_{i-1}$.

Example. $m = 16$, $a = 5$, $c = 3$, $z_0 = 7$, so $z_i = (5 z_{i-1} + 3) \bmod 16$:

$i$	0	1	2	3	4	$\cdots$	16	17
$z_i$	7	6	1	8	11	$\cdots$	7	6
$u_i$	0.438	0.375	0.063	0.5	0.688	$\cdots$	0.438	0.375

The first repeat is $z_{16} = z_0 = 7$, so this LCG has full period $16$.

What if $z_0 = 6$ or $z_0 = 1$ instead? The recursion is the same, so we get the same cycle entered at a different point — e.g. $6, 1, 8, 11, \dots, 7, 6, \dots$ — again with period $16$. In a full-period LCG all $m$ values lie on a single cycle, so any choice of seed produces the entire cycle, in some order; the seed only picks the starting point.

For a non-full-period LCG, the values split into several cycles and the period depends on the seed:

Example. $c = 0$, $a = 13$, $m = 64$, so $z_i = 13 z_{i-1} \bmod 64$.

a. $z_0 = 4$: $\;4 \to 52 \to 36 \to 20 \to 4$, period $4$.
b. $z_0 = 2$: $\;2 \to 26 \to 18 \to 42 \to 34 \to 58 \to 50 \to 10 \to 2$, period $8$.

Notice that in (a) every value is a multiple of $4$, and in (b) every value is $\equiv 2 \pmod 4$: multiplying by $13$ never changes the power of $2$ dividing $z_i$, so the sequence is trapped in its residue class. This is exactly what the theorem below rules out.

3.2.2. Necessary and sufficient conditions for a full-period LCG

Theorem (Hull–Dobell). The LCG $z_i = (a z_{i-1} + c) \bmod m$ has full period $m$ if and only if

1. $c$ and $m$ are relatively prime (their only common divisor is $1$);
2. every prime $q$ that divides $m$ also divides $a - 1$;
3. if $4$ divides $m$, then $4$ divides $a - 1$.

Why three conditions. A full-period LCG must visit every one of $0, 1, \dots, m-1$. The test is to look at the sequence “through a coarser lens”: for a divisor $d$ of $m$, watch only $z_i \bmod d$. If the coarse sequence misses some residue mod $d$, the full sequence misses some value mod $m$. Condition 1 removes the obstruction coming from the increment $c$; condition 2 removes the obstruction coming from the multiplier $a$, one prime $q \mid m$ at a time; condition 3 is needed because the prime $2$ misbehaves in a way no other prime does, so the lens $d = 4$ must be checked separately.

Why each condition is necessary.

• Condition 1. Suppose $c$ and $m$ share a divisor $d > 1$. Modulo $d$ the recursion becomes $z_i \equiv a z_{i-1} \pmod d$, so if $z_0 \equiv 0 \pmod d$ then every $z_i \equiv 0 \pmod d$: the sequence never leaves the multiples of $d$. In the example above $c = 0$, so $\gcd(c, 64) = 64$ and full period is impossible for any seed.
• Condition 2. Let $q$ be a prime dividing $m$. Modulo $q$ the map is $x \mapsto ax + c$. If $a \not\equiv 1 \pmod q$, then $1 - a$ is invertible mod $q$ and the map has a fixed point $x^* \equiv c\,(1-a)^{-1} \pmod q$. Any seed with $z_0 \equiv x^* \pmod q$ stays in that residue class forever, so the period is less than $m$. Hence we need $a \equiv 1 \pmod q$, i.e. $q \mid a - 1$; then the map mod $q$ is $x \mapsto x + c$, which does cycle through all $q$ residues because $\gcd(c, q) = 1$ by condition 1.
• Condition 3. Condition 2 with $q = 2$ only gives $a$ odd, i.e. $a \equiv 1$ or $3 \pmod 4$. If $a \equiv 3 \pmod 4$, then applying the map twice mod $4$ gives $x \mapsto 9x + 3c + c \equiv x + 4c \equiv x \pmod 4$, so the sequence mod $4$ has period at most $2$ and cannot visit all four residues. (Concretely, $a = 3$, $c = 1$, $m = 4$, $z_0 = 0$ gives $0, 1, 0, 1, \dots$) So $4 \mid m$ forces $a \equiv 1 \pmod 4$.

Sufficiency — that the three conditions together guarantee period exactly $m$ — is the harder direction and is what Hull and Dobell proved; it is not needed for the course beyond the statement.

Check on the earlier example ($m = 16$, $a = 5$, $c = 3$): $\gcd(3, 16) = 1$; the only prime dividing $16$ is $2$, and $2 \mid 4 = a - 1$; $4 \mid 16$ and $4 \mid 4$. All three hold, so full period $16$, as observed.

3.2.3. Two types of LCGs

Type 1: Mixed generators (LCGs with $c > 0$)

We want $m$ large (observation 1). The standard choice is

$$m = 2^b,$$

where $b$ is the number of bits available in a machine word. On a 32-bit system one bit is used for the sign, so $b = 31$ and $m = 2^{31}$. The reason for a power of two is speed: computing $x \bmod 2^b$ costs nothing, because it is just keeping the lowest $b$ bits of $x$ (integer overflow does it automatically).

With $m = 2^b$ the only prime dividing $m$ is $2$, so the three Hull–Dobell conditions specialize as follows:

1. $\gcd(c, 2^b) = 1$ $\;\Longleftrightarrow\;$ $c$ is odd.
2. $2 \mid a - 1$ $\;\Longleftrightarrow\;$ $a$ is odd.
3. $4 \mid 2^b$, hence $4 \mid a - 1$ $\;\Longleftrightarrow\;$ $a \equiv 1 \pmod 4$.

Item 3 implies item 2, so item 2 is redundant, but it is a common mistake to stop there: $a = 3, 7, 11, \dots$ are odd yet give short periods. The rule for $m = 2^b$ is therefore simply

$$c \text{ odd} \quad\text{and}\quad a \equiv 1 \pmod 4 .$$

A good mixed LCG: $m = 2^{31}$, $a = 314159269$, $c = 453806245$.
Check: $c$ is odd, and $a - 1 = 314159268 = 4 \cdot 78539817$, so $a \equiv 1 \pmod 4$. Full period $2^{31} \approx 2.1 \times 10^9$.

Type 2: Multiplicative generators (LCGs with $c = 0$)

$$z_i = a\, z_{i-1} \bmod m, \qquad\text{so}\qquad z_i = a^i z_0 \bmod m .$$

1. A multiplicative LCG can no longer have full period. Condition 1 of Hull–Dobell fails, since $\gcd(c, m) = \gcd(0, m) = m \ne 1$. Concretely: $0$ is a fixed point ($a \cdot 0 = 0$), so the sequence can never pass through $0$ without staying there forever. A useful generator must therefore avoid $0$ altogether, and the best it can do is cycle through $\{1, 2, \dots, m-1\}$: the maximum possible period is $m - 1$.
2. Period $m - 1$ is achievable if $a$ and $m$ are chosen carefully.
3. Choose $m$ prime, specifically the largest prime less than $2^b$. (For $b = 31$ this is $2^{31} - 1 = 2147483647$, which happens to be prime.) Primality matters because of the next point.
4. When is the period exactly $m - 1$? Suppose $m$ is prime and $z_0 \ne 0$. The sequence returns to its start at step $l$ when

   $$a^l z_0 \equiv z_0 \pmod m .$$

   Because $m$ is prime and $z_0 \not\equiv 0$, we can cancel $z_0$, leaving

   $$a^l \equiv 1 \pmod m, \quad\text{i.e.}\quad m \mid a^l - 1 .$$

   So the period is the smallest $l \ge 1$ with $a^l \equiv 1 \pmod m$, and it does not depend on the seed. By Fermat’s little theorem $a^{m-1} \equiv 1 \pmod m$ always, so this smallest $l$ divides $m - 1$. The period equals $m - 1$ exactly when the smallest such $l$ is $m - 1$:

   $$\min\{\, l \ge 1 : a^l - 1 \equiv 0 \pmod m \,\} = m - 1 .$$

   An $a$ with this property is called a primitive element modulo $m$ (or primitive root). Hence: if $m$ is prime and $a$ is a primitive element modulo $m$, the multiplicative LCG has period $m - 1$ for every seed $z_0 \ne 0$.

Such generators — multiplicative ($c = 0$) with a prime modulus $m$ and a primitive multiplier $a$ — are called prime modulus multiplicative LCGs (PMMLCGs). The name lists the ingredients: prime modulus is point 3, multiplicative is $c = 0$, and the period $m - 1$ comes from point 4.

Example. $m = 7$ (prime).

• $a = 3$: $3^1 = 3,\ 3^2 = 2,\ 3^3 = 6,\ 3^4 = 4,\ 3^5 = 5,\ 3^6 = 1 \pmod 7$. The first $l$ with $3^l \equiv 1$ is $l = 6 = m - 1$, so $3$ is primitive and the generator $z_i = 3 z_{i-1} \bmod 7$ has period $6$: from $z_0 = 1$ it runs $1, 3, 2, 6, 4, 5, 1, \dots$
• $a = 2$: $2^1 = 2,\ 2^2 = 4,\ 2^3 = 1 \pmod 7$. Here $l = 3 \ne 6$, so $2$ is not primitive; the generator $z_i = 2 z_{i-1} \bmod 7$ splits $\{1,\dots,6\}$ into two cycles of length $3$: $1, 2, 4$ and $3, 6, 5$.

A good PMMLCG. $m = 2^{31} - 1 = 2147483647$ (prime), $a = 7^5 = 16807$ or $a = 630360016$. Both multipliers are primitive elements modulo $m$, so the period is $m - 1 = 2^{31} - 2 \approx 2.1 \times 10^9$.

A bad multiplicative LCG: RANDU. $m = 2^{31}$, $a = 2^{16} + 3 = 65539$, $c = 0$. It was the default generator on IBM mainframes in the 1960s–70s, and many early simulation papers were based on it. Its period is long enough ($2^{29}$, the maximum $m/4$ for a power-of-two modulus), but it has a severe defect: because $a = 2^{16} + 3$,

$$a^2 = 2^{32} + 6 \cdot 2^{16} + 9 \equiv 6a - 9 \pmod{2^{31}},$$

and multiplying by $z_i$ gives the exact linear relation

$$z_{i+2} \equiv 6\, z_{i+1} - 9\, z_i \pmod{2^{31}} .$$

So every value is determined by the two before it via a fixed formula with tiny coefficients. Geometrically, all consecutive triples $(u_i, u_{i+1}, u_{i+2})$ lie on just $15$ parallel planes in the unit cube, instead of filling it — a gross violation of independence (property 2 of a good RNG). Any simulation that uses three successive numbers together (e.g. to generate a point in 3-D) inherits this structure.

The lesson of RANDU: a long period is necessary but not sufficient. Type 1 and Type 2 generators both need to be tested for uniformity and independence, which is the subject of the next section.

3.3. Other Kinds of Generators

More general congruential generators. Everything so far has the shape

$$z_i = g(z_{i-1}, z_{i-2}, \dots, z_0) \bmod m, \qquad u_i = \frac{z_i}{m},$$

for some function $g$ of the previous values. An LCG is the special case $g(z_{i-1}) = a z_{i-1} + c$. Other choices of $g$ give other generators.

Multiple recursive generators (MRGs). Take $g$ linear in the last $q$ values:

$$z_i = (a_1 z_{i-1} + a_2 z_{i-2} + \cdots + a_q z_{i-q}) \bmod m,$$

where $a_1, \dots, a_q$ are integer constants. This is the multiplicative LCG with $q$ lags instead of one ($q = 1$ gives back $z_i = a_1 z_{i-1} \bmod m$).

• Advantage: a huge period. The state of the generator is now the $q$-tuple $(z_{i-1}, \dots, z_{i-q})$, of which there are $m^q$ possibilities. The all-zero tuple is a fixed point (as $0$ was for the multiplicative LCG), so the maximum period is $m^q - 1$, and this is attained for suitable $a_1, \dots, a_q$ when $m$ is prime. With $m \approx 2^{31}$ and $q = 3$ this is about $2^{93} \approx 10^{28}$.
• Disadvantage: $q$ multipliers to choose instead of one, and the theory for choosing them well is harder. The output sequence still looks like an LCG’s — a linear recurrence with lattice structure — just in $q$ dimensions rather than one.

Other choices of $g$.

• Quadratic congruential generator: $g(z_{i-1}) = a z_{i-1}^2 + b z_{i-1} + c$, i.e.

  $$z_i = (a z_{i-1}^2 + b z_{i-1} + c) \bmod m .$$

  This is nonlinear in $z_{i-1}$, but it still depends on the previous value only, so the period is still at most $m$.

• Fibonacci generator: $g(z_{i-1}, z_{i-2}) = z_{i-1} + z_{i-2}$, i.e.

  $$z_i = (z_{i-1} + z_{i-2}) \bmod m .$$

  This is the MRG with $q = 2$ and $a_1 = a_2 = 1$. It is very fast (no multiplication) but statistically poor. To see why, look at three consecutive outputs $u_{i-2}, u_{i-1}, u_i$ with $u_i = (u_{i-1} + u_{i-2}) \bmod 1$. If the sum does not wrap past $1$, then $u_i$ is larger than both $u_{i-1}$ and $u_{i-2}$; if it does wrap, $u_i = u_{i-1} + u_{i-2} - 1$ is smaller than both. So $u_i$ is never the middle value of the three. For i.i.d. uniforms each of the six orderings of three values should occur with probability $1/6$; the Fibonacci generator produces two of them with probability $0$.

Practical issues

1. Period. The period of the RNG should be at least an order of magnitude larger than the total number of random numbers the simulation will consume. Otherwise the sequence wraps around and later parts of the run reuse the same numbers as earlier parts, destroying independence.
2. Streams. A stream is a segment of the RNG’s cycle starting from a chosen seed. Use a different stream for each replication of a simulation, with the seeds spaced far enough apart that the streams do not overlap. Non-overlapping streams from one long-period generator then behave like independent generators.

flowchart LR
    s1["stream 1<br/>seed z₀⁽¹⁾"] --> s2["stream 2<br/>seed z₀⁽²⁾"] --> s3["stream 3<br/>seed z₀⁽³⁾"] --> d[". . ."] --> s1

3. Separate streams for separate inputs. Use different streams to simulate quantities that should be independent — e.g. one stream for interarrival times and another for service times. Besides keeping the inputs independent, this synchronizes the random numbers: when two alternative systems are simulated, each can be fed the same interarrival stream and the same service stream, so any difference in the output is due to the systems and not to the random inputs (the “sharper comparison” mentioned at the start of the chapter).

3.4. Testing RNGs

Two kinds of tests:

• Theoretical tests (number theory). These analyze the generator’s parameters $(a, c, m)$ mathematically and say something about the entire cycle — a global property — without generating a single $u_i$. Example: the spectral test, which measures how far apart the parallel hyperplanes are on which all $k$-tuples $(u_i, \dots, u_{i+k-1})$ of an LCG lie (recall RANDU’s 15 planes).
• Empirical tests (statistics). Generate a finite sample $u_1, \dots, u_n$ and test the null hypothesis

  $$H_0:\ u_1, u_2, \dots \text{ are i.i.d. } U(0,1).$$

  These are only local tests: they examine the part of the cycle that was actually generated, so a defect elsewhere in the cycle goes undetected. $H_0$ has two parts, hence two families of tests: tests for uniformity and tests for independence.

3.4.1. Test for uniformity (frequency test / $\chi^2$ test)

Goal. Check whether the $u_i$ appear to be uniformly distributed over $[0, 1]$.

Steps.

1. Divide $[0, 1]$ into $k$ subintervals of equal length $1/k$.
2. Generate $n$ random numbers $u_1, \dots, u_n$ from the RNG. Rule of thumb: $k \ge 100$ and $n/k \ge 5$ (at least 5 expected observations per interval, so that the approximation in step 4 is accurate).
3. Count the number of $u_i$ falling in the $j$-th interval; call it $f_j$, so $f_1 + \cdots + f_k = n$.

4. Compute the test statistic

   $$\chi^2 = \sum_{j=1}^{k} \frac{\left(f_j - \dfrac{n}{k}\right)^2}{\dfrac{n}{k}} .$$

   Under $H_0$ each interval should receive about $n/k$ points, so $\chi^2$ measures the total squared deviation of the observed counts from that. When $n$ is large, $\chi^2$ approximately follows a chi-square distribution with $k-1$ degrees of freedom.
5. Reject $H_0$ at significance level $\alpha$ if

   $$\chi^2 > \chi^2_{k-1,\,\alpha},$$

   where the threshold $\chi^2_{k-1,\,\alpha}$ is defined by $P\!\left(\chi^2_{k-1} > \chi^2_{k-1,\,\alpha}\right) = \alpha$ (upper $\alpha$-quantile, as with $z_\alpha$ and $t_{\alpha, n-1}$). Large $\chi^2$ means the counts are too uneven to be uniform; only the upper tail is used.

What is the chi-square distribution? If $X_1, \dots, X_d$ are i.i.d. $N(0,1)$, then

$$Q = \sum_{i=1}^{d} X_i^2$$

is a chi-square random variable with $d$ degrees of freedom, written $Q \sim \chi^2_d$.

Why $\chi^2$ is approximately $\chi^2_{k-1}$ under $H_0$.

Assume $H_0$ is true. Each $u_i$ lands in interval $j$ with probability $1/k$, independently of the others, so the count $f_j$ is binomial:

$$f_j \sim \text{Bin}\!\left(n, \tfrac{1}{k}\right), \qquad P(f_j = r) = \binom{n}{r}\left(\tfrac{1}{k}\right)^{r}\left(1 - \tfrac{1}{k}\right)^{n-r},$$

$$\mathbf{E}[f_j] = \frac{n}{k}, \qquad \operatorname{Var}(f_j) = \frac{n}{k}\left(1 - \frac{1}{k}\right).$$

By the CLT (a binomial is a sum of $n$ Bernoullis), the standardized count is approximately standard normal for large $n$:

$$Z_j = \frac{f_j - n/k}{\sqrt{\dfrac{n}{k}\left(1 - \dfrac{1}{k}\right)}} \;\approx\; N(0,1).$$

If the $Z_j$ were independent, $\sum_{j=1}^k Z_j^2$ would be $\chi^2_k$ by the definition above. Writing this sum out,

$$\sum_{j=1}^{k} Z_j^2 = \frac{k}{k-1}\sum_{j=1}^{k}\frac{(f_j - n/k)^2}{n/k} = \frac{k}{k-1}\,\chi^2, \qquad\text{i.e.}\qquad \chi^2 = \frac{k-1}{k}\sum_{j=1}^{k} Z_j^2 .$$

But the $Z_j$ are not independent: the counts satisfy $f_1 + \cdots + f_k = n$, one linear constraint, so only $k-1$ of them are free. The precise result (Pearson) is that this constraint costs exactly one degree of freedom and the factor $(k-1)/k$ is exactly what compensates: $\chi^2 \approx \chi^2_{k-1}$, not $\chi^2_k$. The heuristic to remember is: $k$ approximately normal terms, one constraint, $k-1$ degrees of freedom.

The threshold $\chi^2_{k-1,\alpha}$ cuts off an upper tail of area $\alpha$. Under $H_0$ the statistic lands in that tail with probability only $\alpha$, so landing there is taken as evidence against uniformity.

Example. Twenty numbers are generated from an RNG:

$$0.10,\ 0.32,\ 0.76,\ 0.54,\ 0.08,\ 0.23,\ 0.67,\ 0.89,\ 0.41,\ 0.15,$$

$$0.95,\ 0.60,\ 0.28,\ 0.73,\ 0.49,\ 0.05,\ 0.37,\ 0.82,\ 0.00,\ 0.91 .$$

Are these from a uniform distribution on $[0,1]$, at significance level $\alpha = 0.05$?

$$H_0:\ \text{the numbers are i.i.d. } U(0,1).$$

Step 1. Divide $[0,1]$ into $k$ subintervals. With only $n = 20$ numbers we cannot meet $k \ge 100$; take $k = 4$, which at least satisfies $n/k = 5$. (This is a hand-sized illustration — a real test would use thousands of numbers.)

Step 2. $n = 20$ numbers are already generated.

Step 3. Count how many fall in each subinterval. Expected count per interval: $n/k = 20/4 = 5$.

$j$	interval	numbers in it	$f_j$	$n/k$	$f_j - n/k$
1	$[0, 0.25)$	0.10, 0.08, 0.23, 0.15, 0.05, 0.00	6	5	$+1$
2	$[0.25, 0.5)$	0.32, 0.41, 0.28, 0.49, 0.37	5	5	$0$
3	$[0.5, 0.75)$	0.54, 0.67, 0.60, 0.73	4	5	$-1$
4	$[0.75, 1)$	0.76, 0.89, 0.95, 0.82, 0.91	5	5	$0$
		total	20	20	$0$

Step 4. Compute the test statistic:

$$\chi^2 = \sum_{j=1}^{4} \frac{(f_j - n/k)^2}{n/k} = \frac{(+1)^2 + 0^2 + (-1)^2 + 0^2}{5} = \frac{2}{5} = 0.4 .$$

Step 5. Compare with the threshold. Degrees of freedom $k - 1 = 3$, and from the chi-square table

$$\chi^2_{3,\,0.05} = 7.815 .$$

Since $\chi^2 = 0.4 \le 7.815$, we do not reject $H_0$ at the $5\%$ level: there is no evidence that these numbers are non-uniform.

Remark. “Do not reject” is not “proved uniform.” First, the test only looks at the marginal distribution — a sequence like $0.1, 0.2, \dots, 0.9, 0.1, 0.2, \dots$ would pass the uniformity test easily while being obviously non-random, which is why a test for independence is also needed (§3.4.2). Second, with $n = 20$ the test has very little power; the counts would have to be wildly off to produce a $\chi^2$ above $7.8$.

Serial tests (higher-dimensional $\chi^2$ tests)

The one-dimensional test only checks the marginal distribution. The serial test checks a consequence of the full hypothesis:

If the $u_i$ are i.i.d. $U(0,1)$, then the non-overlapping $d$-tuples

$$(u_1, \dots, u_d),\quad (u_{d+1}, \dots, u_{2d}),\quad \dots$$

are uniformly distributed over the hypercube $[0,1]^d$.

This is an implication, not an equivalence: $A \Rightarrow B$, so failing $B$ disproves $A$. That is enough for a test — if the tuples are not uniform on the cube, the numbers are not i.i.d. $U(0,1)$.

Case $d = 2$.

1. Divide the square $[0,1]^2$ into $k^2$ subsquares of equal size.
2. Generate $n$ pairs $U_1 = (u_1, u_2)$, $U_2 = (u_3, u_4)$, …, $U_n = (u_{2n-1}, u_{2n})$. Note each $u$ is used once: $2n$ numbers give $n$ pairs.
3. Let $f_{ij}$ be the number of $U$’s falling in subsquare $(i,j)$. The expected count per square is $n/k^2$, so

   $$\chi^2 = \sum_{i=1}^{k}\sum_{j=1}^{k} \frac{\left(f_{ij} - \dfrac{n}{k^2}\right)^2}{\dfrac{n}{k^2}} \;\approx\; \chi^2_{k^2 - 1}$$

   ($k^2$ cells, one constraint $\sum_{i,j} f_{ij} = n$, hence $k^2 - 1$ degrees of freedom — same counting as in §3.4.1).
4. Reject $H_0$ at level $\alpha$ if $\chi^2 > \chi^2_{k^2-1,\,\alpha}$; see [[Review of Prob and Stats#2.6.3. Statistical tables ($z$, $t$, $ chi 2$)]].

Why this is an indirect test of independence. The one-dimensional test cannot see dependence at all: a sequence can have a perfectly uniform histogram and still be highly predictable. The serial test looks at successive values jointly, so such structure shows up as clumping in the square.

Limitations.

• Indirect. Passing the serial test does not prove independence; it only means no dependence was detected at the resolution used. RANDU, for instance, passes comfortably in $d = 2$ and fails only in $d = 3$, where its 15 planes appear.
• Costly for large $d$. The number of cells is $k^d$, and the rule $n/k^d \ge 5$ forces $n \ge 5k^d$ — exponential in $d$. With $k = 10$ and $d = 5$ that is already $500{,}000$ tuples, i.e. $2.5$ million random numbers. In practice one uses a small $k$ for larger $d$, or switches to a theoretical test such as the spectral test, which examines all dimensions at once without sampling.

3.4.2. Test for independence (runs test)

What “runs” means. A run up is an unbroken subsequence of $u_i$’s that is monotonically increasing; the run ends as soon as the next number is smaller. The sequence is thus chopped into maximal ascending pieces, and the test looks at how long those pieces are. The name is literal: how far does the sequence “run” upward before turning around?

The idea is that the lengths of these runs have a known distribution when the $u_i$ are i.i.d., and this distribution depends on the order of the numbers, not just on their values. That is exactly what the frequency test cannot see.

Example (the same 20 numbers as in §3.4.1). Mark $\nearrow$ where $u_{i+1} > u_i$ and $|$ where a run ends:

$$0.10 \nearrow 0.32 \nearrow 0.76 \;\big|\; 0.54 \;\big|\; 0.08 \nearrow 0.23 \nearrow 0.67 \nearrow 0.89 \;\big|\; 0.41 \;\big|\; 0.15 \nearrow 0.95 \;\big|\;$$

$$0.60 \;\big|\; 0.28 \nearrow 0.73 \;\big|\; 0.49 \;\big|\; 0.05 \nearrow 0.37 \nearrow 0.82 \;\big|\; 0.00 \nearrow 0.91$$

The run lengths are $3, 1, 4, 1, 2, 1, 2, 1, 3, 2$ (they sum to $n = 20$, as they must).

The counts. Let

$$r_i = \begin{cases} \text{number of runs up of length } i, & i = 1, 2, \dots, 5, \\ \text{number of runs up of length} \ge 6, & i = 6. \end{cases}$$

Lengths $\ge 6$ are pooled into one category because they are rare, so their individual counts would be too small for the approximation to hold — the same reason the frequency test wants $n/k \ge 5$ per cell.

For the example: $r_1 = 4$, $r_2 = 3$, $r_3 = 2$, $r_4 = 1$, $r_5 = r_6 = 0$.

Expected counts. Under $H_0$ the expected number of runs of each length is $\mathbf{E}[r_i] \approx n\,b_i$, where the $b_i$ are universal constants (they do not depend on the distribution, only on the ordering being random):

$$b = \left(\tfrac{1}{6},\; \tfrac{5}{24},\; \tfrac{11}{120},\; \tfrac{19}{720},\; \tfrac{29}{5040},\; \tfrac{1}{840}\right) \approx (0.1667,\; 0.2083,\; 0.0917,\; 0.0264,\; 0.0058,\; 0.0012).$$

Note $\sum_i b_i = 1/2$: a random sequence of length $n$ splits into about $n/2$ runs up. Runs of length $2$ are the most common, not length $1$.

The test statistic.

$$R = \frac{1}{n}\sum_{i=1}^{6}\sum_{j=1}^{6} a_{ij}\,(r_i - n b_i)(r_j - n b_j),$$

where the $a_{ij}$ are constants tabulated in the textbook (Law & Kelton, p. 408); the matrix is symmetric, with entries in the thousands — the first row begins $4529.4,\ 9044.9,\ 13568, \dots$

Why the double sum, unlike the $\chi^2$ statistic of §3.4.1? Because the $r_i$ are strongly correlated with one another: a long run uses up numbers that then cannot start short runs, so an excess of $r_4$ forces a deficit elsewhere. A plain sum of squared deviations $\sum_i (r_i - nb_i)^2 / (nb_i)$ would therefore not have a chi-square distribution. The matrix $(a_{ij})$ is (up to scaling) the inverse of the covariance matrix of the vector $(r_1, \dots, r_6)$, and multiplying by it removes those correlations — this is the general recipe for turning a correlated normal vector into a chi-square variable.

Decision rule. For large $n$ (the textbook recommends $n \ge 4000$), under

$$H_0:\ \text{the } u_i \text{ are i.i.d.}$$

$R$ has approximately a chi-square distribution with 6 degrees of freedom — not $6 - 1 = 5$, because here no constraint like $\sum f_j = n$ ties the six counts together. Reject $H_0$ at significance level $\alpha$ if

$$R > \chi^2_{6,\,\alpha},$$

with the threshold read from [[Review of Prob and Stats#2.6.3. Statistical tables ($z$, $t$, $ chi 2$)]]: e.g. $\chi^2_{6,\,0.05} = 12.592$.

(The 20-number example cannot be carried further: with $n = 20 \ll 4000$ the chi-square approximation is meaningless. It only illustrates how the $r_i$ are counted.)

In practice, perform the runs test first. Two reasons:

1. The frequency test’s null distribution was derived assuming independence — that is what made $f_j$ binomial in §3.4.1. If the numbers are dependent, the $\chi^2$ threshold is not valid, so testing uniformity first can give a meaningless answer.
2. Dependence is the more damaging defect. A generator whose output is uniform but correlated will quietly bias any simulation, whereas mild non-uniformity is often visible and correctable.