The law of large numbers · Probability & Statistics

Everything statistics does rests on one promise: that a sample average tells you something about the true mean. The Law of Large Numbers is that promise, made precise and proved.

The statement

Let X_1, X_2, \dots be independent, identically distributed with E[X_i] = \mu finite. Define the sample mean

$\bar X_n = \frac{1}{n}\sum_{i=1}^{n} X_i$

Weak Law (WLLN). \bar X_n \xrightarrow{p} \mu:

$P\big(|\bar X_n - \mu| > \varepsilon\big) \to 0 \quad \text{for every } \varepsilon > 0$

Strong Law (SLLN). \bar X_n \xrightarrow{a.s.} \mu:

$P\Big(\lim_{n\to\infty}\bar X_n = \mu\Big) = 1$

The two differ exactly in the mode of convergence from the last lesson. The weak law says: at any fixed large n, a big deviation is unlikely. The strong law says: your particular sequence of averages settles down and stays there.

The strong law is the one that matches intuition, and it holds under the same assumption — a finite mean.

Proving the weak law

The weak law falls out of two ingredients.

Chebyshev's inequality. For any random variable with finite variance:

$P\big(|Y - E[Y]| \ge \varepsilon\big) \le \frac{\operatorname{Var}(Y)}{\varepsilon^2}$

The variance of the sample mean. Using linearity and independence:

$E[\bar X_n] = \mu, \qquad \operatorname{Var}(\bar X_n) = \frac{1}{n^2}\sum_{i=1}^n \operatorname{Var}(X_i) = \frac{\sigma^2}{n}$

The \frac{1}{n^2} comes from \operatorname{Var}(aX) = a^2\operatorname{Var}(X); the sum of n identical variances is n\sigma^2; and independence is what allowed the variances to add at all.

Combine them:

$P\big(|\bar X_n - \mu| \ge \varepsilon\big) \le \frac{\sigma^2}{n\varepsilon^2} \longrightarrow 0$

Done. (This proof assumes a finite variance; the weak law actually holds with only a finite mean, but that needs characteristic functions.)

The key quantity: \sigma/\sqrt n

The standard deviation of the sample mean —

$\operatorname{SD}(\bar X_n) = \frac{\sigma}{\sqrt n}$

— is called the standard error, and it is the single most important formula in applied statistics.

Note the \sqrt n. To halve your uncertainty you need four times the data. To reduce it tenfold, a hundredfold increase. This is why precision is expensive, and why polls stop at ~1,000 respondents: the gain from 1,000 to 4,000 is only a halving, at four times the cost.

What it does not say

It says nothing about correcting past deviations. If you flip 10 heads in a row, the LLN does not predict extra tails to compensate. The coin is memoryless (§1). What actually happens is dilution: those 10 excess heads become negligible against n = 10{,}000, not cancelled.

Precisely: the proportion converges, while the absolute difference between heads and tails typically grows like \sqrt n. Believing otherwise is the gambler's fallacy.

It requires a finite mean. For the Cauchy distribution, E[X] doesn't exist and the LLN fails completely — we showed in §1 that \bar X_n has exactly the same distribution as a single observation, for every n. Averaging a million Cauchys is no better than taking one.

It says nothing about how fast. Convergence is guaranteed but the rate is not specified — that's the Central Limit Theorem's job, next lesson.

Worked example

You flip a fair coin. How many flips until you're 95% sure the observed proportion is within 1% of 0.5?

Each flip is Bernoulli(0.5), so \mu = 0.5 and \sigma^2 = 0.25.

Using Chebyshev with \varepsilon = 0.01 and requiring the bound \le 0.05:

$\frac{\sigma^2}{n\varepsilon^2} = \frac{0.25}{n(0.0001)} \le 0.05$

$n \ge \frac{0.25}{0.05 \times 0.0001} = 50{,}000$

So Chebyshev guarantees 50,000 flips suffice.

But that's very conservative. Chebyshev makes no distributional assumptions — it holds for any distribution with that variance, so it must cover the worst case. Using the Normal approximation from the next lesson:

$n \ge \left(\frac{1.96\,\sigma}{\varepsilon}\right)^2 = \left(\frac{1.96 \times 0.5}{0.01}\right)^2 = 9604$

About 9,600 — five times fewer. Chebyshev is the guarantee you can always make; the CLT is the answer you actually use.

Doing it in Python

Watch the convergence, and the \sqrt n rate:

import numpy as np

rng = np.random.default_rng(0)
flips = rng.random(1_000_000) < 0.5
running = np.cumsum(flips) / np.arange(1, len(flips) + 1)

print(f"{'n':>10} {'proportion':>12} {'|error|':>10} {'sigma/sqrt(n)':>15}")
for n in (10, 100, 1_000, 10_000, 100_000, 1_000_000):
    p = running[n - 1]
    print(f"{n:>10,} {p:>12.5f} {abs(p - 0.5):>10.5f} {0.5/np.sqrt(n):>15.5f}")

The error column tracks the standard error column closely — that's the LLN converging at exactly the rate \sigma/\sqrt n predicts.

Now the gambler's fallacy, dismantled numerically:

import numpy as np

