import subprocess, sys
try:
    import pingouin as pg
except ImportError:   # Colab does not ship pingouin, and 0.7 needs a newer pandas than Colab has
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "pingouin<0.7"], check=False)
    import pingouin as pg

import numpy as np, pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False,
    "axes.grid":True,"grid.alpha":0.22,"axes.titleweight":"bold","axes.titlesize":12.5,"axes.titlelocation":"left"})
GOLD, AMBER, DARK, RD, GR, MUT = "#b45309", "#d97706", "#78350f", "#dc2626", "#16a34a", "#94a3b8"
TARGET = 340.0   # grams printed on the label

BASE = "https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
fn = "capstone-coffee-fill-weight.xlsx"
try: raw = pd.read_excel("../../data/" + fn, sheet_name="fills")
except FileNotFoundError: raw = pd.read_excel(BASE + fn, sheet_name="fills")
print("loaded", raw.shape[0], "rows")
raw.head()

loaded 78 rows

print(raw.dtypes)
print()
print("a continuous outcome (fill_g) vs one target value  ->  one-sample t-test family")

bag_id     object
line       object
shift      object
fill_g    float64
dtype: object

a continuous outcome (fill_g) vs one target value  ->  one-sample t-test family

s = raw.fill_g
print(f"rows {len(raw)}   unique bag_id {raw.bag_id.nunique()}   "
      f"fill_g recorded on {s.notna().sum()}")
print(f"fill_g, raw and uncleaned:  median {s.median():.2f} g   "
      f"quartiles {s.quantile(.25):.2f} to {s.quantile(.75):.2f}")
print(f"                            min {s.min():.2f} g   max {s.max():.2f} g   "
      f"mean {s.mean():.2f} g")
display(pd.crosstab(raw.line, raw["shift"], margins=True, margins_name="all"))

rows 78   unique bag_id 76   fill_g recorded on 75
fill_g, raw and uncleaned:  median 341.40 g   quartiles 340.00 to 342.75
                            min 0.00 g   max 3410.00 g   mean 373.56 g

fig, ax = plt.subplots(1, 2, figsize=(11, 3.4))
ax[0].hist(raw.fill_g.dropna(), bins=40, color=GOLD, alpha=0.85)
ax[0].set_title("Raw weights: impossible values")
ax[0].set_xlabel("fill_g (grams)"); ax[0].set_ylabel("bags")
plaus = raw.fill_g[(raw.fill_g > 300) & (raw.fill_g < 380)]
ax[1].hist(plaus, bins=20, color=AMBER, alpha=0.9)
ax[1].axvline(TARGET, color=RD, lw=2, ls="--", label="label 340 g")
ax[1].set_title("Zoomed in: the real distribution")
ax[1].set_xlabel("fill_g (grams)"); ax[1].legend()
plt.tight_layout(); plt.show()

before = len(raw)
df = raw.drop_duplicates()
after_dupes = len(df)
df = df.dropna(subset=["fill_g"])
after_na = len(df)
bad = df[(df.fill_g <= 300) | (df.fill_g >= 380)]
df = df[(df.fill_g > 300) & (df.fill_g < 380)].copy()
print(f"start rows          : {before}")
print(f"after drop_duplicates: {after_dupes}   (removed {before-after_dupes} duplicate rows)")
print(f"after dropna fill_g : {after_na}   (removed {after_dupes-after_na} blank weights)")
print("impossible values removed:", sorted(bad.fill_g.tolist()))
print(f"clean analysis sample: n = {len(df)}")

start rows          : 78
after drop_duplicates: 76   (removed 2 duplicate rows)
after dropna fill_g : 73   (removed 3 blank weights)
impossible values removed: [0.0, 34.1, 3410.0]
clean analysis sample: n = 70

x = df.fill_g.values
n = len(x)
mean, sd = x.mean(), x.std(ddof=1)
se = sd/np.sqrt(n)
print(f"n = {n}    mean = {mean:.3f} g    sd = {sd:.3f} g    se = {se:.3f} g")
print(f"observed gap from target: {mean-TARGET:+.3f} g")

n = 70    mean = 341.229 g    sd = 2.105 g    se = 0.252 g
observed gap from target: +1.229 g

