05 - Central Limit Theorem and Law of Large Numbers

A research institute wants to estimate the share of voters in favour of a proposal, with a margin of error of three percentage points and 95% confidence. How many people must it interview? The answer is 1,068, and reaching it takes two results that are the heart of the theory: the law of large numbers, which guarantees that the observed frequency approaches the true probability, and the central limit theorem, which says how closely it approaches and with what shape of error.

Here you state them, test them on examples where the error can be computed exactly (coins, dice, sums of exponential variables), see where they fail (the Cauchy distribution) and close with the poll, checking with an exact calculation whether those 1,068 respondents are really enough. You need the continuous variables, in particular the normal, and the concept of independence from article 01.

What "the mean converges" means

Let X₁, X₂, … be independent variables with the same distribution (called i.i.d.), with mean μ and variance σ², and let X̄_n be the mean of the first n. This mean is itself a random variable, with mean μ and variance σ²/n (the variance of a sum of independent variables is the sum of the variances, then you divide by n²). "Converges to μ" can mean different things.

Convergence in probability requires that, for every margin ε > 0, the probability that X̄_n is farther than ε from μ tends to zero. Almost sure convergence is stronger: it concerns the single sequence of observations and says that with probability 1 the running means end up staying close to μ. Convergence in distribution, finally, says nothing about the values but only about their law: the distribution function of X_n tends to that of the limit at the points where the latter is continuous.

The weak law of large numbers asserts the convergence in probability of X̄_n to μ, and with finite variance it is proved in one line with Chebyshev's inequality, P(|Y − E[Y]| ≥ ε) ≤ Var(Y)/ε² (theorem 8.1 of chapter 8 of Grinstead and Snell): applied to X̄_n it gives P(|X̄_n − μ| ≥ ε) ≤ σ²/(nε²), which tends to zero.

The strong law, due to Kolmogorov, asserts almost sure convergence and holds under the sole assumption that the mean exists; an elementary proof by Etemadi (1981) exists, which I do not reproduce here. The difference between the two laws is concrete: the weak one guarantees that at every large n you are probably close to μ, the strong one that your single data series will not move away again, from some point onwards.

The first program, which uses the count of successes of the binomial distribution, compares three numbers for the fair coin (μ = 1/2, σ² = 1/4): the exact probability that the frequency of heads is at least 0.05 away from the true value, computed by summing the binomial distribution, Chebyshev's upper bound and the normal approximation supplied by the central limit theorem. Then it follows a single trajectory, the running mean of a million rolls of a die (μ = 3.5, σ² = 35/12), compared with the typical error σ/√n.

import numpy as np
from scipy.stats import binom, norm

# Mean of n fair coin tosses: how likely is it that the frequency of heads departs from 1/2 by at least eps?
p, eps = 0.5, 0.05
print(f"{'n':>5} {'exact':>9} {'Chebyshev':>10} {'CLT':>9}")
for n in (100, 400, 1600):
    k = np.arange(n + 1)
    esatta = binom.pmf(k, n, p)[np.abs(k / n - p) >= eps - 1e-12].sum()   # exact sum over the binomial
    cheb = p * (1 - p) / (n * eps**2)                                    # sigma^2 / (n eps^2)
    clt = 2 * (1 - norm.cdf(eps * n**0.5 / (p * (1 - p)) ** 0.5))        # normal approximation
    print(f"{n:5d} {esatta:9.5f} {min(cheb, 1):10.5f} {clt:9.5f}")

# One trajectory: running mean of die rolls, with typical error sigma/sqrt(n), sigma^2 = 35/12
rng = np.random.default_rng(101)
lanci = rng.integers(1, 7, size=10**6)
media = np.cumsum(lanci) / np.arange(1, 10**6 + 1)
sigma = (35 / 12) ** 0.5
print(f"{'n':>8} {'running mean':>15} {'|mean-3.5|':>12} {'sigma/sqrt(n)':>14}")
for n in (10, 100, 10**3, 10**4, 10**5, 10**6):
    print(f"{n:8d} {media[n - 1]:15.4f} {abs(media[n - 1] - 3.5):12.4f} {sigma / n**0.5:14.4f}")

# Output:
#     n     exact  Chebyshev       CLT
#   100   0.36820    1.00000   0.31731
#   400   0.05104    0.25000   0.04550
#  1600   0.00007    0.06250   0.00006
#        n    running mean   |mean-3.5|  sigma/sqrt(n)
#       10          3.5000       0.0000         0.5401
#      100          3.4000       0.1000         0.1708
#     1000          3.4360       0.0640         0.0540
#    10000          3.4993       0.0007         0.0171
#   100000          3.5047       0.0047         0.0054
#  1000000          3.5023       0.0023         0.0017

