🎲Randomness and Estimation

Random Numbers and Monte Carlo Simulation

Pseudo-random generators and their flaws, sampling distributions, numpy.random.Generator, confidence intervals

Random numbers fall mainly in the planes, 1968

In the 1960s, one of the most widely distributed random-number generators in the world was RANDU, part of IBM's Scientific Subroutine Package for the System/360 family of computers. It was simple and fast: multiply the previous integer by 65539 and keep the remainder after dividing by

2^{31}

. Its numbers passed the usual statistical checks: they were spread evenly between 0 and 1, with the right mean and variance. Thousands of simulations — of nuclear reactors, of queues, of physical systems — were built on it. In 1968 George Marsaglia of Boeing Scientific Research Laboratories published a short paper with a memorable title: Random numbers fall mainly in the planes. He proved that the successive outputs of every generator of RANDU's type, taken in groups of

n

, lie on a small number of parallel hyperplanes. For RANDU the result was spectacular: every triple of consecutive numbers lies on one of just 15 planes in the unit cube. Any simulation that used three consecutive numbers as a point in space — a direction, a position, three random loads — had been sampling only those planes. Results computed with RANDU went on appearing in the scientific literature for years afterwards.

Monte Carlo methods (lesson 8-5) need streams of numbers that behave as if independent and uniformly distributed. A computer is deterministic, so it produces pseudo-random numbers: a recurrence

x_{k+1} = f(x_k)

whose outputs pass statistical tests of randomness. Reproducibility is then a feature: the same seed gives the same stream, so a simulation can be repeated exactly — which is why every lab in this course sets one.

💡Linear congruential generators and their lattice

The classic generator is the linear congruential generator

x_{k+1} = (a\,x_k + c) \bmod m

u_k = x_k/m

. For RANDU,

a = 65539 = 2^{16} + 3

c = 0

m = 2^{31}

. Squaring the multiplier modulo

2^{31}

gives

a^2 \equiv 6a - 9

, so

x_{k+2} = 6x_{k+1} - 9x_k \pmod{2^{31}} \quad\Longrightarrow\quad 9u_k - 6u_{k+1} + u_{k+2} \in \mathbb Z.

