⚙️ Setup¶
import numpy as np
rng = np.random.default_rng(461)
print("ready")
ready
ests = np.array([rng.normal(50,10,size=20).mean() for _ in range(100_000)])
print(f"E[sample mean] = {ests.mean():.3f} (true 50)")
print(f"bias = {ests.mean()-50:+.3f} -> unbiased")
E[sample mean] = 49.988 (true 50) bias = -0.012 -> unbiased
Answer: the average of the estimates is 50, so the bias is essentially 0: the sample mean is unbiased. On average it neither over- nor under-shoots the true mean, whatever the sample size.
n=4
var_n = np.mean([rng.normal(0,3,size=n).var(ddof=0) for _ in range(200_000)])
var_nm1 = np.mean([rng.normal(0,3,size=n).var(ddof=1) for _ in range(200_000)])
print(f"true variance = 9")
print(f"E[/n] = {var_n:.3f} (biased low; factor (n-1)/n = {(n-1)/n})")
print(f"E[/(n-1)]= {var_nm1:.3f} (unbiased)")
true variance = 9 E[/n] = 6.757 (biased low; factor (n-1)/n = 0.75) E[/(n-1)]= 8.979 (unbiased)
Answer: dividing by n averages to about 6.75 (= 9 × 3/4), biased low, while dividing by n − 1 averages to 9, unbiased. Bessel's correction (n − 1) compensates for using the sample mean in place of the unknown true mean.
data = rng.exponential(1/0.4, size=5000)
mle = 1/data.mean()
print(f"MLE rate = 1/mean = {mle:.3f} (true 0.4)")
MLE rate = 1/mean = 0.399 (true 0.4)
Answer: the exponential MLE is rate = 1/mean ≈ 0.4. Maximizing the log-likelihood n·log(λ) − λ·Σx gives λ = 1/x̄ exactly, a textbook maximum-likelihood result.
for n in [10, 100, 1000, 10000]:
ests = [rng.normal(5,2,size=n).mean() for _ in range(2000)]
print(f"n={n:6d}: sd of estimate = {np.std(ests):.4f} (theory 2/sqrt(n) = {2/np.sqrt(n):.4f})")
n= 10: sd of estimate = 0.6429 (theory 2/sqrt(n) = 0.6325) n= 100: sd of estimate = 0.2019 (theory 2/sqrt(n) = 0.2000) n= 1000: sd of estimate = 0.0632 (theory 2/sqrt(n) = 0.0632) n= 10000: sd of estimate = 0.0190 (theory 2/sqrt(n) = 0.0200)
Answer: the spread of the estimate shrinks like 2/√n, from about 0.63 at n = 10 to 0.02 at n = 10,000, converging on the true mean. This is consistency: with enough data the estimator is essentially certain to be close.
true_mu = 20
ests = np.array([0.9*rng.normal(true_mu,4,size=50).mean() for _ in range(100_000)])
bias = ests.mean() - true_mu
var = ests.var()
mse = ((ests - true_mu)**2).mean()
print(f"bias = {bias:.3f}, bias^2 = {bias**2:.3f}")
print(f"variance = {var:.3f}")
print(f"bias^2 + variance = {bias**2+var:.3f} vs MSE = {mse:.3f}")
bias = -2.001, bias^2 = 4.005 variance = 0.261 bias^2 + variance = 4.266 vs MSE = 4.266
Answer: shrinking the mean by 10% introduces a bias of about −2 (bias² ≈ 4), and the identity holds: MSE = bias² + variance. This decomposition is the exact form of the bias-variance trade-off, sometimes a little bias buys a large drop in variance and a lower overall MSE, which is the whole rationale for regularization.