With 100 tosses the exact probability of a deviation of at least five percentage points is 0.368, and Chebyshev's bound is 1: true but it tells you nothing. With 400 tosses the exact one is 0.0510 against a Chebyshev of 0.25, with 1600 it is 0.00007 against 0.0625. The bound is a guarantee valid for any distribution with finite variance, and precisely for this reason it is loose (here, by almost a thousand times); the normal approximation (0.317, 0.0455, 0.00006) follows the exact one very closely, and its error is the subject of the next section.

The die's trajectory shows the strong law at work: the deviation from the mean 3.5 oscillates (0.100 at n = 100, 0.064 at 1,000, 0.0007 at 10,000) but stays of the order of σ/√n, which halves each time n quadruples, that is it drops by a factor of 10 only by multiplying the observations by a hundred.

When the law of large numbers does not hold: the Cauchy

The two laws require the mean to exist. The standard Cauchy distribution, with density 1/(π(1 + x²)), has tails so heavy that the integral of |x| diverges: it has no expected value. In this case the sample mean of n observations is still a standard Cauchy, whatever n is: averaging is of no use.

The quartiles of the standard Cauchy are −1 and +1, so the program estimates the quartiles and the median of the means of n observations, over 2,000 repetitions, for n equal to 1, 100 and 10,000, and compares them with those of a distribution that does have a mean, the exponential with mean 1.

import numpy as np

# Standard Cauchy: no mean. The mean of n observations is still standard Cauchy, for every n.
# The quartiles of the standard Cauchy are -1 and +1 (exact).
rng = np.random.default_rng(103)
R = 2000                                          # replicates for each n
print(f"{'n':>7} {'q25 of the means':>17} {'median':>9} {'q75':>8}   (expected -1, 0, +1; typical error about 0.06)")
for n in (1, 100, 10_000):
    medie = rng.standard_cauchy((R, n)).mean(axis=1)
    q25, q50, q75 = np.quantile(medie, [0.25, 0.5, 0.75])
    print(f"{n:7d} {q25:16.3f} {q50:9.3f} {q75:8.3f}")

# Comparison with a distribution that has a mean: Exp(1), mean 1. The quartiles of the means tighten around 1.
for n in (1, 100, 10_000):
    medie = rng.exponential(1.0, (R, n)).mean(axis=1)
    q25, q50, q75 = np.quantile(medie, [0.25, 0.5, 0.75])
    print(f"Exp  n={n:6d}: q25 {q25:.4f}  median {q50:.4f}  q75 {q75:.4f}")

# Output:
#       n  q25 of the means    median      q75   (expected -1, 0, +1; typical error about 0.06)
#       1           -0.996    -0.046    0.907
#     100           -1.058    -0.034    1.040
#   10000           -1.070    -0.036    1.011
# Exp  n=     1: q25 0.2980  median 0.7037  q75 1.3764
# Exp  n=   100: q25 0.9278  median 0.9951  q75 1.0655
# Exp  n= 10000: q25 0.9942  median 1.0006  q75 1.0071

For the Cauchy the quartiles of the means stay at about −1 and +1 and the median at about zero with n = 1, 100 and 10,000 (for example −1.070 and 1.011 at n = 10,000), without any narrowing: the estimation error at 10,000 observations is as large as that of a single one.

For the exponential, on the other hand, the interval between the quartiles shrinks from (0.298; 1.376) to (0.928; 1.066) and finally to (0.994; 1.007), that is the width falls like 1/√n. The message is operational: before averaging quantities with extreme tails (insurance payouts, response times, returns in crashes) check that the mean exists; if it does not, or has infinite variance, the theorem that follows does not apply.

The central limit theorem and its error

The law of large numbers says that X̄_n approaches μ; the central limit theorem (of Lindeberg-Lévy) says how: if the X_i are i.i.d. with finite and positive variance, the standardized sum (S_n − nμ)/(σ√n), where S_n = X₁ + … + X_n, has a distribution function that tends to Φ, that of the standard normal, at every point.

The starting shape of the distribution does not matter, as long as the variance is finite. The theorem is a limit and does not say how large the error is at fixed n. The Berry-Esseen inequality (Berry 1941) does: the maximum distance between the distribution function of the standardized sum and Φ is at most Cρ/(σ³√n), where ρ = E[|X − μ|³] and C is a universal constant, for which Shevtsova (2011) proved the value 0.4748 in the i.i.d. case.

To measure the error without simulating you need a case in which the law of the sum is known: the sum of n exponentials of mean 1 has exactly a Gamma distribution of shape n, so the Kolmogorov distance is computed with the exact distribution function, on a fine grid. The program computes it for n from 1 to 1,000, computes ρ by numerical integration (the closed form is 12/e − 2), compares it with the Berry-Esseen bound, and finally simulates 200,000 sums with n = 30 to see the same distance from a sample.

import numpy as np
from scipy.stats import gamma, norm, kstest
from scipy.integrate import quad

