The law of large numbers

The opening note showed the share of heads in coin flips drifting toward one half and called a probability a long-run frequency. The law of large numbers turns that observation into a theorem: the average of independent copies of a random variable converges to its expectation. It is what makes a probability measurable, a validation score meaningful and a simulation a method of computation. This note states the law in its two forms, explains the two kinds of convergence they use, and turns the law into Monte Carlo estimation.

Notation. $X_1, X_2, \dots$ are i.i.d. with mean $\mu$ (and, where needed, variance $\sigma^2$); $\bar X_n = (X_1 + \dots + X_n)/n$. An event happens almost surely if it has probability 1.

The weak law#

For every $\varepsilon > 0$,

$$P\big(\lvert \bar X_n - \mu \rvert \ge \varepsilon\big) \to 0 \quad \text{as } n \to \infty .$$

The average converges in probability to $\mu$: for large $n$, a deviation of any fixed size is unlikely. With a finite variance the proof is one line of the previous note: Chebyshev’s inequality bounds the probability by $\sigma^2/(n\varepsilon^2)$. Khinchin showed that a finite mean is enough, with no assumption on the variance.

The strong law#

The weak law speaks about one large $n$ at a time. It leaves open whether a particular infinite sequence of flips could keep making excursions away from the mean forever, rarer and rarer but never stopping. The strong law rules this out:

$$P\big(\bar X_n \to \mu\big) = 1 .$$

The average converges almost surely: with probability 1 the sequence of averages, as a sequence of numbers, has the limit $\mu$. Kolmogorov proved it under the same assumption as Khinchin, a finite mean. Almost sure convergence implies convergence in probability, not the other way round.

For practice the distinction rarely matters: both laws say that one long simulation, or one large validation set, gives an average close to the expectation. It matters for algorithms that run a single trajectory forever, such as stochastic gradient descent or a Markov chain, where the question is whether that one trajectory converges.

The Borel–Cantelli lemmas#

Optional depth: the tool that separates the two kinds of convergence, and the idea of the simplest proof of the strong law. The rest of the note does not depend on it.

The difference is easiest to see on indicators. Let $A_1, A_2, \dots$ be independent events and $X_n = \mathbb 1_{A_n}$. If $P(A_n) \to 0$, then $X_n \to 0$ in probability. Whether $X_n \to 0$ almost surely, that is, whether the events eventually stop happening, is decided by the Borel–Cantelli lemmas: if $\sum_n P(A_n) < \infty$, only finitely many $A_n$ happen, almost surely; if the events are independent and $\sum_n P(A_n) = \infty$, infinitely many happen, almost surely.

Figure

Convergence in probability and almost surely

With $P(A_n) = 1/n$ the expected number of events up to $10^5$ is $1 + \frac12 + \dots + \frac{1}{10^5} \approx 12.1$, and they never stop: the chance of at least one between $10^4$ and $10^5$ is $1 - \frac{10^4}{10^5} = 0.9$. With $P(A_n) = 1/n^2$ the expected total is $\pi^2/6 \approx 1.64$, almost all of it in the first steps. The simplest proof of the strong law applies the first lemma to the events $\{\lvert \bar X_n - \mu \rvert \ge \varepsilon\}$, whose probabilities fall fast enough when the tails are light, for example with a finite fourth moment; with only a finite mean the proof needs a more delicate argument.

Monte Carlo estimation#

Integrals become expectations by choosing a distribution. The area of the quarter disk $x^2 + y^2 \le 1$ in the unit square is $P(U_1^2 + U_2^2 \le 1) = \pi/4$ for independent uniform $U_1, U_2$, so $\pi$ is four times the share of random points that land inside:

Figure

π by Monte Carlo

With a million points the estimate is $3.1413$, with standard error $0.0016$: Monte Carlo is slow, three correct digits for a million samples, and each extra digit costs a hundred times more. Its strength is that the rate $1/\sqrt n$ does not depend on the dimension. A grid with 10 points per axis needs $10^{10}$ points in 10 dimensions; a Monte Carlo estimate of the volume of the 10-dimensional unit ball, $\pi^5/120 = 2.55$, from a million random points of the cube $[-1, 1]^{10}$ is $2.53$. The relative standard error, 2%, is large (the realized error of this run is 0.9%) because only one point in 400 lands inside the ball: estimating the probability of a rare event by plain sampling is expensive, which is the problem importance sampling solves.

The same principle runs through machine learning. A loss on held-out data is a Monte Carlo estimate of the risk, a minibatch gradient a Monte Carlo estimate of the full gradient, a validation accuracy a Monte Carlo estimate of the probability of a correct answer, and the reparametrized loss of a variational autoencoder a Monte Carlo estimate of an expectation over latent codes.

When averages do not converge#

The law needs a finite mean, and its speed depends on the tails.

• No mean. The Cauchy distribution of the expectation note has none, and the average of $n$ Cauchy samples is again standard Cauchy, exactly as spread out as a single sample. Its averages over $10^2$, $10^4$ and $10^6$ samples in the code below are $-0.01$, $-0.33$ and $2.86$: more data does not help.
• A mean but no variance. A Pareto distribution with tail exponent $1.5$ has mean 3 and infinite variance. Its averages do converge, but slowly and in jumps, whenever a single huge value arrives: $3.22$, $3.74$, $3.04$ for the same sample sizes. Incomes, file sizes and the number of followers of an account behave like this, and so do averages computed from them.
• Dependence. Averages of dependent values, such as the successive states of a Markov chain, still converge under conditions studied in the Markov chains note, but the effective number of samples is usually smaller than $n$.