sw = stats.shapiro(x)
ks = stats.kstest(x, stats.norm(loc=mean, scale=sd).cdf)
print(f"Shapiro-Wilk : W = {sw.statistic:.3f}, p = {sw.pvalue:.3f}   -> {'normal enough' if sw.pvalue>0.05 else 'not normal'}")
print(f"Kolmogorov-Smirnov: D = {ks.statistic:.3f}, p = {ks.pvalue:.3f}   -> {'normal enough' if ks.pvalue>0.05 else 'not normal'}")
print(f"skew = {stats.skew(x):.2f}   kurtosis = {stats.kurtosis(x):.2f}")

fig, ax = plt.subplots(1, 2, figsize=(11, 3.4))
ax[0].hist(x, bins=18, color=AMBER, alpha=0.9, density=True)
xs = np.linspace(x.min(), x.max(), 200)
ax[0].plot(xs, stats.norm.pdf(xs, mean, sd), color=DARK, lw=2, label="normal fit")
ax[0].axvline(TARGET, color=RD, lw=2, ls="--", label="target 340 g")
ax[0].set_title("Clean weights with a normal curve overlaid"); ax[0].set_xlabel("fill_g (grams)"); ax[0].legend(fontsize=9)
stats.probplot(x, dist="norm", plot=ax[1])
ax[1].get_lines()[0].set(marker="o", markerfacecolor=GOLD, markeredgecolor="white", markersize=6)
ax[1].get_lines()[1].set(color=DARK, lw=2)
ax[1].set_title("Q-Q: the points hug the line")
plt.tight_layout(); plt.show()

Shapiro-Wilk : W = 0.973, p = 0.142   -> normal enough
Kolmogorov-Smirnov: D = 0.080, p = 0.735   -> normal enough
skew = -0.58   kurtosis = 1.01

t, p = stats.ttest_1samp(x, TARGET)
print(f"one-sample t-test:  t({n-1}) = {t:.3f},  p = {p:.2e}")

res = pg.ttest(x, TARGET)
print("\npingouin one-sample t-test:")
print(res[["T","dof","alternative","cohen_d","power"]].round(4).to_string(index=False))
ci = res["CI95"].iloc[0]
d  = float(res["cohen_d"].iloc[0])
print(f"\n95% CI for the mean fill weight: [{ci[0]:.2f}, {ci[1]:.2f}] g")
print(f"Cohen's d = {d:.3f}  ({'small' if d<0.5 else 'medium' if d<0.8 else 'large'} effect)")

w, pw = stats.wilcoxon(x - TARGET)
print(f"\nWilcoxon signed-rank (robustness): W = {w:.1f}, p = {pw:.2e}  -> same conclusion")

one-sample t-test:  t(69) = 4.884,  p = 6.47e-06

pingouin one-sample t-test:
     T  dof alternative  cohen_d  power
4.8841   69   two-sided   0.5838 0.9978

95% CI for the mean fill weight: [340.73, 341.73] g
Cohen's d = 0.584  (medium effect)

Wilcoxon signed-rank (robustness): W = 404.5, p = 7.50e-06  -> same conclusion

# The result in one picture: the sample mean and its 95% CI against the target.
tcrit = stats.t.ppf(0.975, n-1)
lo, hi = mean - tcrit*se, mean + tcrit*se
fig, ax = plt.subplots(figsize=(8.4, 3.4))
ax.hist(x, bins=18, color=AMBER, alpha=0.55)
ax.axvline(TARGET, color=RD, lw=2.4, ls="--", label=f"target {TARGET:.0f} g")
ax.axvline(mean, color=DARK, lw=2.4, label=f"sample mean {mean:.2f} g")
ax.axvspan(lo, hi, color=GOLD, alpha=0.20, label=f"95% CI [{lo:.2f}, {hi:.2f}]")
ax.set_title("The 95% CI sits entirely above 340 g: a real, if small, overfill")
ax.set_xlabel("fill_g (grams)"); ax.set_ylabel("bags"); ax.legend(fontsize=9)
plt.tight_layout(); plt.show()

# A percentile bootstrap: resample the data with replacement, recompute, repeat.
# Chapter 75 introduced this; it is how we get an interval when no formula exists.
def boot_ci(data, stat, B=10000, seed=11, ci=95):
    rng = np.random.default_rng(seed)
    data = np.asarray(data)
    reps = np.array([stat(rng.choice(data, size=len(data), replace=True)) for _ in range(B)])
    half = (100 - ci) / 2
    return np.percentile(reps, half), np.percentile(reps, 100 - half), reps

