Lesson 2.2: Law of Large Numbers, Sampling Distributions and Confidence Intervals

Written lesson

Setup: the patient population

Python · Colab cell

!pip install -q numpy scipy polars matplotlib seaborn
import numpy as np, scipy.stats as stats, polars as pl
import matplotlib.pyplot as plt, seaborn as sns

rng = np.random.default_rng(101)
sns.set_theme(style="whitegrid")
N = 50_000
population = pl.DataFrame({
    "age": np.clip(rng.normal(45, 12, N), 18, 90),
    "cost": rng.lognormal(7.5, 1.0, N),          # annual health-care cost, very skewed
})
TRUE_MEAN_AGE = population["age"].mean()          # the population value we try to recover
print(f"True mean age: {TRUE_MEAN_AGE:.3f}")

The Law of Large Numbers

The LLN says that as a sample grows, its mean gets closer to the population mean. That is why a pilot with 30 participants gives volatile impact estimates while a national roll-out does not.

Python · Colab cell

sample_sizes = [10, 30, 100, 300, 1_000, 3_000, 10_000]
means = [population.sample(n=n, seed=i)["age"].mean() for i, n in enumerate(sample_sizes)]

for n, m in zip(sample_sizes, means):
    print(f"N = {n:>6,}: mean age = {m:6.2f}   (error {m - TRUE_MEAN_AGE:+.2f})")

plt.figure(figsize=(9, 4))
plt.plot(sample_sizes, means, marker="o", label="Sample mean")
plt.axhline(TRUE_MEAN_AGE, color="red", ls="--", label="Population mean")
plt.xscale("log"); plt.xlabel("Sample size (log scale)"); plt.ylabel("Mean age")
plt.title("Law of Large Numbers"); plt.legend(); plt.show()

Sampling distributions and the standard error

One sample gives one estimate. Draw many samples of the same size and you get a whole distribution of estimates: the sampling distribution. Its standard deviation is the standard error (SE), the typical distance between an estimate and the truth.

Python · Colab cell

ages = population["age"].to_numpy()

def sampling_distribution(n, reps=2_000):
    return np.array([rng.choice(ages, size=n, replace=False).mean() for _ in range(reps)])

for n in [25, 100, 400]:
    est = sampling_distribution(n)
    formula_se = ages.std() / np.sqrt(n)
    print(f"n={n:>3}: simulated SE = {est.std():.3f} | formula sigma/sqrt(n) = {formula_se:.3f}")

Quadrupling the sample halves the standard error: precision grows with the square root of n. That rule drives every power calculation you will do in Module 4.

The Central Limit Theorem

Health-care costs are extremely skewed. Yet the mean of a sample of costs is approximately normal once n is moderately large. That is the Central Limit Theorem (CLT), and it is why normal-based confidence intervals work for so many estimators.

Python · Colab cell

costs = population["cost"].to_numpy()
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
sns.histplot(costs, bins=60, ax=axes[0]).set(title="Individual costs (skewed)")
for ax, n in zip(axes[1:], [5, 200]):
    m = [rng.choice(costs, n).mean() for _ in range(3_000)]
    sns.histplot(m, bins=50, ax=ax).set(title=f"Means of samples of n={n}")
plt.tight_layout(); plt.show()

Confidence intervals and coverage

A 95% confidence interval is estimate ± 1.96 × SE. The "95%" is a property of the procedure: if you repeated the study many times, about 95% of the intervals would contain the truth. We can check that directly.

Python · Colab cell

def coverage(n, reps=2_000, use_t=False):
    hits = 0
    for _ in range(reps):
        s = rng.choice(ages, n, replace=False)
        se = s.std(ddof=1) / np.sqrt(n)
        crit = stats.t.ppf(0.975, df=n - 1) if use_t else 1.96
        hits += (s.mean() - crit * se <= TRUE_MEAN_AGE <= s.mean() + crit * se)
    return hits / reps

for n in [5, 30, 200]:
    print(f"n={n:>3}: coverage (z) = {coverage(n):.3f} | coverage (t) = {coverage(n, use_t=True):.3f}")

With tiny samples the normal approximation under-covers; the t-distribution corrects for it. Coverage simulations like this one are how methodologists test every estimator in this course, and you will reuse the idea for permutation tests and clustered standard errors.

Exercise

Repeat the coverage check for the mean cost (the skewed variable) with n = 10, 50 and 500. Where does the CLT approximation break down, and what does that imply for evaluating a programme on health spending with a small sample?

Show a solution

Python · Colab cell

TRUE_MEAN_COST = costs.mean()
for n in [10, 50, 500]:
    hits = 0
    for _ in range(2_000):
        s = rng.choice(costs, n, replace=False)
        se = s.std(ddof=1) / np.sqrt(n)
        hits += abs(s.mean() - TRUE_MEAN_COST) <= stats.t.ppf(0.975, n - 1) * se
    print(f"n={n:>3}: coverage = {hits / 2000:.3f}")
# Coverage is well below 95% for small n: skewed outcomes need larger samples, a log transform,
# or a bootstrap before you trust a confidence interval.

Checklist and recap

I can demonstrate the LLN and the CLT with a simulation.
I know the standard error is the spread of the sampling distribution, and that it shrinks with √n.
I can explain what "95% confidence" promises and check coverage myself.