# X ~ Exp(1): mean 1, variance 1, skewness 2. S_n = X_1+...+X_n has EXACT law Gamma(n, 1).
rho, _ = quad(lambda x: abs(x - 1) ** 3 * np.exp(-x), 0, np.inf, limit=200)   # E|X - mu|^3
print(f"rho = E|X-1|^3 = {rho:.6f}  (formula 12/e - 2 = {12 / np.e - 2:.6f})")
C = 0.4748                                       # Berry-Esseen constant for i.i.d. sums (Shevtsova 2011)
z = np.linspace(-6, 6, 120001)
print(f"{'n':>5} {'exact distance':>16} {'B-E bound':>16} {'dist*sqrt(n)':>13}")
for n in (1, 5, 30, 200, 1000):
    # F_n(z) = P((S_n - n)/sqrt(n) <= z) = Gamma(n,1).cdf(n + z sqrt(n))
    Fn = gamma.cdf(n + z * n**0.5, a=n)
    dist = np.max(np.abs(Fn - norm.cdf(z)))
    print(f"{n:5d} {dist:16.5f} {C * rho / n**0.5:16.5f} {dist * n**0.5:13.4f}")

# Confirmation by simulation at n = 30: 200000 sums, Kolmogorov-Smirnov distance from the standard normal
rng = np.random.default_rng(107)
n, M = 30, 200_000
exact30 = np.max(np.abs(gamma.cdf(n + z * n**0.5, a=n) - norm.cdf(z)))
S = rng.exponential(1.0, size=(M, n)).sum(axis=1)
Zs = (S - n) / n**0.5
D = kstest(Zs, norm.cdf).statistic
print(f"n=30, {M} simulated sums: D = {D:.5f}, exact distance {exact30:.5f}, difference {D - exact30:+.5f} (typical noise 1/sqrt(M) = {1 / M**0.5:.5f})")
terzo = Zs**3
print(f"third moment of the standardized sums: {terzo.mean():.3f} +/- {terzo.std() / M**0.5:.3f} (exact 2/sqrt(30) = {2 / 30**0.5:.3f})")

# Output:
# rho = E|X-1|^3 = 2.414553  (formula 12/e - 2 = 2.414553)
#     n   exact distance        B-E bound  dist*sqrt(n)
#     1          0.15866          1.14643        0.1587
#     5          0.05963          0.51270        0.1333
#    30          0.02429          0.20931        0.1330
#   200          0.00940          0.08106        0.1330
#  1000          0.00421          0.03625        0.1330
# n=30, 200000 simulated sums: D = 0.02781, exact distance 0.02429, difference +0.00352 (typical noise 1/sqrt(M) = 0.00224)
# third moment of the standardized sums: 0.341 +/- 0.010 (exact 2/sqrt(30) = 0.365)

The exact distance falls like 1/√n, from 0.159 at n = 1 to 0.0243 at n = 30 and to 0.0042 at n = 1,000, and the Berry-Esseen bound always covers it but with a margin of almost nine times (for n = 1 it is 1.146, so it tells you nothing). The product of the distance by √n stabilizes at 0.1330: it is the first term of the Edgeworth correction, γ/(6√(2π)) with γ = 2 the skewness of the exponential, a result that I do not derive here, and that explains why the slowness of the approximation comes from the asymmetry.

The simulation at n = 30 gives D = 0.0278 against an exact distance of 0.0243: the difference (0.0035) is about one and a half times the typical noise, because the sample distance is the maximum of noisy deviations and on average cannot lie below the true one.

The third moment of the simulated sums is 0.341 ± 0.010 against the exact 2/√30 = 0.365: it is two and a half times the standard error, a fluctuation of this seed (five other seeds, tried separately, give values between 0.353 and 0.391 around the same exact value). The rule "with n ≥ 30 the normal holds" is therefore a habit, not a theorem: for an asymmetric distribution like this one the 2.4% error in the distribution function at n = 30 may be too much.

Worked exercise: how many respondents are needed

Back to the poll. Each respondent answers yes with unknown probability p; the frequency p̂ of yeses among n people has mean p and variance p(1 − p)/n, and by the central limit theorem it is approximately normal. The interval p̂ ± m contains p with probability about 95% when m = 1.96√(p(1 − p)/n), where 1.96 is the 97.5% quantile of the standard normal, that is n = 1.96²p(1 − p)/m².

Since p is not known, you use the worst case p = 1/2, which maximizes p(1 − p), with m = 0.03: n = 1,067.07, to be rounded up to 1,068.

The program checks this answer with the exact calculation of the coverage probability, summing the binomial distribution, for four true values of p and for two kinds of interval: the fixed-margin one of ±0.03 and the Wald one, whose width is estimated from the observed frequency, p̂ ± 1.96√(p̂(1 − p̂)/n). Finally it simulates 200,000 polls with p = 1/2 to confirm the coverage.

