import matplotlib
matplotlib.use("Agg") # headless: render to buffers, never to a screen
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(42) # one generator, pinned, drives the whole notebook
plt.rcParams["figure.figsize"] = (9, 5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3nb1.4 — The LLN and the CLT, Shown Not Asserted
Week 1, Chapter 1.4. Companion notebook.
In Ch 1.3 we watched a sampling distribution tighten as N grew, but we never said why it concentrates or what shape it settles into. This notebook pays both debts by simulation:
- The Weak Law of Large Numbers (WLLN) — Sam rolls a fair die (\mu = 3.5) and we watch the running mean stop wandering and home in on the truth.
- The Central Limit Theorem (CLT) — Priya averages a violently skewed population (the exponential) and we watch the standardized mean morph into the standard normal bell as N = 1, 2, 5, 30, 100.
- When the CLT is slow — or fails — Devon’s fat-tailed crypto returns, modeled as Student’s t with \nu = 5 (slow), \nu = 2 (infinite variance), and \nu = 1 (Cauchy, where the LLN itself dies).
Everything is reproducible: one seeded generator, no animation (we use static matplotlib panels so the notebook runs headless), and scipy.stats for the reference curves.
1. The Weak Law of Large Numbers
The result in one sentence. As the sample size N grows, the sample mean \bar{x} stops wandering and converges in probability to the true population mean \mu: \bar{x} \xrightarrow{p} \mu.
The cleanest way to see this is not ten thousand histograms — it is to follow a single sample as it grows and plot the running mean after each new draw. Sam rolls a fair die, where the population mean is
\mu = \frac{1+2+3+4+5+6}{6} = 3.5.
We draw one die, then a second, then a third, and after each new roll recompute the average of everything so far.
N = 5000
rolls = rng.integers(1, 7, size=N) # uniform on {1,...,6}; integers high-end is exclusive
mu_die = 3.5
running_mean = np.cumsum(rolls) / np.arange(1, N + 1)
fig, ax = plt.subplots()
ax.plot(running_mean, lw=1.2, label="running sample mean")
ax.axhline(mu_die, color="crimson", ls="--", label=r"true $\mu = 3.5$")
ax.set_xscale("log") # log-x so the early wildness is visible
ax.set_xlabel("number of rolls (N)")
ax.set_ylabel("running sample mean")
ax.set_title("WLLN: a single running mean settling onto $\\mu$")
ax.legend()
fig.tight_layout()
fig.savefig("nb14_wlln_single.png", dpi=110)
plt.close(fig)
print(f"After {N} rolls the running mean is {running_mean[-1]:.4f} (true mu = {mu_die}).")The curve starts jagged — after three rolls the average could be anything from 1 to 6 — and then, as N climbs into the hundreds and thousands, flattens and presses against the dashed line. It never lands exactly on 3.5 and stays there; it keeps making ever-smaller corrections forever. The room it has to wander in keeps closing.
A sanity check worth running: the destination is a property of the population, but the route is random. Let us overlay several independent paths and confirm they all funnel toward the same line, with the funnel narrowing like 1/\sqrt{N}.
n_paths = 8
Npath = 5000
x_axis = np.arange(1, Npath + 1)
fig, ax = plt.subplots()
for _ in range(n_paths):
r = rng.integers(1, 7, size=Npath)
rm = np.cumsum(r) / x_axis
ax.plot(rm, lw=0.9, alpha=0.7)
# theoretical 1/sqrt(N) funnel: sd of the running mean is sigma/sqrt(N)
sigma_die = np.sqrt(np.mean((np.arange(1, 7) - mu_die) ** 2)) # population sd of a die
funnel = sigma_die / np.sqrt(x_axis)
ax.plot(x_axis, mu_die + 2 * funnel, "k--", lw=1.2, label=r"$\mu \pm 2\sigma/\sqrt{N}$")
ax.plot(x_axis, mu_die - 2 * funnel, "k--", lw=1.2)
ax.axhline(mu_die, color="crimson", lw=1.0)
ax.set_xscale("log")
ax.set_xlabel("number of rolls (N)")
ax.set_ylabel("running sample mean")
ax.set_title(f"{n_paths} seeds funnel toward $\\mu = 3.5$")
ax.legend()
fig.tight_layout()
fig.savefig("nb14_wlln_funnel.png", dpi=110)
plt.close(fig)
print(f"Population sd of a fair die: sigma = {sigma_die:.4f}")Why it works: Chebyshev does the heavy lifting
We do not need the full proof. Chebyshev’s inequality says no random variable strays far from its own mean very often: for any random variable Y with mean \mu_Y, variance \sigma_Y^2, and any \epsilon > 0,
\Pr\!\big(|Y - \mu_Y| \ge \epsilon\big) \le \frac{\sigma_Y^2}{\epsilon^2}.
Apply it to Y = \bar{x}. From Ch 1.3 we know \operatorname{Var}(\bar{x}) = \sigma^2/N, so
\Pr\!\big(|\bar{x} - \mu| \ge \epsilon\big) \le \frac{\sigma^2}{N\,\epsilon^2} \xrightarrow[N\to\infty]{} 0.
The probability of a meaningful miss is squeezed to zero. Let us check the bound empirically: simulate many sample means at each N, measure the true miss-probability, and confirm it sits under the Chebyshev ceiling.
eps = 0.3 # tolerance: how far from mu counts as a "miss"
reps_cheb = 20_000
Ns = [10, 30, 100, 300, 1000]
print(f"{'N':>6} {'P(|xbar-mu|>=eps)':>20} {'Chebyshev bound':>18}")
for Nc in Ns:
draws = rng.integers(1, 7, size=(reps_cheb, Nc))
xbar = draws.mean(axis=1)
miss_prob = np.mean(np.abs(xbar - mu_die) >= eps)
bound = sigma_die**2 / (Nc * eps**2)
print(f"{Nc:>6} {miss_prob:>20.5f} {min(bound, 1.0):>18.5f}")2. The Central Limit Theorem
The LLN tells us where \bar{x} ends up. It says nothing about the shape of the cloud of possible \bar{x} values around that destination. That shape is the CLT’s job.
The result in one sentence. Standardize the sample mean — subtract its mean and divide by its standard deviation — and as N grows its distribution becomes the standard normal, no matter what shape the raw population had (as long as it has finite mean and finite variance):
z_N = \frac{\bar{x} - \mu}{\sigma/\sqrt{N}} \xrightarrow{d} N(0, 1).
To be convinced, we start from a population that looks nothing like a bell. Priya studies insurance claims, which are violently right-skewed: most policies cost nothing in a month, a few cost a fortune. A clean stand-in is the exponential distribution with scale 1, so \mu = 1 and \sigma = 1. First, look at the raw population.
mu_exp, sigma_exp = 1.0, 1.0 # exponential(scale=1): mean = sd = scale
raw = rng.exponential(scale=1.0, size=50_000)
fig, ax = plt.subplots()
ax.hist(raw, bins=80, density=True, alpha=0.6, color="steelblue")
xg = np.linspace(0, 8, 300)
ax.plot(xg, stats.expon.pdf(xg), "k-", lw=1.5, label="exponential pdf")
ax.axvline(mu_exp, color="crimson", ls="--", label=r"$\mu = 1$")
ax.set_xlim(0, 8)
ax.set_xlabel("value")
ax.set_ylabel("density")
ax.set_title("Priya's raw population: strongly right-skewed (exponential)")
ax.legend()
fig.tight_layout()
fig.savefig("nb14_exp_population.png", dpi=110)
plt.close(fig)
print(f"raw mean = {raw.mean():.4f}, raw sd = {raw.std():.4f} (both should be ~1)")See it: the bell emerges, panel by panel
Here is the experiment at the heart of this notebook. For a fixed N we draw N exponential values, compute \bar{x}, standardize it to z_N = (\bar{x}-\mu)/(\sigma/\sqrt{N}), and store that one number. Repeat 50{,}000 times to map the sampling distribution of z_N, then histogram it against the standard normal pdf. Finally redo the whole thing for N = 1, 2, 5, 30, 100 and watch the histogram morph from the skewed exponential into a perfect bell.
We lay the five N values out as separate panels (subplots) rather than a true animation, so the notebook renders headless.
def standardized_means(sampler, mu, sigma, N, reps):
# Draw `reps` samples of size N, return the standardized sample means.
samples = sampler(size=(reps, N))
xbar = samples.mean(axis=1)
return (xbar - mu) / (sigma / np.sqrt(N))
reps = 50_000
grid = np.linspace(-4, 4, 300)
normal_pdf = stats.norm.pdf(grid)
Ns_clt = [1, 2, 5, 30, 100]
fig, axes = plt.subplots(1, len(Ns_clt), figsize=(18, 3.6), sharey=True)
for ax, N in zip(axes, Ns_clt):
z = standardized_means(lambda size: rng.exponential(scale=1.0, size=size),
mu_exp, sigma_exp, N, reps)
ax.hist(z, bins=70, range=(-4, 4), density=True, alpha=0.55, color="steelblue")
ax.plot(grid, normal_pdf, "k--", lw=1.5)
ax.set_title(f"N = {N}")
ax.set_xlabel(r"$(\bar{x}-\mu)/(\sigma/\sqrt{N})$")
axes[0].set_ylabel("density")
fig.suptitle("CLT: standardized mean of an exponential population vs. N(0,1) (dashed)", y=1.03)
fig.tight_layout()
fig.savefig("nb14_clt_exponential.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("Exponential CLT panels saved.")Read the panels left to right — this sequence of pictures is the CLT:
- N = 1. A single standardized exponential draw: a sharp wall on the left, a long right tail. Nothing bell-like.
- N = 2. Averaging two draws softens the wall and pulls the peak toward center, but it is still visibly lopsided.
- N = 5. The skew is fading; a bell is trying to form.
- N = 30. Histogram and dashed normal are nearly on top of each other. This is the origin of the folk rule “N \ge 30 is enough.”
- N = 100. Indistinguishable from the normal curve by eye.
We can quantify how fast the bell arrives with skewness, which should decay toward 0 as the standardized mean normalizes. Theory predicts the skewness of the sample mean falls like \text{skew}(x)/\sqrt{N}; for the exponential, \text{skew}(x) = 2.
print(f"{'N':>6} {'empirical skew':>16} {'theory 2/sqrt(N)':>18}")
for N in Ns_clt:
z = standardized_means(lambda size: rng.exponential(scale=1.0, size=size),
mu_exp, sigma_exp, N, reps)
print(f"{N:>6} {stats.skew(z):>16.4f} {2/np.sqrt(N):>18.4f}")3. When the CLT is slow — or fails entirely
Every assumption in the CLT is load-bearing. The one that breaks most dramatically in finance is the innocent clause “finite variance.” Devon pulls daily crypto returns and reasons: “Returns are returns; the CLT kicks in by N = 30, so my sample mean is approximately normal.” He is about to discover that fat tails wreck this.
A fat-tailed distribution makes extreme values far more probable than a normal of the same width allows. A clean knob for fat tails is the Student’s t with \nu degrees of freedom, which finance routinely uses for returns. There is a hierarchy of bad behavior:
- \nu = 5 — finite variance, so the CLT eventually holds, but slowly. The true sd is \sigma = \sqrt{\nu/(\nu-2)}.
- \nu = 2 — the mean exists but the variance is infinite. The CLT’s hypothesis is violated; there is no finite \sigma to standardize by.
- \nu = 1 — the Cauchy distribution. Even the mean fails to exist; the sample mean of N draws has the same distribution as a single draw. The LLN itself fails.
Let us watch the bell refuse to form. We rerun the Section 2 experiment with \nu = 5 at N = 30 and N = 100 — the sizes where the exponential had already converged.
nu5 = 5
sigma_t5 = np.sqrt(nu5 / (nu5 - 2)) # true sd of t_5 = sqrt(5/3)
Ns_fat = [30, 100]
fig, axes = plt.subplots(1, len(Ns_fat), figsize=(13, 4.2), sharey=True)
for ax, N in zip(axes, Ns_fat):
z = standardized_means(lambda size: rng.standard_t(df=nu5, size=size),
0.0, sigma_t5, N, reps)
ax.hist(z, bins=80, range=(-4, 4), density=True, alpha=0.55, color="darkorange")
ax.plot(grid, normal_pdf, "k--", lw=1.5)
ax.set_title(f"Student-t(nu=5), N = {N}")
ax.set_xlabel(r"$(\bar{x}-\mu)/(\sigma/\sqrt{N})$")
axes[0].set_ylabel("density")
fig.suptitle("Slow CLT: heavy tails keep too much probability past $\\pm 3$", y=1.0)
fig.tight_layout()
fig.savefig("nb14_t5_slow.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print(f"t_5 true sd = {sigma_t5:.4f}")Compare this to Priya’s panels. With the exponential, N = 30 already matched the bell. Here, at N = 30, the histogram is too peaked in the middle and too heavy in the tails — the occasional monster draw keeps throwing \bar{x} far out, and even at N = 100 it is visibly off. The only thing that changed is the weight in the tails.
We can put a number on “too heavy in the tails” with the fraction of standardized means landing beyond \pm 3 — for a true N(0,1) that is only about 0.27\%.
def tail_frac(z):
return np.mean(np.abs(z) > 3)
print(f"Normal N(0,1) reference: P(|Z|>3) = {2*stats.norm.sf(3):.5f}\n")
print(f"{'N':>6} {'exp tail>3':>14} {'t5 tail>3':>14}")
for N in Ns_fat:
z_exp = standardized_means(lambda size: rng.exponential(scale=1.0, size=size),
mu_exp, sigma_exp, N, reps)
z_t5 = standardized_means(lambda size: rng.standard_t(df=nu5, size=size),
0.0, sigma_t5, N, reps)
print(f"{N:>6} {tail_frac(z_exp):>14.5f} {tail_frac(z_t5):>14.5f}")The variance disappears (\nu = 2), then the mean disappears (Cauchy, \nu = 1)
At \nu = 5 the CLT is merely slow. Push to \nu = 2 and the variance is infinite: there is no finite \sigma, so the standardization itself is ill-defined and the histogram of sample means never settles onto the bell — it stays heavy-tailed however large N gets. We standardize by the sample spread just to draw something, and watch it stay non-normal.
nu2 = 2
fig, axes = plt.subplots(1, len(Ns_fat), figsize=(13, 4.2), sharey=True)
for ax, N in zip(axes, Ns_fat):
samples = rng.standard_t(df=nu2, size=(reps, N))
xbar = samples.mean(axis=1)
# No finite population sigma exists; standardize by the empirical spread of xbar
# just so the histogram is on a comparable scale. The shape is the point.
z = (xbar - np.median(xbar)) / xbar.std()
ax.hist(z, bins=120, range=(-6, 6), density=True, alpha=0.55, color="firebrick")
ax.plot(grid, normal_pdf, "k--", lw=1.5)
ax.set_title(f"Student-t(nu=2), infinite variance, N = {N}")
ax.set_xlabel("empirically scaled mean")
axes[0].set_ylabel("density")
fig.suptitle("CLT fails: infinite variance keeps the tails fat at every N", y=1.0)
fig.tight_layout()
fig.savefig("nb14_t2_infvar.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("nu=2 panels saved. Tails stay heavy no matter how large N is.")The pathological extreme is \nu = 1, the Cauchy distribution, where even the mean fails to exist. The most vivid way to see the LLN die is the running-mean plot from Section 1: instead of settling onto a line, the running average of Cauchy draws jumps around forever, leaping to wild new levels whenever a giant draw arrives. We put a well-behaved die path next to a Cauchy path for contrast.
Ncau = 5000
xc = np.arange(1, Ncau + 1)
die_path = np.cumsum(rng.integers(1, 7, size=Ncau)) / xc
cauchy_draws = rng.standard_cauchy(size=Ncau)
cauchy_path = np.cumsum(cauchy_draws) / xc
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 4.5))
ax1.plot(die_path, lw=1.0, color="steelblue")
ax1.axhline(3.5, color="crimson", ls="--")
ax1.set_xscale("log")
ax1.set_title("Fair die: running mean settles onto $\\mu = 3.5$ (LLN works)")
ax1.set_xlabel("N"); ax1.set_ylabel("running mean")
ax2.plot(cauchy_path, lw=1.0, color="firebrick")
ax2.axhline(0.0, color="crimson", ls="--")
ax2.set_xscale("log")
ax2.set_title("Cauchy ($\\nu=1$): running mean never settles (LLN fails)")
ax2.set_xlabel("N"); ax2.set_ylabel("running mean")
fig.tight_layout()
fig.savefig("nb14_cauchy_lln_fail.png", dpi=110)
plt.close(fig)
print(f"die final running mean = {die_path[-1]:.4f} (-> 3.5)")
print(f"cauchy final running mean = {cauchy_path[-1]:.4f} (meaningless: it keeps jumping)")
print(f"cauchy max |running mean| over the path = {np.max(np.abs(cauchy_path)):.2f}")4. What Devon should take away
The point is not to give up; it is to be honest about what the tools promise. “The distribution has a mean” buys you, at most, the hope that the LLN works. It does not buy you a fast CLT, and if the variance is infinite it does not buy you the CLT at all. In practice: plot the raw returns and look at the tails before trusting any average; quantify fat-tailedness rather than assuming normality; prefer robust summaries (the median is not hostage to one monster day); demand far larger samples before trusting a normal approximation; and use standard errors that do not silently assume thin tails.
The cheapest, most valuable habit: before trusting an average, plot the data and look at the tails.
Your Turn
How large must N be? The exponential needed about N = 30. Write a loop that, for a given population, finds the smallest N at which the empirical skewness of the standardized mean drops below, say, 0.1. Run it for the exponential, then for a more sharply skewed population (e.g.
rng.gamma(shape=0.2, scale=1.0, size=...), whose skewness is 2/\sqrt{0.2} \approx 4.5). How much bigger does N have to be?Race the tails. Run the exponential and the Student-t(\nu = 5) experiments side by side and find, for each, the N at which the fraction of standardized means beyond \pm 3 first falls to the normal value (\approx 0.0027). Report the ratio — your personal, simulated price of heavy tails.
Test a Pareto tail. The Pareto distribution with tail index \alpha has finite variance only when \alpha > 2, finite mean only when \alpha > 1. Use
(rng.pareto(a=alpha, size=...) + 1)and repeat the CLT panel experiment for \alpha = 3 (finite variance — bell forms), \alpha = 1.5 (infinite variance — bell refuses), and \alpha = 0.8 (infinite mean — even the running mean refuses to settle). Predict each outcome before you plot, then check.
Hint: you already have standardized_means(sampler, mu, sigma, N, reps) — pass a new sampler lambda and the matching \mu, \sigma.