⚙️ Setup¶
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(40)
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False})
TEAL="#0d9488"; PINK="#db2777"; AMBER="#d97706"
print("ready")
ready
# a skewed population: exponential with mean 2, sd 2
pop_mean, pop_sd = 2.0, 2.0
n = 30
means = np.array([rng.exponential(2.0, size=n).mean() for _ in range(50_000)])
print(f"population: mean = {pop_mean:.3f}, sd = {pop_sd:.3f}")
print(f"sample means: mean = {means.mean():.3f}, sd = {means.std():.3f}")
print(f"theory: SE = sd/sqrt(n) = {pop_sd/np.sqrt(n):.3f}")
population: mean = 2.000, sd = 2.000 sample means: mean = 1.997, sd = 0.364 theory: SE = sd/sqrt(n) = 0.365
fig,ax=plt.subplots(figsize=(7,3.2))
xs=np.linspace(0,8,200)
ax.plot(xs, np.exp(-xs/2)/2, color=AMBER, lw=2, label="population (skewed)")
ax.hist(means, bins=60, density=True, color=TEAL, alpha=0.65, label=f"means of n={n}")
ax.axvline(pop_mean, color=PINK, ls="--", lw=2, label="population mean = 2")
ax.set_xlabel("value"); ax.set_ylabel("density"); ax.set_title("Sample means cluster tightly around the population mean"); ax.legend()
plt.tight_layout(); plt.show()
The population is heavily skewed, yet the sample means form a tight, symmetric mound centered exactly on the population mean of 2. The means vary far less than individual values, their spread is the standard error, sd/sqrt(n) = 0.365.
fig,axes=plt.subplots(1,4,figsize=(11,2.8))
for ax,n in zip(axes,[1,2,5,30]):
m=np.array([rng.exponential(2.0,size=n).mean() for _ in range(40_000)])
ax.hist(m,bins=45,density=True,color=TEAL,alpha=0.8)
ax.set_title(f"n = {n}"); ax.set_yticks([]); ax.set_xlim(0,6)
axes[0].set_ylabel("density")
fig.suptitle("Means of an exponential population become normal as n grows", y=1.05)
plt.tight_layout(); plt.show()
At n=1 the distribution of "means" is just the skewed population itself. By n=5 it is already more symmetric, and by n=30 it is a clean bell. This is the Central Limit Theorem: averages of many independent draws are approximately normal, which is why the normal distribution is everywhere and why so much of inference can lean on it.
sigma = 2.0
for n in [25, 100, 400, 1600]:
print(f"n={n:5d}: SE = sigma/sqrt(n) = {sigma/np.sqrt(n):.4f}")
print("\nQuadrupling n (25 -> 100 -> 400 -> 1600) halves the SE each step.")
n= 25: SE = sigma/sqrt(n) = 0.4000 n= 100: SE = sigma/sqrt(n) = 0.2000 n= 400: SE = sigma/sqrt(n) = 0.1000 n= 1600: SE = sigma/sqrt(n) = 0.0500 Quadrupling n (25 -> 100 -> 400 -> 1600) halves the SE each step.
ns=np.arange(10,1001)
fig,ax=plt.subplots(figsize=(7,3))
ax.plot(ns, sigma/np.sqrt(ns), color=TEAL, lw=2.2)
for n in [50,200,800]: ax.scatter(n, sigma/np.sqrt(n), color=PINK, zorder=3)
ax.set_xlabel("sample size n"); ax.set_ylabel("standard error"); ax.set_title("SE shrinks like 1/sqrt(n): diminishing returns")
plt.tight_layout(); plt.show()
The standard error falls steeply at first, then crawls: going from n=50 to n=200 helps a lot, but n=800 to n=950 barely moves it. This diminishing return is the central economic fact of data collection, and the reason "just get more data" eventually stops paying off.
true_mu, sigma, n = 2.0, 2.0, 50
SE = sigma/np.sqrt(n)
covered = 0; trials = 20_000
for _ in range(trials):
xbar = rng.exponential(2.0, size=n).mean()
lo, hi = xbar - 1.96*SE, xbar + 1.96*SE
if lo <= true_mu <= hi: covered += 1
print(f"95% CIs that captured the true mean: {covered/trials*100:.1f}%")
95% CIs that captured the true mean: 95.3%
About 95% of the intervals contain the true mean, exactly as advertised. This is the bridge the CLT builds: it lets us attach a margin of error, x_bar +/- 1.96 x SE, to a single sample, the foundation of the confidence intervals and hypothesis tests in the chapters ahead.
# per-example gradient component: true value g with noise
true_g, noise_sd = 0.50, 3.0
for batch in [8, 32, 128, 512]:
ests = np.array([(true_g + rng.normal(0, noise_sd, size=batch)).mean() for _ in range(5_000)])
print(f"batch={batch:4d}: gradient estimate SE = {ests.std():.4f} (theory {noise_sd/np.sqrt(batch):.4f})")
# a test-accuracy estimate is also a sample mean: SE = sqrt(p(1-p)/n)
p, n_test = 0.92, 2000
print(f"\ntest accuracy 0.92 on {n_test} examples: SE = {np.sqrt(p*(1-p)/n_test):.4f}, margin +/- {1.96*np.sqrt(p*(1-p)/n_test):.3f}")
batch= 8: gradient estimate SE = 1.0643 (theory 1.0607) batch= 32: gradient estimate SE = 0.5299 (theory 0.5303) batch= 128: gradient estimate SE = 0.2674 (theory 0.2652) batch= 512: gradient estimate SE = 0.1335 (theory 0.1326) test accuracy 0.92 on 2000 examples: SE = 0.0061, margin +/- 0.012
A mini-batch gradient is a noisy estimate of the true gradient, and its noise shrinks like 1/sqrt(batch), pure standard error. That is why large batches give smoother updates (and why the noise of small batches can actually help generalization). The same formula puts an honest margin of error on a reported accuracy: 0.92 on 2,000 examples is really 0.92 +/- 0.012.