⚙️ Setup¶
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
BLUE="#2563eb"; DEEP="#1d4ed8"; LIGHT="#60a5fa"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
"axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
"axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
rng = np.random.default_rng(70)
sigma_true = 12.0; n = 45
sample = rng.normal(100, sigma_true, n)
xbar = sample.mean(); se = sigma_true/np.sqrt(n); z = stats.norm.ppf(0.975)
print(f"x-bar = {xbar:.2f}, sigma known = {sigma_true}, SE = {se:.3f}")
print(f"95% z-interval: [{xbar-z*se:.2f}, {xbar+z*se:.2f}] (z* = {z:.3f})")
x-bar = 101.61, sigma known = 12.0, SE = 1.789 95% z-interval: [98.11, 105.12] (z* = 1.960)
Clean, but it assumes we know σ. In practice we only have the sample standard deviation s, and plugging s into the z-interval understates the uncertainty for small samples. The fix is the t-distribution.
s = sample.std(ddof=1)
se_t = s/np.sqrt(n)
tcrit = stats.t.ppf(0.975, df=n-1)
print(f"sample sd s = {s:.2f}, SE = s/sqrt(n) = {se_t:.3f}")
print(f"t* (df={n-1}) = {tcrit:.3f} vs z* = {z:.3f}")
print(f"95% t-interval: [{xbar-tcrit*se_t:.2f}, {xbar+tcrit*se_t:.2f}]")
# scipy one-liner agrees
ci = stats.t.interval(0.95, df=n-1, loc=xbar, scale=se_t)
print(f"scipy stats.t.interval: [{ci[0]:.2f}, {ci[1]:.2f}]")
sample sd s = 10.90, SE = s/sqrt(n) = 1.624 t* (df=44) = 2.015 vs z* = 1.960 95% t-interval: [98.34, 104.89] scipy stats.t.interval: [98.34, 104.89]
dfs = [2,5,10,30,100,1000]
print(f"{'df':>5} | {'t*':>7} | {'z*':>6}")
for d in dfs: print(f"{d:>5} | {stats.t.ppf(0.975,d):>7.3f} | {z:>6.3f}")
print("t* shrinks toward z* as df grows: by df=1000 they are nearly identical")
df | t* | z*
2 | 4.303 | 1.960
5 | 2.571 | 1.960
10 | 2.228 | 1.960
30 | 2.042 | 1.960
100 | 1.984 | 1.960
1000 | 1.962 | 1.960
t* shrinks toward z* as df grows: by df=1000 they are nearly identical
The t critical value exceeds 1.96 for small samples (2.776 at df=4) and converges to z as the sample grows. By n ≈ 30–60 the difference is tiny. Using the t-interval keeps small-sample conclusions honest; using z when σ is unknown overstates precision.
pop = rng.gamma(4, 6, 1_000_000) # right-skewed population
MU = pop.mean()
def covers(n):
s = rng.choice(pop, n, replace=False); xb=s.mean(); se=s.std(ddof=1)/np.sqrt(n)
lo,hi = stats.t.interval(0.95, n-1, xb, se); return lo<=MU<=hi
for n in [10, 30, 100]:
cov = np.mean([covers(n) for _ in range(4000)])
print(f"n={n:>4}: t-interval coverage = {cov*100:.1f}% (target 95%)")
print("\nskewed data + tiny n undercovers; by n=30+ the CLT makes the t-interval reliable")
n= 10: t-interval coverage = 93.5% (target 95%)
n= 30: t-interval coverage = 94.0% (target 95%)
n= 100: t-interval coverage = 94.7% (target 95%) skewed data + tiny n undercovers; by n=30+ the CLT makes the t-interval reliable
With a skewed population and n = 10 the interval under-covers (the CLT has not kicked in), but by n = 30 coverage is close to 95% and by n = 100 it is on target. The rule of thumb: the t-interval is trustworthy when the data is roughly symmetric, or when n is at least ~30.
A logistics team exports 180 completed shipments (confidence-intervals-for-a-mean--delivery_times.xlsx) and wants a confidence interval for the average door-to-door delivery time. Sigma is unknown, so we use the t-interval, after checking the shape.
try: ship = pd.read_excel("../../data/confidence-intervals-for-a-mean--delivery_times.xlsx", sheet_name="Shipments")
except FileNotFoundError: ship = pd.read_excel(BASE+"confidence-intervals-for-a-mean--delivery_times.xlsx", sheet_name="Shipments")
print("loaded:", ship.shape)
h = ship["delivery_hours"]
n = len(h); xbar = h.mean(); s = h.std(ddof=1); se = s/np.sqrt(n)
tcrit = stats.t.ppf(0.975, n-1)
lo, hi = stats.t.interval(0.95, n-1, xbar, se)
print(f"n = {n} shipments")
print(f"mean delivery time = {xbar:.2f} h, sd = {s:.2f}, SE = {se:.3f}")
print(f"t* (df={n-1}) = {tcrit:.3f}")
print(f"95% t-interval for the MEAN: [{lo:.2f} h, {hi:.2f} h]")
loaded: (180, 6) n = 180 shipments mean delivery time = 35.85 h, sd = 10.11, SE = 0.753 t* (df=179) = 1.973 95% t-interval for the MEAN: [34.37 h, 37.34 h]
fig,ax=plt.subplots(1,2,figsize=(11,3.2))
ax[0].hist(h, bins=24, color=LIGHT, alpha=0.85); ax[0].axvline(xbar,color=DEEP,lw=2)
ax[0].set_title("Roughly symmetric, so CLT-safe"); ax[0].set_xlabel("hours")
ax[1].errorbar([0],[xbar],yerr=[[xbar-lo],[hi-xbar]],fmt="o",color=BLUE,capsize=8,lw=2,ms=8)
ax[1].set_xlim(-1,1); ax[1].set_xticks([]); ax[1].set_ylabel("mean delivery time (h)")
ax[1].set_title(f"95% CI: [{lo:.1f}, {hi:.1f}] h")
plt.tight_layout(); plt.show()
The average delivery time is about 35.9 hours, with a 95% interval of roughly 34.4 to 37.3 hours. The histogram is reasonably symmetric and n = 180 is large, so the conditions hold and the interval is trustworthy. The team can promise "about a day and a half on average" with quantified confidence, not a bare point.