def boot_ci_pairs(a, b, stat, B=10000, seed=11, ci=95):
    """Bootstrap for two paired/related arrays: resample ROW indices together."""
    rng = np.random.default_rng(seed)
    a, b = np.asarray(a), np.asarray(b)
    idx = np.arange(len(a))
    reps = np.array([stat(a[i], b[i]) for i in
                     (rng.choice(idx, size=len(idx), replace=True) for _ in range(B))])
    reps = reps[np.isfinite(reps)]
    half = (100 - ci) / 2
    return np.percentile(reps, half), np.percentile(reps, 100 - half), reps

def boot_ci_groups(groups, stat, B=10000, seed=11, ci=95):
    """Bootstrap for k independent groups: resample within each group."""
    rng = np.random.default_rng(seed)
    groups = [np.asarray(g) for g in groups]
    reps = []
    for _ in range(B):
        rs = [rng.choice(g, size=len(g), replace=True) for g in groups]
        v = stat(rs)
        if np.isfinite(v):
            reps.append(v)
    reps = np.array(reps)
    half = (100 - ci) / 2
    return np.percentile(reps, half), np.percentile(reps, 100 - half), reps

# 1. The overfill itself. The CI for the mean shifts straight onto the overfill scale.
overfill = mean - TARGET
over_lo, over_hi = lo - TARGET, hi - TARGET
print(f"overfill per bag       {overfill:+.2f} g    95% CI [{over_lo:+.2f}, {over_hi:+.2f}] g")

# 2. The same interval expressed as the thing a plant manager budgets for.
print(f"per million bags       {overfill*1000:,.0f} kg   95% CI [{over_lo*1000:,.0f}, {over_hi*1000:,.0f}] kg")

# 3. An interval for the effect size, which is almost always reported bare.
d_lo, d_hi, _ = boot_ci(x, lambda s: (s.mean() - TARGET) / s.std(ddof=1))
print(f"Cohen's d              {d:.2f}         95% CI [{d_lo:.2f}, {d_hi:.2f}]")

# 4. Margin of error: half the width of the interval for the mean.
print(f"\nmargin of error on the mean fill: +/- {(hi-lo)/2:.2f} g at 95% confidence")
print(f"to halve it we would need roughly {n*4} bags, since precision improves with the square root of n")

overfill per bag       +1.23 g    95% CI [+0.73, +1.73] g
per million bags       1,229 kg   95% CI [727, 1,730] kg

Cohen's d              0.58         95% CI [0.32, 0.94]

margin of error on the mean fill: +/- 0.50 g at 95% confidence
to halve it we would need roughly 280 bags, since precision improves with the square root of n

overfill = mean - TARGET
print(f"average overfill: {overfill:.2f} g per bag")
print(f"over a 1,000,000-bag run that is about {overfill*1_000_000/1000:,.0f} kg of extra coffee given away")
print("\nStatistically significant AND practically meaningful at scale,")
print("but every bag is still ABOVE the 340 g label, so no customer is shorted.")

average overfill: 1.23 g per bag
over a 1,000,000-bag run that is about 1,229 kg of extra coffee given away

Statistically significant AND practically meaningful at scale,
but every bag is still ABOVE the 340 g label, so no customer is shorted.

	bag_id	line	shift	fill_g
0	QC-001	Line 1	Night	342.2
1	QC-002	Line 2	Morning	343.4
2	QC-003	Line 2	Morning	342.2
3	QC-004	Line 2	Afternoon	337.9
4	QC-005	Line 1	Morning	343.7

Is the Coffee Underfilled?

Step 1 · Goal and hypotheses¶

Step 2 · Know the data types¶

Step 3 · Describe the sample¶

Step 4 · Visualize (first look, warts and all)¶

Step 5 · Prepare and clean¶

Step 6 · Name the design¶

Step 7 · Check the assumption¶

Step 8 · Decide, and fix if needed¶

Step 9 · Choose the test¶

Step 10 · Run the test and report fully¶

Step 11 · Estimate, do not just test¶

Step 12 · Interpret in plain language¶

Step 13 · Ethics, bias, and limits¶