The Central Limit Theorem

Delivery times are skewed: most orders arrive in 20 to 30 minutes, a few take an hour. Rainfall is skewed. Incomes are skewed. Yet the methods in the next lessons, confidence intervals and tests, use the normal curve almost everywhere. How can that work? The answer is the central limit theorem (CLT), probably the most important result in statistics.

In this lesson you will see the CLT happen in simulations, learn exactly what it says (and what it does not say), and use it to calculate the probability that an average falls in a range.

What the theorem says

Take random samples of size n from any population that has a mean μ and a finite standard deviation σ. Compute the mean x̄ of each sample. Then, as n grows:

1. the sample means are centred on the population mean:   mean of x̄ = μ
2. their spread shrinks with the square root of n:        SD of x̄ = σ / √n
3. their distribution gets closer and closer to normal:   x̄ ≈ N(μ, σ²/n)

Points 1 and 2 are true for any n. Point 3 is the surprising part: it holds whatever the shape of the population, as long as n is large enough. The collection of all possible sample means is called the sampling distribution of the mean.

Watching it happen

This simulation invents a strongly right-skewed population: delivery times of 100,000 orders, with a minimum of 10 minutes and a long tail. It then draws 10,000 samples for each sample size and measures their means.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
# Invented, strongly right-skewed population: delivery times (minutes) of 100,000 orders
population = rng.exponential(scale=20, size=100_000) + 10
mu, sigma = population.mean(), population.std()
print(f"Population: mean {mu:.2f}, SD {sigma:.2f}, skew {stats.skew(population):.2f}")

for n in [1, 2, 5, 30, 100]:
    means = rng.choice(population, size=(10_000, n)).mean(axis=1)
    print(f"n = {n:>3}: mean of means {means.mean():6.2f}   SD of means {means.std():5.2f}"
          f"   sigma/sqrt(n) {sigma / np.sqrt(n):5.2f}   skew {stats.skew(means):5.2f}")

Population: mean 29.99, SD 19.85, skew 1.96
n =   1: mean of means  29.99   SD of means 19.68   sigma/sqrt(n) 19.85   skew  1.89
n =   2: mean of means  29.95   SD of means 14.03   sigma/sqrt(n) 14.04   skew  1.38
n =   5: mean of means  30.00   SD of means  8.90   sigma/sqrt(n)  8.88   skew  0.86
n =  30: mean of means  30.03   SD of means  3.65   sigma/sqrt(n)  3.62   skew  0.38
n = 100: mean of means  29.99   SD of means  1.99   sigma/sqrt(n)  1.99   skew  0.23

Read the table column by column:

• Mean of means stays at about 30 for every n: sample means are centred on μ.
• SD of means matches σ / √n closely: 19.85 / √30 = 3.62, and the simulation gives 3.65.
• Skew falls from almost 2 (strongly skewed) to 0.23 (nearly symmetric) as n grows.

rng.choice(population, size=(10_000, n)) draws a 10,000 x n table of values, so each row is one sample; .mean(axis=1) averages each row. Sampling here is with replacement, which makes no practical difference when the population is this large.

Seeing the shape change

Numbers are convincing, but pictures are more so. Plot the sample means for three sample sizes.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(7)
population = rng.exponential(scale=20, size=100_000) + 10

fig, axes = plt.subplots(1, 3, figsize=(10, 3))
for ax, n in zip(axes, [1, 5, 30]):
    means = rng.choice(population, size=(10_000, n)).mean(axis=1)
    ax.hist(means, bins=40, color="#7c3aed")
    ax.set_title(f"Means of samples, n = {n}")
    ax.set_xlabel("Minutes")
plt.show()

With n = 1 the histogram is just the skewed population: a wall on the left at 10 minutes and a long tail to the right. With n = 5 the hump moves towards the middle. With n = 30 it is close to a bell curve, much narrower, centred on 30 minutes. Notice that the x-axis ranges shrink too: that is σ / √n at work.

Dice: from flat to bell-shaped

A single die is the opposite of bell-shaped: all six faces are equally likely. Average several dice, and the CLT takes over again.

import numpy as np

rng = np.random.default_rng(1)
for n in [1, 2, 10, 50]:
    avg = rng.integers(1, 7, size=(100_000, n)).mean(axis=1)
    within = (np.abs(avg - 3.5) < 0.5).mean()
    print(f"Average of {n:>2} dice: SD {avg.std():.3f} (theory {np.sqrt(35 / 12 / n):.3f}),"
          f" share within 0.5 of 3.5: {within:.3f}")

Average of  1 dice: SD 1.709 (theory 1.708), share within 0.5 of 3.5: 0.000
Average of  2 dice: SD 1.209 (theory 1.208), share within 0.5 of 3.5: 0.167
Average of 10 dice: SD 0.539 (theory 0.540), share within 0.5 of 3.5: 0.592
Average of 50 dice: SD 0.242 (theory 0.242), share within 0.5 of 3.5: 0.957

