import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
AMBER="#d97706"; TEAL="#0d9488"; INK="#1a2138"; GRID="#e6e9f2"; PINK="#db2777"
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/"
try:    df = pd.read_csv("../../data/logistics_freight_operations.csv")
except FileNotFoundError: df = pd.read_csv(BASE+"logistics_freight_operations.csv")
print("loaded:", df.shape)
df.head()

loaded: (1000, 4)

w = df["individual_package_weight_lbs"]
print(f"sample of weights: n={len(w)}, mean={w.mean():.2f}, sd={w.std():.2f}, skew={stats.skew(w):.3f} (right-skewed)")
# fit a gamma by method of moments and generate a large population
a = w.mean()**2/w.var(); b = w.var()/w.mean()
rng = np.random.default_rng(57)
POP = rng.gamma(a, b, size=1_000_000)
MU, SIG = POP.mean(), POP.std()
print(f"generated population: N=1,000,000, mu={MU:.3f}, sigma={SIG:.3f}, skew={stats.skew(POP):.3f}")

sample of weights: n=1000, mean=15.07, sd=4.00, skew=0.527 (right-skewed)
generated population: N=1,000,000, mu=15.073, sigma=4.001, skew=0.532

fig,ax=plt.subplots(figsize=(7,3.2))
ax.hist(POP,bins=80,density=True,color=AMBER,alpha=0.6,label="population (1,000,000 packages)")
ax.hist(w,bins=40,density=True,histtype="step",color=TEAL,lw=2,label="original 1,000-package sample")
ax.axvline(MU,color=PINK,ls="--",lw=2,label=f"population mean {MU:.1f} lb")
ax.set_xlabel("package weight (lb)"); ax.set_ylabel("density"); ax.set_title("A skewed population, grown from the data"); ax.legend()
plt.tight_layout(); plt.show()

def sample_means(n, reps=40_000): return POP[rng.integers(0,len(POP),size=(reps,n))].mean(axis=1)
fig,axes=plt.subplots(1,4,figsize=(12,2.8))
for ax,n in zip(axes,[1,5,30,100]):
    means=sample_means(n)
    ax.hist(means,bins=45,density=True,color=AMBER,alpha=0.75)
    xs=np.linspace(means.min(),means.max(),200); ax.plot(xs,stats.norm.pdf(xs,MU,SIG/np.sqrt(n)),color=TEAL,lw=2)
    ax.set_title(f"n = {n}\nSE = {SIG/np.sqrt(n):.2f}"); ax.set_yticks([])
fig.suptitle("Distribution of the sample mean: skewed at n=1, normal by n=30", y=1.06, fontweight="bold")
plt.tight_layout(); plt.show()

for n in [1,5,30,100]:
    means=sample_means(n)
    print(f"n={n:3d}: mean of means {means.mean():.3f} (mu {MU:.3f}), SE {means.std():.3f} (sigma/sqrt(n) {SIG/np.sqrt(n):.3f}), skew {stats.skew(means):+.3f}")

n=  1: mean of means 15.075 (mu 15.073), SE 3.977 (sigma/sqrt(n) 4.001), skew +0.501
n=  5: mean of means 15.063 (mu 15.073), SE 1.782 (sigma/sqrt(n) 1.789), skew +0.241
n= 30: mean of means 15.070 (mu 15.073), SE 0.730 (sigma/sqrt(n) 0.730), skew +0.086
n=100: mean of means 15.073 (mu 15.073), SE 0.400 (sigma/sqrt(n) 0.400), skew +0.059

n = 50; SE = SIG/np.sqrt(n)
truck = rng.choice(POP, n, replace=False)
xbar = truck.mean()
lo, hi = xbar - 1.96*SE, xbar + 1.96*SE
print(f"one truckload of {n}: sample mean = {xbar:.2f} lb")
print(f"95% CI for fleet mean weight = [{lo:.2f}, {hi:.2f}]  (true mu = {MU:.2f})")
# coverage check
cov=np.mean([(lambda m: m-1.96*SE<=MU<=m+1.96*SE)(rng.choice(POP,n,replace=False).mean()) for _ in range(5000)])
print(f"coverage of 95% CIs over 5000 truckloads = {cov*100:.1f}%")

one truckload of 50: sample mean = 14.52 lb
95% CI for fleet mean weight = [13.41, 15.63]  (true mu = 15.07)
coverage of 95% CIs over 5000 truckloads = 94.8%

	truck_id	route_id	driver_id	individual_package_weight_lbs
0	TRK_4000	RT_I10	DRV_124	19.73
1	TRK_4001	RT_I95	DRV_127	15.54
2	TRK_4002	RT_I95	DRV_147	10.08
3	TRK_4003	RT_I80	DRV_148	8.69
4	TRK_4004	RT_I80	DRV_125	17.86

⚙️ Setup & data¶