rng = np.random.default_rng(1)
flips = rng.random(1_000_000) < 0.5
heads = np.cumsum(flips)
n = np.arange(1, len(flips) + 1)

print(f"{'n':>10} {'proportion':>12} {'heads - tails':>15}")
for k in (100, 10_000, 1_000_000):
    print(f"{k:>10,} {heads[k-1]/k:>12.5f} {2*heads[k-1] - k:>15,}")

print("\nThe PROPORTION converges to 0.5.")
print("The absolute SURPLUS does not shrink -- it typically grows like sqrt(n).")

And the failure case — Cauchy, where there's no mean to converge to:

import numpy as np

rng = np.random.default_rng(2)
normal = rng.standard_normal(1_000_000)
cauchy = rng.standard_cauchy(1_000_000)

n_run = np.arange(1, 1_000_001)
normal_run = np.cumsum(normal) / n_run
cauchy_run = np.cumsum(cauchy) / n_run

print(f"{'n':>10} {'normal mean':>14} {'cauchy mean':>14}")
for k in (10, 1_000, 100_000, 1_000_000):
    print(f"{k:>10,} {normal_run[k-1]:>14.5f} {cauchy_run[k-1]:>14.5f}")

print("\nNormal settles toward 0. Cauchy never settles -- no finite mean, no LLN.")

Your turn

1. A die is rolled 6,000 times. Roughly how many 6s do you expect, and what is the standard error of the proportion?

2. Chebyshev bound: how many samples to be 99% sure \bar X is within 0.1 of \mu, if \sigma = 2?

3. You've flipped 10 heads in a row. What's the probability the next flip is heads?

Solutions

1. Each roll is a 6 with probability 1/6, so

$E[\text{number of 6s}] = 6000 \times \tfrac16 = 1000$

For the proportion, the per-roll variance is p(1-p) = \frac16 \cdot \frac56 = \frac{5}{36}, so

$\text{SE} = \sqrt{\frac{p(1-p)}{n}} = \sqrt{\frac{5/36}{6000}} \approx 0.00481$

So the proportion is typically within about 0.5 percentage points of 0.1667 — equivalently, the count is typically within about $0.00481 \times 6000 \approx 29$ of 1,000.

2. Chebyshev with \varepsilon = 0.1, requiring the bound \le 0.01:

$\frac{\sigma^2}{n\varepsilon^2} = \frac{4}{n(0.01)} \le 0.01$

$n \ge \frac{4}{0.01 \times 0.01} = 40{,}000$

Again very conservative. The Normal-based answer would be \left(\frac{2.576 \times 2}{0.1}\right)^2 \approx 2654 — fifteen times smaller. Chebyshev's strength is that it needs no distributional assumption at all; that's also its weakness.

3. Exactly 0.5.

Flips are independent — the coin has no memory of the previous ten. The LLN does not say that tails become more likely to "balance things out".

What it says is that those 10 excess heads become irrelevant: after a million flips the proportion is \frac{500{,}005}{1{,}000{,}010} \approx 0.5000. The surplus was never cancelled, just diluted by the volume of later data.

If the question had been "is this coin fair?", 10 heads in a row is weak evidence against it — P = (0.5)^{10} \approx 0.001 under fairness. That's a hypothesis test (§5), and a completely different question from what the next flip does.

Check yourself in code

Simulate coin flips and confirm the running proportion converges to 0.5 while the absolute surplus of heads does not shrink.

Print exactly this:

n=100 proportion 0.47
n=10000 proportion 0.4953
n=1000000 proportion 0.50037
proportion converged: True
surplus shrank: False

Use default_rng(1) and 1000000 flips. Print each proportion rounded to 5 decimal places. Report proportion converged: True if the final error from 0.5 is under 0.001, and surplus shrank: False if the absolute head-minus-tail surplus at n = 10^6 is not smaller than at n = 100.

import numpy as np

rng = np.random.default_rng(1)
flips = rng.random(1_000_000) < 0.5
heads = np.cumsum(flips)

for k in (100, 10_000, 1_000_000):
    print(f"n={k} proportion", round(heads[k - 1] / k, 5))

# Report whether the final proportion is within 0.001 of 0.5, and whether
# the absolute surplus |2*heads - n| at n=10^6 is smaller than at n=100.

import numpy as np

rng = np.random.default_rng(1)
flips = rng.random(1_000_000) < 0.5
heads = np.cumsum(flips)

for k in (100, 10_000, 1_000_000):
    print(f"n={k} proportion", round(heads[k - 1] / k, 5))

final_error = abs(heads[-1] / 1_000_000 - 0.5)
print("proportion converged:", bool(final_error < 0.001))

surplus_small = abs(2 * heads[99] - 100)
surplus_big = abs(2 * heads[-1] - 1_000_000)
print("surplus shrank:", bool(surplus_big < surplus_small))

The sample mean converges to the true mean — in probability (weak law) and almost surely (strong law), provided the mean exists at all. The rate is governed by the standard error \sigma/\sqrt n, which is why quadrupling your data only halves your uncertainty. And nothing here corrects past deviations; they are diluted, never cancelled.

Next: the theorem that says not just that \bar X_n converges, but exactly what its distribution looks like on the way.

23. The law of large numbers