S022 · Statistics — CLT, Hypothesis Testing, Confidence Intervals

"""stats_playground.py — CLT, t-test, bootstrap, p-hacking demo.
 
Run with:  uv run stats_playground.py
Deps:     numpy, scipy, matplotlib (only for the CLT figure — comment out if headless)
"""
from __future__ import annotations
import numpy as np
from scipy import stats
 
rng = np.random.default_rng(42)
 
 
# ------------------------------------------------------------------
# 1) Central Limit Theorem — averages of an ugly distribution
# ------------------------------------------------------------------
def demo_clt(n_per_sample: int = 30, n_samples: int = 10_000) -> None:
    """Sample from an exponential (very skewed) and take the mean 10,000 times."""
    means = rng.exponential(scale=2.0, size=(n_samples, n_per_sample)).mean(axis=1)
    print(f"[CLT] population mean=2.0  sample-of-means mean={means.mean():.3f}")
    print(f"[CLT] theoretical SE=σ/√n={2.0/np.sqrt(n_per_sample):.3f}  observed={means.std():.3f}")
    # If you have matplotlib, uncomment to see the bell curve:
    # import matplotlib.pyplot as plt
    # plt.hist(means, bins=60); plt.title("Distribution of sample means"); plt.show()
 
 
# ------------------------------------------------------------------
# 2) Two-sample t-test — is treatment better than control?
# ------------------------------------------------------------------
def demo_ttest() -> None:
    control = rng.normal(loc=100, scale=15, size=200)      # e.g. session length in seconds
    treatment = rng.normal(loc=104, scale=15, size=200)    # small true lift of 4s
    t, p = stats.ttest_ind(treatment, control, equal_var=False)   # Welch's t-test
    diff = treatment.mean() - control.mean()
    print(f"[t-test] observed diff={diff:+.2f}s  t={t:.2f}  p={p:.4f}  -> {'REJECT H0' if p<0.05 else 'fail to reject'}")
 
 
# ------------------------------------------------------------------
# 3) Bootstrap 95% CI for the difference in means
# ------------------------------------------------------------------
def bootstrap_ci(a: np.ndarray, b: np.ndarray, n_boot: int = 10_000, alpha: float = 0.05) -> tuple[float, float, float]:
    diffs = np.empty(n_boot)
    n_a, n_b = len(a), len(b)
    for i in range(n_boot):
        ra = rng.choice(a, size=n_a, replace=True)
        rb = rng.choice(b, size=n_b, replace=True)
        diffs[i] = rb.mean() - ra.mean()
    lo, hi = np.quantile(diffs, [alpha/2, 1 - alpha/2])
    return diffs.mean(), lo, hi
 
 
def demo_bootstrap() -> None:
    control = rng.normal(loc=100, scale=15, size=200)
    treatment = rng.normal(loc=104, scale=15, size=200)
    est, lo, hi = bootstrap_ci(control, treatment)
    print(f"[bootstrap] est diff={est:+.2f}s  95% CI=[{lo:+.2f}, {hi:+.2f}]  contains 0? {lo <= 0 <= hi}")
 
 
# ------------------------------------------------------------------
# 4) p-hacking demo — test 20 metrics with no real effect
# ------------------------------------------------------------------
def demo_phacking(n_metrics: int = 20, n_per_group: int = 100) -> None:
    """Run 20 A/B tests where the truth is that NOTHING changed. Count false positives."""
    false_positives = 0
    for i in range(n_metrics):
        a = rng.normal(size=n_per_group)
        b = rng.normal(size=n_per_group)     # same distribution — truth: no effect
        _, p = stats.ttest_ind(a, b)
        if p < 0.05:
            false_positives += 1
            print(f"  metric #{i:02d}: p={p:.3f}  <-- 'SIGNIFICANT' (but it isn't)")
    print(f"[p-hacking] {false_positives}/{n_metrics} tests came out ‘significant’ despite zero true effect.")
    print("            Expected number of false positives at α=0.05: n·α =", n_metrics * 0.05)
 
 
# ------------------------------------------------------------------
# 5) Peeking — early-stopping inflates false positives
# ------------------------------------------------------------------
def demo_peeking(n_final: int = 2000, n_peeks: int = 20, n_experiments: int = 1000) -> None:
    hits = 0
    for _ in range(n_experiments):
        a = rng.normal(size=n_final)
        b = rng.normal(size=n_final)          # again, truth: no effect
        step = n_final // n_peeks
        for k in range(step, n_final + 1, step):
            _, p = stats.ttest_ind(a[:k], b[:k])
            if p < 0.05:
                hits += 1
                break                          # ‘we won! ship it!’
    rate = hits / n_experiments
    print(f"[peeking] false-positive rate with {n_peeks} peeks: {rate:.1%}  (should be ~5%!)")
 
 
if __name__ == "__main__":
    demo_clt()
    demo_ttest()
    demo_bootstrap()
    demo_phacking()
    demo_peeking()