""" 04-09 Hypothesis Testing Toolkit Chapters 10-01, 10-02, 11-01, 11-02, 11-03 Run from the project folder: python python/analysis.py """ import os import numpy as np import pandas as pd from scipy import stats import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt fill = pd.read_csv("data/fill_volume.csv")["ml"].to_numpy() methods = pd.read_csv("data/two_methods.csv") pairs = pd.read_csv("data/before_after.csv") def report(name, H0, H1, stat_name, stat, df, pval, crit, alpha, effect, ci, conclusion): print("\n" + "=" * 74) print(name) print("=" * 74) print(f"H0: {H0}") print(f"H1: {H1}") print(f"alpha = {alpha:.2f}\n") df_txt = "-" if df is None else f"{df:.3f}" print(f"{stat_name:<10} = {stat:.5f} df = {df_txt}") print(f"critical = {' / '.join(f'{c:.5f}' for c in np.atleast_1d(crit))}") print(f"p-value = {pval:.6f}") print(f"decision = {'REJECT H0' if pval <= alpha else 'FAIL TO REJECT H0'}") print(f"effect = {effect}") print(f"95% CI = ({ci[0]:.4f}, {ci[1]:.4f})") print(f"\n{conclusion}") def interp_d(d): a = abs(d) return ("negligible" if a < 0.2 else "small" if a < 0.5 else "medium" if a < 0.8 else "large") def n_outliers(x): q1, q3 = np.percentile(x, [25, 75]) iqr = q3 - q1 return int(((x < q1 - 1.5 * iqr) | (x > q3 + 1.5 * iqr)).sum()) # ════════════════════════════════════════════════════════════════════════════ # TEST 1 — ONE-SAMPLE t # ════════════════════════════════════════════════════════════════════════════ print("\n>>> ASSUMPTION CHECK — fill volume") print(f"n = {fill.size}, Shapiro-Wilk p = {stats.shapiro(fill).pvalue:.4f}, " f"outliers = {n_outliers(fill)}") n1 = fill.size se1 = fill.std(ddof=1) / np.sqrt(n1) r1 = stats.ttest_1samp(fill, 250) d1 = (fill.mean() - 250) / fill.std(ddof=1) ci1 = r1.confidence_interval(0.95) report("TEST 1 — ONE-SAMPLE t : fill volume vs the 250 ml target", "mu = 250", "mu != 250 (two-tailed: over- AND under-filling matter)", "t", r1.statistic, n1 - 1, r1.pvalue, [-stats.t.ppf(0.975, n1 - 1), stats.t.ppf(0.975, n1 - 1)], 0.05, f"Cohen's d = {d1:.4f} ({interp_d(d1)})", (ci1.low, ci1.high), f"The mean fill volume is {fill.mean():.2f} ml. There " f"{'IS' if r1.pvalue <= 0.05 else 'is NOT'} sufficient evidence at\n" f"the 5% level that the machine differs from its 250 ml target.") print(f"\nnonparametric check (Wilcoxon signed-rank): " f"p = {stats.wilcoxon(fill - 250).pvalue:.5f}") # ════════════════════════════════════════════════════════════════════════════ # TEST 2 — TWO-SAMPLE t # ════════════════════════════════════════════════════════════════════════════ A = methods.loc[methods["method"] == "A", "score"].to_numpy() B = methods.loc[methods["method"] == "B", "score"].to_numpy() print("\n>>> ASSUMPTION CHECK — two methods") print(f"n_A = {A.size} (Shapiro p = {stats.shapiro(A).pvalue:.4f}), " f"n_B = {B.size} (Shapiro p = {stats.shapiro(B).pvalue:.4f})") lev = stats.levene(A, B) print(f"variance ratio = {A.var(ddof=1)/B.var(ddof=1):.3f}, " f"Levene p = {lev.pvalue:.4f} -> " + ("variances are compatible; either test is defensible" if lev.pvalue > 0.05 else "variances DIFFER; use Welch")) r2 = stats.ttest_ind(A, B, equal_var=False) # Welch — the safe default ci2 = r2.confidence_interval(0.95) sp2 = ((A.size - 1) * A.var(ddof=1) + (B.size - 1) * B.var(ddof=1)) / (A.size + B.size - 2) d2 = (A.mean() - B.mean()) / np.sqrt(sp2) report("TEST 2 — TWO-SAMPLE t (Welch) : method A vs method B", "mu_A = mu_B", "mu_A != mu_B", "t", r2.statistic, r2.df, r2.pvalue, [-stats.t.ppf(0.975, r2.df), stats.t.ppf(0.975, r2.df)], 0.05, f"Cohen's d = {d2:.4f} ({interp_d(d2)})", (ci2.low, ci2.high), f"Method A averaged {A.mean():.2f} and method B {B.mean():.2f}, a difference " f"of {A.mean()-B.mean():.2f}.\nThere " f"{'IS' if r2.pvalue <= 0.05 else 'is NOT'} sufficient evidence at the 5% " f"level of a real difference.") print(f"\npooled t-test for comparison: " f"p = {stats.ttest_ind(A, B, equal_var=True).pvalue:.5f}") print(f"Mann-Whitney U (nonparametric): p = {stats.mannwhitneyu(A, B).pvalue:.5f}") # ════════════════════════════════════════════════════════════════════════════ # TEST 3 — PAIRED t # ════════════════════════════════════════════════════════════════════════════ d = pairs["before_ms"].to_numpy() - pairs["after_ms"].to_numpy() # + = improvement print("\n>>> ASSUMPTION CHECK — the DIFFERENCES (not the originals)") print(f"n pairs = {d.size}, Shapiro-Wilk p = {stats.shapiro(d).pvalue:.4f}, " f"outliers = {n_outliers(d)}") print(f"correlation between before and after: " f"{np.corrcoef(pairs['before_ms'], pairs['after_ms'])[0,1]:.4f} " "-> pairing is worth it") r3 = stats.ttest_rel(pairs["before_ms"], pairs["after_ms"], alternative="greater") r3two = stats.ttest_rel(pairs["before_ms"], pairs["after_ms"]) ci3 = r3two.confidence_interval(0.95) d3 = d.mean() / d.std(ddof=1) report("TEST 3 — PAIRED t : reaction time before vs after training", "mu_d <= 0", "mu_d > 0 (one-tailed: training is claimed to REDUCE the time)", "t", r3.statistic, r3.df, r3.pvalue, stats.t.ppf(0.95, r3.df), 0.05, f"Cohen's d = {d3:.4f} ({interp_d(d3)})", (ci3.low, ci3.high), f"Mean improvement {d.mean():.2f} ms. There " f"{'IS' if r3.pvalue <= 0.05 else 'is NOT'} sufficient evidence at the\n" f"5% level that the training reduces reaction time.") ind = stats.ttest_ind(pairs["before_ms"], pairs["after_ms"], equal_var=False, alternative="greater") print("\nTHE COST OF IGNORING THE PAIRING") print(f" paired t = {r3.statistic:7.4f}, p = {r3.pvalue:.6f}") print(f" independent t = {ind.statistic:7.4f}, p = {ind.pvalue:.6f} <- much weaker") print(f"Wilcoxon signed-rank (nonparametric): " f"p = {stats.wilcoxon(d, alternative='greater').pvalue:.6f}") # ════════════════════════════════════════════════════════════════════════════ # TEST 4 — Z-TEST FOR A PROPORTION # ════════════════════════════════════════════════════════════════════════════ poll_path = os.path.join("..", "04-08-confidence-interval-estimator", "data", "poll_responses.csv") if os.path.exists(poll_path): poll = pd.read_csv(poll_path)["supports"] x = int((poll == "Yes").sum()) n4 = int(poll.size) p0 = 0.50 print("\n>>> ASSUMPTION CHECK — proportion") print(f"n*p0 = {n4*p0:.0f} and n*(1-p0) = {n4*(1-p0):.0f} (both >= 5)") se4 = np.sqrt(p0 * (1 - p0) / n4) # p0, NOT p-hat ph = x / n4 z4 = (ph - p0) / se4 p4 = 2 * stats.norm.sf(abs(z4)) se_ci = np.sqrt(ph * (1 - ph) / n4) # p-hat for the INTERVAL ci4 = (ph - 1.959964 * se_ci, ph + 1.959964 * se_ci) report("TEST 4 — Z-TEST FOR A PROPORTION : is support different from 50%?", "p = 0.50", "p != 0.50", "z", z4, None, p4, [-1.959964, 1.959964], 0.05, f"p-hat - p0 = {ph-p0:+.4f} ({100*(ph-p0):.1f} percentage points)", ci4, f"{x} of {n4} respondents ({100*ph:.1f}%) said Yes. There " f"{'IS' if p4 <= 0.05 else 'is NOT'} sufficient\n" f"evidence at the 5% level that support differs from 50%.") print("\nNOTE the two different standard errors:") print(f" TEST uses p0: SE = {se4:.6f}") print(f" INTERVAL uses p-hat: SE = {se_ci:.6f}") else: print("\n(Test 4 skipped: poll_responses.csv lives in project 04-08.)") # ── PLOTS ─────────────────────────────────────────────────────────────────── fig, axes = plt.subplots(2, 2, figsize=(11, 8.5)) axes[0, 0].boxplot(fill, vert=False, showmeans=True) axes[0, 0].axvline(250, color="#b4122e", lw=2) axes[0, 0].axvline(fill.mean(), color="#0FA3A3", lw=2, ls="--") axes[0, 0].set(title="Test 1: fill volume vs 250 ml target", xlabel="ml", yticks=[]) axes[0, 1].boxplot([A, B], tick_labels=["A", "B"]) axes[0, 1].set(title="Test 2: method A vs B", ylabel="score") for b, a in zip(pairs["before_ms"], pairs["after_ms"]): axes[1, 0].plot([0, 1], [b, a], "o-", color="#5B2A86", alpha=0.7) axes[1, 0].set(xticks=[0, 1], xticklabels=["Before", "After"], ylabel="reaction time (ms)", title="Test 3: each line is one subject") axes[1, 1].hist(d, bins=8, color="#0B7A7A", edgecolor="white") axes[1, 1].axvline(0, color="#b4122e", lw=2) axes[1, 1].axvline(d.mean(), color="#0FA3A3", lw=2, ls="--") axes[1, 1].set(title="Test 3: differences (before - after)", xlabel="ms improved") plt.tight_layout() plt.savefig("testing_plots.png", dpi=110) print("\nWrote testing_plots.png")