What changed#

	Before	With the law of large numbers
Probability as frequency	an observation	a theorem: $\bar X_n \to \mu$ almost surely
Convergence	intuitive	in probability (weak law) and almost surely (strong law), separated by Borel–Cantelli
Expectations and integrals	computed analytically	estimated by averaging samples, error $\sigma/\sqrt n$ in any dimension
Failure modes		no mean (Cauchy), infinite variance (slow), dependence (fewer effective samples)

Limitations#

• The law says nothing about the rate. How close $\bar X_n$ is to $\mu$ for a given $n$ needs the inequalities of the previous note or the central limit theorem of the next.
• Identically distributed is an assumption. If the distribution drifts, as user behaviour does between training and deployment, the average converges to the average of a changing target, or to nothing.
• Monte Carlo is slow for rare events and high precision. Relative error grows as the target probability shrinks; variance reduction, importance sampling and quasi-random sequences exist for that reason.

What you should remember#

• Weak law: $P(\lvert \bar X_n - \mu \rvert \ge \varepsilon) \to 0$ (convergence in probability); strong law: $P(\bar X_n \to \mu) = 1$ (almost sure convergence). Both need only a finite mean.
• Borel–Cantelli: $\sum P(A_n) < \infty$ means that almost surely only finitely many $A_n$ happen; for independent events, $\sum P(A_n) = \infty$ means that almost surely infinitely many happen.
• Monte Carlo: $E[g(X)] \approx \frac1n \sum g(X_i)$ converges when $E\lvert g(X) \rvert < \infty$; with a finite variance its standard error is $\sigma/\sqrt n$, $\sigma$ the standard deviation of $g(X)$, at a rate independent of the dimension.
• Without a mean (Cauchy) averages do not settle; with infinite variance they settle slowly and in jumps.

Code#

The two kinds of convergence, two Monte Carlo integrals, and averages that misbehave:

import numpy as np, math

rng = np.random.default_rng(0)

# independent events with P(A_n) = 1/n and 1/n^2, n up to 100,000: how many happen, 1,000 runs each
N = 100_000
n = np.arange(1, N + 1)
c1 = [int((rng.random(N) < 1 / n).sum()) for _ in range(1000)]
c2 = [int((rng.random(N) < 1 / n**2).sum()) for _ in range(1000)]
print(round(float(np.mean(c1)), 2), round(float(np.sum(1 / n)), 2), round(float(np.mean(c2)), 2), round(math.pi**2 / 6, 2))   # 12.05 12.09 1.68 1.64
print(1 - 10_000 / 100_000, round(1 - float(np.prod(1 - 1 / n[100:] ** 2)), 4))   # 0.9 0.0099: an event after 10^4 (1/n); after 100 (1/n^2)

# Monte Carlo: pi from the unit square, and the volume of the 10-dimensional unit ball
u = rng.random((2, 1_000_000))
print(round(4 * float(((u**2).sum(axis=0) <= 1).mean()), 4), round(4 * math.sqrt(math.pi / 4 * (1 - math.pi / 4) / 1e6), 4))   # 3.1413 0.0016
x = rng.uniform(-1, 1, (1_000_000, 10))
share = float(((x**2).sum(axis=1) <= 1).mean())
print(round(share * 2**10, 3), round(math.pi**5 / 120, 3), round(math.sqrt(share * (1 - share) / 1e6) / share, 3))   # 2.527 2.55 0.02

# no mean: Cauchy averages wander; a finite mean with infinite variance (Pareto, mean 3) converges slowly
c = rng.standard_cauchy(10**6)
p = rng.pareto(1.5, 10**6) + 1
print([round(float(c[:k].mean()), 2) for k in (10**2, 10**4, 10**6)], [round(float(p[:k].mean()), 2) for k in (10**2, 10**4, 10**6)])
# [-0.01, -0.33, 2.86] [3.22, 3.74, 3.04]

Next#

The law of large numbers says where the average goes; the $1/\sqrt n$ rule says how far it typically is. Neither says what shape its fluctuations have. The next note shows that the shape is almost always the same, the normal distribution, whatever the distribution of the individual terms, and turns that fact into error bars.

References#

• J. K. Blitzstein, J. Hwang. Introduction to Probability, 2nd ed. CRC Press, 2019. Chapter 10 (the law of large numbers).
• W. Feller. An Introduction to Probability Theory and Its Applications, Vol. 1, 3rd ed. Wiley, 1968. Chapter X (the law of large numbers), Chapter VIII (the Borel–Cantelli lemmas).
• A. N. Shiryaev. Probability-1, 3rd ed. Springer, 2016. Chapter 2, §10 (kinds of convergence).
• A. N. Shiryaev. Probability-2, 3rd ed. Springer, 2019. Chapter 4, §1–3 (the zero–one laws and the strong law of large numbers).
• C. P. Robert, G. Casella. Monte Carlo Statistical Methods, 2nd ed. Springer, 2004. Chapter 3 (Monte Carlo integration).