The variance of one die is 35/12, as you computed in Random Variables, so the SD of the average of n dice is √(35/12 / n). One die can never be within 0.5 of 3.5 (it shows whole numbers), but the average of 50 dice lands there 95.7% of the time.

Using the CLT to calculate probabilities

Because sample means are approximately normal, you can answer questions about averages with scipy.stats.norm, using σ / √n as the SD. This quantity is so important that it has its own name, the standard error, and its own lesson.

Suppose delivery times have mean 30 minutes and SD 20 minutes (like the population above). A restaurant will review its partner if the average of its next 36 deliveries is over 35 minutes. How likely is that, if nothing has changed?

import numpy as np
from scipy import stats

mu, sigma, n = 30, 20, 36
se = sigma / np.sqrt(n)
print(f"SE = {sigma} / sqrt({n}) = {se:.3f}")
print(f"CLT approximation: P(mean > 35) = {stats.norm.sf(35, mu, se):.4f}")

rng = np.random.default_rng(5)
sims = (rng.exponential(20, size=(200_000, n)) + 10).mean(axis=1)
print(f"Simulation:        P(mean > 35) = {(sims > 35).mean():.4f}")

SE = 20 / sqrt(36) = 3.333
CLT approximation: P(mean > 35) = 0.0668
Simulation:        P(mean > 35) = 0.0746

The average of 36 deliveries has SD 20 / 6 = 3.33 minutes, so 35 minutes is 1.5 standard errors above 30, and the normal curve gives 6.7%. The simulation from the true skewed population gives 7.5%. The CLT answer is close but not exact, because at n = 36 the sampling distribution of such a skewed population still has a little right skew (you saw skew 0.38 at n = 30). For less skewed data the match is better.

How large must n be?

A common rule of thumb is n ≥ 30, but it is only a rough guide:

• For roughly symmetric populations (marks, heights), n = 10 or even less can be enough.
• For strongly skewed populations (incomes, rainfall, insurance claims), you may need n in the hundreds.
• For proportions, a common check is that n x p and n x (1 - p) are both at least about 10.

When in doubt, simulate, exactly as you did above.

What the CLT does not say

• It does not say the data become normal. Individual delivery times stay skewed, however many you collect. Only the averages become normal.
• It does not fix bias. If your sample is biased, its mean is centred on the wrong value, with a nice bell shape around it.
• It needs a finite variance. Some distributions have tails so heavy that their variance is infinite, and then averages never settle down.

The standard example of the last point is the Cauchy distribution:

import numpy as np

rng = np.random.default_rng(3)
for n in [10, 1_000, 100_000]:
    means = [rng.standard_cauchy(n).mean() for _ in range(5)]
    print(f"n = {n:>7,}:", "  ".join(f"{m:7.2f}" for m in means))

n =      10:    0.47    -0.63     0.52    -8.77     0.89
n =   1,000:    8.29     0.07     0.41    -0.55     0.11
n = 100,000:   -1.92    -0.19    -1.33     0.09    -0.67

Even with 100,000 values, the means still jump around, because now and then one enormous value appears and drags the average with it. Real data rarely behaves this badly, but very heavy-tailed data (some financial returns, for example) can come close, so it is worth knowing the CLT has limits.

Common mistakes

• Saying the CLT makes the data normal. It is about the distribution of sample means (and sums), not the raw values.
• Treating n = 30 as a magic number. Skewed populations need more.
• Confusing σ and σ / √n. σ describes individual values; σ / √n describes sample means.
• Forgetting that the samples must be random and independent. The theorem says nothing about biased or clustered samples.

Key points

• Sample means are centred on μ and have standard deviation σ / √n.
• As n grows, the distribution of sample means approaches a normal distribution, whatever the population's shape.
• This is why normal-based confidence intervals and tests work for averages of skewed data.
• The larger the skew, the larger n must be; simulation is the honest check.
• The CLT needs random, independent samples and a finite variance.

import numpy as np
from scipy import stats

rng = np.random.default_rng(11)
population = rng.gamma(shape=2, scale=6, size=50_000)
sigma = population.std()
print(f"Population mean {population.mean():.2f}, SD {sigma:.2f}, skew {stats.skew(population):.2f}")

for n in [4, 16, 64]:
    means = rng.choice(population, size=(20_000, n)).mean(axis=1)
    print(f"n = {n:>2}: mean {means.mean():.2f}, SD {means.std():.3f}, "
          f"sigma/sqrt(n) {sigma / np.sqrt(n):.3f}, skew {stats.skew(means):.2f}")

Population mean 11.96, SD 8.45, skew 1.39
n =  4: mean 11.96, SD 4.242, sigma/sqrt(n) 4.227, skew 0.68
n = 16: mean 11.94, SD 2.110, sigma/sqrt(n) 2.114, skew 0.33
n = 64: mean 11.96, SD 1.054, sigma/sqrt(n) 1.057, skew 0.16