import numpy as np
from math import ceil
from scipy.stats import binom

# Margin of error +/-3 points, confidence 95%: n >= z^2 p(1-p) / m^2, worst case p = 1/2.
z, m = 1.959964, 0.03                             # z = 97.5% quantile of the standard normal
n = ceil(z**2 * 0.25 / m**2)
print(f"minimum n (worst case p=1/2): {z**2 * 0.25 / m**2:.2f} -> {n}")

k = np.arange(n + 1)
phat = k / n
print(f"{'true p':>7} {'coverage +/-3%':>16} {'Wald coverage':>15}")
for p in (0.5, 0.3, 0.1, 0.02):
    pesi = binom.pmf(k, n, p)
    fisso = pesi[np.abs(phat - p) <= m].sum()                              # interval phat +/- 0.03
    semi = z * np.sqrt(phat * (1 - phat) / n)                              # Wald: width estimated from phat
    wald = pesi[(phat - semi <= p) & (p <= phat + semi)].sum()
    print(f"{p:7.2f} {fisso:16.4f} {wald:15.4f}")

rng = np.random.default_rng(109)
R = 200_000
x = rng.binomial(n, 0.5, size=R) / n
cop = np.mean(np.abs(x - 0.5) <= m)
esatta = binom.pmf(k, n, 0.5)[np.abs(phat - 0.5) <= m].sum()
tol = 3 * (esatta * (1 - esatta) / R) ** 0.5
print(f"p=0.5, {R} simulated polls: coverage {cop:.4f} (exact {esatta:.4f} +/- {tol:.4f})")
assert abs(cop - esatta) < tol

# Output:
# minimum n (worst case p=1/2): 1067.07 -> 1068
#  true p   coverage +/-3%   Wald coverage
#    0.50           0.9533          0.9462
#    0.30           0.9674          0.9511
#    0.10           0.9988          0.9445
#    0.02           1.0000          0.9290
# p=0.5, 200000 simulated polls: coverage 0.9532 (exact 0.9533 +/- 0.0014)

With 1,068 respondents and p = 1/2 the exact coverage of the fixed margin is 0.9533, in line with the 95% required (the excess is due to the rounding of n), and the simulation gives 0.9532, within the tolerance of ±0.0014. The fixed margin is prudent: for p = 0.3 the coverage rises to 0.9674 and for p = 0.1 to 0.9988, because at that p the variance is smaller.

The Wald interval, on the other hand, shows where the normal approximation weakens: 0.9462 at p = 0.5, 0.9511 at 0.3, 0.9445 at 0.1 and only 0.9290 at p = 0.02, because when success is rare the distribution of p̂ is asymmetric and the estimate of the variance from the observed frequency is unreliable. This is a limit of the model, not of the calculation: all of this holds for a simple random sample with honest answers; a real poll has other sources of error, such as non-response and sample selection, which this calculation does not measure.

What is left out

This article treats the case of i.i.d. variables with finite variance. When the variables are not identically distributed you need Lindeberg's condition (1922), which requires that no term dominates the sum; the Lindeberg-Feller theorem makes it a necessary and sufficient condition.

I do not cover either the Edgeworth corrections beyond the first term, the versions for vectors (which lead to a multivariate normal, as in article 04), or the laws of large numbers for dependent variables. The proof of the theorem, by characteristic functions, requires tools of analysis that are not needed here. The next step is statistical inference, where the central limit theorem justifies confidence intervals.

Frequently asked questions

Does the law of large numbers say that after many heads tails must come up?

No. The tosses are independent and the coin has no memory: the frequency approaches 1/2 because future tosses dilute the accumulated deviation, not because they compensate for it. The absolute deviation can keep growing like √n even though the frequency converges.

From how many observations does the normal hold?

It depends on the shape of the distribution. For the exponential, the distance from the normal at 30 observations is 0.024; for symmetric distributions it is much smaller. The Berry-Esseen bound gives a guarantee, but a prudent one.

Why does the mean of Cauchy data not converge?

Because the distribution has no expected value: the extreme values are frequent enough to dominate the sum. The mean of n observations has exactly the same distribution as a single one.

What you can do now

You can tell convergence in probability, almost sure convergence and convergence in distribution apart, prove the weak law with Chebyshev and recognize a case, the Cauchy, in which neither law holds. You can state the central limit theorem, measure its error with an exact calculation and compare it with the Berry-Esseen bound, and size a sample, then check the coverage instead of trusting the formula alone. The next step is to exploit these results to estimate a parameter from data: estimators, their properties and confidence intervals.

The same series also includes Hypothesis Testing: t-test, z-test, p-value and Exploratory Data Analysis (EDA) with Python. The full list is on the series page, Probability and Statistics.