import numpy as np
import matplotlib
matplotlib.use("Agg") # non-interactive backend: render to memory, never to a window
import matplotlib.pyplot as plt
rng = np.random.default_rng(42) # one generator, one seed -> one reproducible story
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3nb1.3 — The Sampling-Distribution Explorer
Companion to Chapter 1.3 — Estimators and Their Sampling Distributions.
Priya, interning at a regional insurer, was handed one spreadsheet of N=200 wildfire claims and asked for the average payout. She averaged them, got $8,410, and wrote it on a sticky note. The honest question the chapter asked was: how much should she trust that number? If the universe had handed her a different 200 claims, how far would the sticky note have moved?
In real life Priya draws once. She can never see the cloud of values her estimate could have taken. But in a simulation we are gods: we own the true population, so we can re-run her season ten thousand times, collect ten thousand sticky notes, and look at the whole cloud. That cloud is the sampling distribution, and this notebook builds it by brute force.
Here is the plan, and you can already guess the punchlines from the chapter:
- Build Priya’s true population (lognormal losses, \mu=\$8{,}000) and take her one sample.
- Re-draw the season 10,000 times and watch the sampling distribution of \bar{x} appear.
- Confirm it centers on \mu (unbiased) and its spread is \sigma/\sqrt N — measured, not just derived.
- Overlay several N to see the 1/\sqrt N shrink.
- Score the plain mean against a trimmed mean using MSE = bias² + variance.
- Your Turn: switch on a heavy tail and watch the tidy story wobble.
0. Setup
We fix the random seed so every run of this notebook tells the same story — reproducibility is non-negotiable. We also force matplotlib’s non-interactive Agg backend so the notebook runs cleanly start-to-finish without a display.
1. Priya’s “true” population — and her one sample
Insurance losses are the textbook case of a right-skewed distribution: a flood of small claims and a thin tail of catastrophic ones. The standard rough model is the lognormal — a variable whose logarithm is normal.
A lognormal with log-mean m and log-sd s has actual mean \mathbb{E}[x]=\exp(m+s^2/2). We want the population mean to be exactly \mu=\$8{,}000, so we pick the skewness knob s=0.9 and then solve for m:
m=\ln(8000)-\tfrac{s^2}{2}\quad\Longrightarrow\quad \exp(m+\tfrac{s^2}{2})=8000\ \text{exactly}.
Because we built the population, we know \mu and \sigma in closed form — the luxury Priya never has.
# --- population knobs ---
sigma_log = 0.9 # skewness dial for the lognormal
mu_log = np.log(8000) - sigma_log**2 / 2 # solved so the mean is exactly 8000
# closed-form population mean and sd of a lognormal (these are the "truth" we get to see)
mu_true = np.exp(mu_log + sigma_log**2 / 2)
sigma_true = np.sqrt((np.exp(sigma_log**2) - 1) * np.exp(2 * mu_log + sigma_log**2))
def draw_claims(n):
# n i.i.d. claim amounts from Priya's lognormal population
return rng.lognormal(mean=mu_log, sigma=sigma_log, size=n)
print(f"True population mean mu = ${mu_true:,.2f}")
print(f"True population sd sigma = ${sigma_true:,.2f}")Now Priya’s reality: one draw of 200 claims, one estimate. Run this cell and note the number. That single \bar{x} is all she will ever get to see — everything below is the part of the universe hidden from her.
N_PRIYA = 200
priyas_sample = draw_claims(N_PRIYA)
print(f"Priya's single estimate: x-bar = ${priyas_sample.mean():,.2f}")
print(f"(the true mu she is chasing: ${mu_true:,.2f})")
print(f"\nA peek at her first 8 claims (note the skew — small claims, occasional big one):")
print(np.round(np.sort(priyas_sample)[-8:], 0)) # the 8 largest, sorted2. Playing god: 10,000 alternate seasons
Here is the conceptual core of the chapter, turned into three lines of code. We re-run Priya’s season n_reps times. Each repetition draws a fresh sample of size N, computes its \bar{x}, records it, and throws the sample away. The array of recorded \bar{x} values is the sampling distribution — we are watching the random variable \bar{x} do its thing over and over, instead of seeing it just once.
def sampling_distribution(n, n_reps=10_000, estimator=np.mean):
# n_reps estimates, each from its own fresh sample of size n.
# estimator: a function array -> scalar (default the sample mean).
return np.array([estimator(draw_claims(n)) for _ in range(n_reps)])
xbars_200 = sampling_distribution(200)
print(f"Collected {xbars_200.size:,} alternate-universe estimates for N=200.")3. The two summaries: centered on \mu, spread \sigma/\sqrt N
The chapter derived two facts about this cloud before ever drawing it:
- Unbiased: \mathbb{E}[\bar{x}]=\mu. So the mean of our 10,000 estimates should sit on \mu.
- Spread: \operatorname{Var}(\bar{x})=\sigma^2/N, i.e. the standard error is \sigma/\sqrt N.
Let’s check both against the simulation. This is the moment the formulas stop being algebra and become something we measured.
se_theory = sigma_true / np.sqrt(200) # sigma / sqrt(N)
var_theory = sigma_true**2 / 200 # sigma^2 / N
print("CENTER (unbiasedness check)")
print(f" true mu .......................... ${mu_true:,.2f}")
print(f" mean of 10,000 x-bars ............ ${xbars_200.mean():,.2f} <- lands on mu")
print()
print("SPREAD (Var(x-bar) = sigma^2 / N check)")
print(f" empirical Var(x-bar) ............. {xbars_200.var():,.0f}")
print(f" theory sigma^2 / N .............. {var_theory:,.0f}")
print(f" empirical se (std of x-bars) ..... ${xbars_200.std():,.2f}")
print(f" theory sigma / sqrt(N) .......... ${se_theory:,.2f}")
print()
ratio = xbars_200.var() / var_theory
print(f" ratio empirical/theory variance .. {ratio:.3f} (want ~1.0)")
assert 0.9 < ratio < 1.1, "empirical variance should match sigma^2/N"
print(" PASS: empirical Var(x-bar) matches sigma^2/N.")Two lines of that output deserve applause. The mean of the cloud sits essentially on $8,000 — unbiasedness, seen. And the empirical variance of the 10,000 estimates matches \sigma^2/N to within a percent or two — the \sigma/\sqrt N formula, measured. The formula and the histogram are two views of the same object.
Now the part you cannot get from a formula: the shape. Let’s overlay where Priya’s lone estimate fell against the full cloud she never got to see.
fig, ax = plt.subplots()
ax.hist(xbars_200, bins=50, density=True, color="#4C72B0", edgecolor="white")
ax.axvline(mu_true, color="black", lw=2, label=f"true mu = ${mu_true:,.0f}")
ax.axvline(priyas_sample.mean(), color="crimson", lw=2, ls="--",
label=f"Priya's one estimate = ${priyas_sample.mean():,.0f}")
ax.set_xlabel("estimate x-bar ($)")
ax.set_ylabel("density")
ax.set_title("Sampling distribution of x-bar (N=200): the cloud Priya never sees")
ax.legend()
fig.tight_layout()
fig.savefig("nb13_fig1_one_sample.png", dpi=90)
plt.close(fig)
print("saved nb13_fig1_one_sample.png")
print("Priya's red dashed line is just ONE draw from the blue cloud -- a typical, unremarkable member.")4. Watching 1/\sqrt N with your own eyes
The spread shrinks like 1/\sqrt N — not 1/N. That square root is the most consequential in statistics, and it is a little disappointing: to halve your standard error you must quadruple your sample. Let’s overlay the sampling distribution for several sample sizes and watch the cloud collapse onto \mu.
Ns = [10, 50, 200, 1000]
colors = ["#C44E52", "#DD8452", "#4C72B0", "#55A868"]
fig, ax = plt.subplots(figsize=(9, 5))
print(f"{'N':>6} | {'empirical se':>14} | {'theory sigma/sqrtN':>20}")
print("-" * 48)
for n, c in zip(Ns, colors):
xb = sampling_distribution(n)
se_emp = xb.std()
se_th = sigma_true / np.sqrt(n)
print(f"{n:>6} | ${se_emp:>12,.0f} | ${se_th:>18,.0f}")
ax.hist(xb, bins=60, density=True, alpha=0.55, color=c,
label=f"N={n} (se=${se_emp:,.0f})")
ax.axvline(mu_true, color="black", lw=2, label=f"mu = ${mu_true:,.0f}")
ax.set_xlim(3000, 16000)
ax.set_xlabel("estimate x-bar ($)")
ax.set_ylabel("density")
ax.set_title("Bigger N -> narrower cloud, same center (the 1/sqrt(N) shrink)")
ax.legend()
fig.tight_layout()
fig.savefig("nb13_fig2_overlay_N.png", dpi=90)
plt.close(fig)
print("\nsaved nb13_fig2_overlay_N.png")Read the table down the rows. Going from N=10 to N=1000 is a 100-fold jump in data, so 1/\sqrt N predicts the standard error should fall by \sqrt{100}=10-fold — and it does. The center never budges, because unbiasedness holds at every N. You are watching \bar{x}\xrightarrow{p}\mu: the cloud being squeezed onto the truth. That squeeze is consistency, and at small N you can also see the histogram is still visibly right-skewed — it remembers the lognormal underneath — while at large N it relaxes into a clean bell. Why the bell? That’s the Central Limit Theorem, the headline of Chapter 1.4.
Let’s confirm the 1/\sqrt N rate numerically rather than by eye.
# If se = sigma / sqrt(N), then doubling-and-doubling N should scale se by 1/sqrt(ratio).
N_a, N_b = 100, 400 # N_b is 4x N_a
se_a = sampling_distribution(N_a).std()
se_b = sampling_distribution(N_b).std()
print(f"se at N={N_a} ............ ${se_a:,.2f}")
print(f"se at N={N_b} ............ ${se_b:,.2f}")
print(f"observed ratio se_a/se_b . {se_a/se_b:.3f}")
print(f"predicted sqrt(N_b/N_a) .. {np.sqrt(N_b/N_a):.3f} (4x the data -> 2x sharper)")
print("\nQuadrupling N roughly halves the standard error -- the price of precision.")5. Two recipes, one scorecard: MSE = bias² + variance
The sample mean is not the only estimator of \mu. The chapter floated a 5%-trimmed mean: drop the largest and smallest 5% of claims, then average. On a right-skewed loss distribution this lops off the volatile catastrophe tail — which should cut the variance. But it also systematically discards the heavy right tail, so it lands low on average — that’s bias.
Which recipe is better? Neither bias alone nor variance alone can say. The honest scorecard is the mean squared error, and it splits exactly into the two costs:
\operatorname{MSE}(\hat\mu)=\mathbb{E}\!\big[(\hat\mu-\mu)^2\big]=\operatorname{Var}(\hat\mu)+\operatorname{Bias}(\hat\mu)^2.
We’ll estimate all three quantities by simulation for each recipe and let the arithmetic decide.
from scipy import stats
def trimmed_mean_5(x):
# mean after dropping the top and bottom 5% of the sample
return stats.trim_mean(x, 0.05)
N = 200
xbar_dist = sampling_distribution(N, estimator=np.mean)
trim_dist = sampling_distribution(N, estimator=trimmed_mean_5)
def score(dist, name):
# estimate bias, variance, MSE of an estimator from its simulated distribution
bias = dist.mean() - mu_true
var = dist.var()
mse = np.mean((dist - mu_true) ** 2) # direct MSE
print(f"{name:<22} bias=${bias:>8,.1f} var={var:>12,.0f} "
f"bias^2={bias**2:>12,.0f} MSE={mse:>12,.0f}")
# verify the decomposition MSE = var + bias^2 holds numerically
assert abs(mse - (var + bias**2)) < 1e-6 * mse, "MSE decomposition must hold"
return mse
print(f"Scoring two estimators of mu at N={N} (truth mu=${mu_true:,.0f}):\n")
mse_mean = score(xbar_dist, "plain mean")
mse_trim = score(trim_dist, "5%-trimmed mean")
print()
print("Check: for each recipe, MSE == var + bias^2 (PASS, asserted above).")
winner = "trimmed mean" if mse_trim < mse_mean else "plain mean"
print(f"Lower MSE wins -> {winner} (MSE {min(mse_mean, mse_trim):,.0f} "
f"vs {max(mse_mean, mse_trim):,.0f}).")Look carefully — the result is more honest than the textbook’s optimistic hypothetical. The plain mean is unbiased (bias ≈ 0), so its entire MSE is variance. The trimmed mean does exactly what we predicted on the variance side: it cuts the variance roughly in half by lopping off the volatile catastrophe tail. But on a tail this heavy, throwing away the right tail also drags its center more than $1,000 below \mu. That bias is large, and MSE charges you its square — so the bias² term blows up and the trimmed mean loses decisively here.
That is not a failure of the experiment; it is the lesson. The chapter said plainly that whether the trimmed mean beats the mean “is an empirical question — exactly what you’ll measure in nb1.3,” and you just measured it: for this lognormal, the plain mean wins. A little bias is a bargain only when the variance it buys outweighs the bias² you pay — and on a strongly skewed loss distribution, trimming for the population mean pays too much bias. (Tellingly, the trimmed mean is a much better recipe if your target is the median or a less skewed population — try sigma_log = 0.3 in your turn and watch the gap close.)
This is the bias–variance tradeoff in the flesh, and it is not a curiosity: it is the engine behind shrinkage estimators, regularization, and essentially all of machine learning (you’ll meet it again in Week 4). Let’s see both recipes’ clouds side by side so the tradeoff is visible, not just tabulated.
fig, ax = plt.subplots(figsize=(9, 5))
ax.hist(xbar_dist, bins=55, density=True, alpha=0.55, color="#4C72B0",
label=f"plain mean (MSE={mse_mean:,.0f})")
ax.hist(trim_dist, bins=55, density=True, alpha=0.55, color="#55A868",
label=f"5%-trimmed (MSE={mse_trim:,.0f})")
ax.axvline(mu_true, color="black", lw=2, label=f"mu = ${mu_true:,.0f}")
ax.axvline(trim_dist.mean(), color="#55A868", lw=2, ls="--",
label=f"trimmed center = ${trim_dist.mean():,.0f} (biased low)")
ax.set_xlabel("estimate ($)")
ax.set_ylabel("density")
ax.set_title("Plain vs trimmed mean: the trimmed cloud is shifted left but much tighter")
ax.legend()
fig.tight_layout()
fig.savefig("nb13_fig3_mean_vs_trimmed.png", dpi=90)
plt.close(fig)
print("saved nb13_fig3_mean_vs_trimmed.png")
print("Green cloud: nudged left of the black truth line (bias) but visibly narrower (less variance).")Your Turn
You now have a working sampling-distribution machine. Turn the dials.
1. Crank up the skew — and foreshadow Ch 1.4. Re-run sections 3–4 with sigma_log = 1.8 (a far heavier right tail; remember to re-solve mu_log so the mean stays $8,000). Does the empirical variance still match \sigma^2/N? At N=10, does the sampling distribution still look skewed? How large must N get before it looks like a bell? You’re stress-testing the Central Limit Theorem before you’ve met it.
2. A genuinely heavy tail (a peek at where formulas break). Replace draw_claims with draws from a Student-t with 2 degrees of freedom (rng.standard_t(2, size=n), shifted to mean $8,000). This distribution has infinite variance. Re-run the \sigma/\sqrt N check — what goes wrong, and why can’t the tidy Var(x-bar)=sigma^2/N story survive? (This is exactly the kind of assumption-failure Chapter 1.4 will care about.)
3. Add a third estimator. Add the median to the section-5 contest (np.median). Where does it land on the bias–variance scorecard relative to the mean and the trimmed mean? On this lognormal, is the median biased for \mu? Should it be?
4. The price of precision (the chapter’s check question). Using \sigma=\$12{,}000, solve for the N that drives the standard error below $100. Confirm your algebra by simulating at that N and reading off the empirical se. Is collecting that many claims realistic in one season?
Whatever you change, keep rng = np.random.default_rng(42) at the top so your experiment stays reproducible — and so a future you can rebuild the exact same alternate universes.