Every triple satisfies one of 15 integer equations — 15 planes. Its one-dimensional statistics are fine (mean 0.503, variance 0.0832 against 0.5 and 0.0833); only the joint distribution reveals the flaw. Modern generators — the Mersenne Twister (1998), PCG (2014, NumPy's default since 2019), Philox and xoshiro — are designed and tested against such structure in hundreds of dimensions.

Derivation

RANDU's planes in three lines

Square the multiplier exactly:

a^2 = (2^{16} + 3)^2 = 2^{32} + 6 \cdot 2^{16} + 9 = 2^{32} + 6(2^{16} + 3) - 9 = 2^{32} + 6a - 9.

Modulo

2^{31}

the term

2^{32} = 2 \cdot 2^{31}

vanishes, so

a^2 \equiv 6a - 9

. Two steps of the generator multiply by

a^2

, hence

x_{k+2} \equiv a^2x_k \equiv 6\,a x_k - 9x_k \equiv 6x_{k+1} - 9x_k \pmod{2^{31}}.

Divide by

2^{31}

u_{k+2} - 6u_{k+1} + 9u_k

is a whole number. With every

u

between 0 and 1 that whole number lies between

-6

and

10

, so only a handful of values — 15 in practice — are possible, and each defines a plane

9x - 6y + z = \text{integer}

in the unit cube. A single unlucky choice of multiplier, made for speed, put every triple of a world-standard generator on 15 sheets.

Any one who considers arithmetical methods of producing random digits is, of course, in a state of sin.

– John von Neumann, Various Techniques Used in Connection with Random Digits (1951)

Sampling other distributions. Uniform numbers are turned into whatever the model needs:

Inverse transform: if $U$ is uniform and $F$ is a cumulative distribution, $X = F^{-1}(U)$ has distribution $F$ . Wind speeds with a Weibull distribution, $F(v) = 1 - e^{-(v/c)^k}$ , come from $v = c(-\ln(1-U))^{1/k}$ .
Box–Muller (1958): two uniforms give two independent standard normals, $\sqrt{-2\ln U_1}\cos 2\pi U_2$ and $\sqrt{-2\ln U_1}\sin 2\pi U_2$ .
In practice, use numpy.random.default_rng(seed) and its methods (normal, weibull, gumbel, multivariate_normal, …), which implement fast, tested algorithms.

Check samples against the target distribution with a goodness-of-fit test such as Kolmogorov–Smirnov (scipy.stats.kstest).

Derivation

Why inverse transform and Box–Muller work

Inverse transform. Let

U

be uniform on

(0, 1)

, so

P(U \le u) = u

, and let

F

be increasing. Then

P\big(F^{-1}(U) \le x\big) = P\big(U \le F(x)\big) = F(x),

which is exactly the statement that

F^{-1}(U)

has distribution

F

. For the Weibull, solve

u = 1 - e^{-(v/c)^k}

for

v

e^{-(v/c)^k} = 1 - u

(v/c)^k = -\ln(1 - u)

, so

v = c\,(-\ln(1 - u))^{1/k}

Box–Muller. Two independent standard normals

Z_1, Z_2

have joint density

\frac1{2\pi}e^{-(z_1^2 + z_2^2)/2}

, which depends only on the radius

r

. So their angle is uniform on

[0, 2\pi)

— that is

2\pi U_2

. For the radius, the probability of landing outside radius

\sqrt s

is a one-line polar integral:

P(R^2 > s) = \int_{\sqrt s}^{\infty}\frac{1}{2\pi}e^{-r^2/2}\,2\pi r\,dr = e^{-s/2}.

Inverting that by the rule above (with

1 - U

replaced by

U_1

, which is just as uniform) gives

R^2 = -2\ln U_1

. Point

(R\cos 2\pi U_2,\ R\sin 2\pi U_2)

is then a pair of independent normals.

📐Wind power is not average wind cubed

Wind speed at a site follows a Weibull distribution with

k = 2

c = 8

m/s (mean 7.09 m/s). Turbine power is proportional to

v^3

. From 100,000 sampled speeds, the mean of

v^3

is 674.5 while the cube of the mean speed is 353.4: designing with the average wind would predict only 52% of the true average power. Nonlinear models of uncertain inputs must be averaged over the distribution, not evaluated at the average — the most common error in uncertainty analysis.

Derivation

The 52%, exactly — and the confidence interval for a probability

The wind. The Weibull distribution has moments

E[v^n] = c^n\,\Gamma(1 + n/k)

(a standard integral). With

k = 2

c = 8

E[v] = 8\,\Gamma(1.5) = 7.09 \text{ m/s}, \qquad E[v^3] = 512\,\Gamma(2.5) = 680.6, \qquad \frac{(E[v])^3}{E[v^3]} = \frac{356.4}{680.6} = 0.52.

Because

v^3

curves upwards, the strong gusts add far more power than the calm spells remove: the mean of the cube exceeds the cube of the mean (Jensen's inequality). The 100,000-sample estimates, 674.5 and 353.4, give the same 52%.

The interval. Counting failures among

N

independent trials, each failing with probability

p

, gives a binomial count

K

with variance

Np(1 - p)

. The estimate

\hat p = K/N

therefore has variance

p(1 - p)/N

(dividing by

N

divides a variance by

N^2

). For large

N

the count is nearly normal, and 95% of a normal distribution lies within 1.96 standard deviations — which gives the interval in the box above.

💡Confidence intervals

A Monte Carlo estimate is itself random. For a probability

p

estimated from

N

samples, the standard error is

\sqrt{p(1-p)/N}

, and an approximate 95% confidence interval is

\hat p \pm 1.96\sqrt{\hat p(1-\hat p)/N}

. For a steel tie with normally distributed resistance and a Gumbel-distributed extreme load,

N = 10^4, 10^5, 10^6

give

p_f = 0.0074 \pm 0.0017

0.0070 \pm 0.0005

0.0076 \pm 0.0002

; numerical integration of

P(R < s)f_S(s)

gives 0.00742. Always report the interval: an estimate of

10^{-4}

from 2 failures in 20,000 samples is barely an estimate at all. Rare events need importance sampling or structural-reliability methods (FORM/SORM) rather than brute force.

💡Where you will use this

Structural and Geotechnical reliability: failure probabilities with uncertain loads, strengths and soil properties. Rocket Propulsion and Guidance: Monte Carlo dispersion analysis — thousands of simulated flights with sampled thrust, mass and wind. Orbital Mechanics: collision probabilities. GNSS and Navigation: simulated sensor noise to test filters. Every course: whenever a number is uncertain, sampling is the honest way to see what follows.

📚Further reading

Donald Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms (3rd ed., Addison-Wesley, 1997), chapter 3, 'Random Numbers' — linear congruential generators, the spectral test and transformation methods.
George Marsaglia, 'Random numbers fall mainly in the planes', Proceedings of the National Academy of Sciences 61(1), 25–28 (1968).
Melissa O'Neill, PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation, Harvey Mudd College technical report HMC-CS-2014-0905 (2014) — the generator behind numpy.random.default_rng.

↓But there is a problem

Monte Carlo shows how uncertainty in the inputs becomes uncertainty in the answer — by brute force, thousands of runs at a time. Often you need the answer in one calculation, every second, on board. On 10 February 2009 a screening service that checks the whole public satellite catalogue for close approaches predicted that two satellites over Siberia would pass 584 m apart. It was not considered a dangerous number. The prediction was missing one thing: how uncertain it was.