⚙️ Setup¶
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(45)
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False})
VIOLET="#7c3aed"; PINK="#db2777"; TEAL="#0d9488"
print("ready")
ready
normal = rng.normal(0,1,size=200_000)
skewed = rng.exponential(1.0,size=200_000)
for name,d in [("Normal",normal),("Exponential",skewed)]:
print(f"{name:12s}: mean={d.mean():+.3f}, var={d.var():.3f}, skew={stats.skew(d):+.3f}, kurtosis={stats.kurtosis(d):+.3f}")
Normal : mean=+0.000, var=1.000, skew=-0.000, kurtosis=+0.002 Exponential : mean=+0.996, var=0.999, skew=+2.032, kurtosis=+6.351
The normal is symmetric (skew ≈ 0) with the baseline tail weight (excess kurtosis ≈ 0). The exponential is strongly right-skewed (skew ≈ 2) with heavier tails (positive excess kurtosis). The first four moments capture center, spread, lopsidedness, and tail heaviness, the four ways a distribution can differ in shape.
rate = 2.0
X = rng.exponential(1/rate, size=500_000)
def M(t): return np.mean(np.exp(t*X)) # empirical MGF
h = 1e-3
mean_from_mgf = (M(h) - M(-h)) / (2*h) # M'(0)
EX2_from_mgf = (M(h) - 2*M(0) + M(-h)) / h**2 # M''(0)
print(f"M'(0) = {mean_from_mgf:.4f} (true mean 1/rate = {1/rate})")
print(f"M''(0) = {EX2_from_mgf:.4f} (true E[X^2] = 2/rate^2 = {2/rate**2})")
print(f"variance = M''(0) - M'(0)^2 = {EX2_from_mgf - mean_from_mgf**2:.4f} (true 1/rate^2 = {1/rate**2})")
M'(0) = 0.5015 (true mean 1/rate = 0.5) M''(0) = 0.5040 (true E[X^2] = 2/rate^2 = 0.5) variance = M''(0) - M'(0)^2 = 0.2525 (true 1/rate^2 = 0.25)
Differentiating the MGF at 0 recovers the mean (0.5) and the second moment (0.5), and from them the variance (0.25). The MGF is a single function that encodes every moment, and because it uniquely identifies a distribution, it is the classic tool for proving that a sum of independent variables has a particular distribution (their MGFs simply multiply).
X = rng.exponential(1.0, size=500_000) # mean 1, sd 1
mu, sd = X.mean(), X.std()
for k in [2, 3, 4]:
actual = np.mean(np.abs(X-mu) >= k*sd)
print(f"k={k}: P(|X-mu| >= {k} sd) actual = {actual:.4f} Chebyshev bound 1/k^2 = {1/k**2:.4f}")
k=2: P(|X-mu| >= 2 sd) actual = 0.0499 Chebyshev bound 1/k^2 = 0.2500 k=3: P(|X-mu| >= 3 sd) actual = 0.0184 Chebyshev bound 1/k^2 = 0.1111 k=4: P(|X-mu| >= 4 sd) actual = 0.0067 Chebyshev bound 1/k^2 = 0.0625
Chebyshev's bound 1/k² holds for every distribution, here comfortably (the exponential's real tails are much lighter than the worst case). The power of these inequalities is that they need almost no assumptions, just a mean (Markov) or a variance (Chebyshev), making them the universal safety net of probability.
def tail_prob(n, eps, trials=20000):
means = rng.random((trials, n)).mean(axis=1) # uniform[0,1], mu=0.5
return np.mean(np.abs(means - 0.5) >= eps)
eps = 0.1
for n in [25, 100, 400]:
emp = tail_prob(n, eps)
bound = 2*np.exp(-2*n*eps**2)
print(f"n={n:4d}: P(|mean-0.5|>={eps}) empirical = {emp:.4f} Hoeffding bound = {bound:.4f}")
n= 25: P(|mean-0.5|>=0.1) empirical = 0.0824 Hoeffding bound = 1.2131 n= 100: P(|mean-0.5|>=0.1) empirical = 0.0004 Hoeffding bound = 0.2707 n= 400: P(|mean-0.5|>=0.1) empirical = 0.0000 Hoeffding bound = 0.0007
As n grows, the chance the sample mean is more than 0.1 from the truth collapses, and Hoeffding's exponential bound tracks (and safely exceeds) it. This concentration is a sharper, finite-sample cousin of the Law of Large Numbers, and it is the mathematical reason a model evaluated on enough data gives a trustworthy estimate.
def hoeffding_n(eps, delta): return int(np.ceil(np.log(2/delta) / (2*eps**2)))
for eps, delta in [(0.05,0.05),(0.02,0.05),(0.01,0.01)]:
print(f"to be within +/-{eps} with {100*(1-delta):.0f}% confidence: need n >= {hoeffding_n(eps,delta):,} samples")
# method of moments: estimate exponential rate from the sample mean (lambda = 1/mean)
data = rng.exponential(1/3.0, size=2000)
print(f"\nmethod of moments: rate_hat = 1/mean = {1/data.mean():.3f} (true 3.0)")
to be within +/-0.05 with 95% confidence: need n >= 738 samples to be within +/-0.02 with 95% confidence: need n >= 4,612 samples to be within +/-0.01 with 99% confidence: need n >= 26,492 samples method of moments: rate_hat = 1/mean = 3.035 (true 3.0)
Hoeffding inverts into a sample-complexity rule: to pin a test accuracy to within ±0.02 at 95% confidence you need about 4,600 examples, a concrete answer to "how much test data is enough". These bounds are the backbone of PAC learning and generalization theory. The same moment ideas give the method of moments, a quick estimator (here rate = 1/mean), and moment matching is how models like GANs and MMD are trained to mimic a target distribution.