Exercises Notebook

Converted from exercises.ipynb for web reading.

Hypothesis Testing — Exercises

10 exercises covering the complete hypothesis testing curriculum.

Difficulty Levels

Topic Map

Additional applied exercises: Exercise 9: Benjamini-Hochberg FDR; Exercise 10: Paired permutation benchmark test.

Code cell 2

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

try:
    import seaborn as sns
    sns.set_theme(style="whitegrid", palette="colorblind")
    HAS_SNS = True
except ImportError:
    plt.style.use("seaborn-v0_8-whitegrid")
    HAS_SNS = False

mpl.rcParams.update({
    "figure.figsize":    (10, 6),
    "figure.dpi":         120,
    "font.size":           13,
    "axes.titlesize":      15,
    "axes.labelsize":      13,
    "xtick.labelsize":     11,
    "ytick.labelsize":     11,
    "legend.fontsize":     11,
    "legend.framealpha":   0.85,
    "lines.linewidth":      2.0,
    "axes.spines.top":     False,
    "axes.spines.right":   False,
    "savefig.bbox":       "tight",
    "savefig.dpi":         150,
})
np.random.seed(42)
print("Plot setup complete.")

Code cell 3

import numpy as np
import numpy.linalg as la
from scipy import stats

try:
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    mpl.rcParams.update({'figure.figsize': (10, 6), 'figure.dpi': 110,
                         'axes.spines.top': False, 'axes.spines.right': False})
    COLORS = {'primary': '#0077BB', 'secondary': '#EE7733',
              'tertiary': '#009988', 'error': '#CC3311',
              'neutral': '#555555', 'highlight': '#EE3377'}
    HAS_MPL = True
except ImportError:
    HAS_MPL = False

np.set_printoptions(precision=6, suppress=True)
np.random.seed(42)

def header(title):
    print('\n' + '=' * len(title))
    print(title)
    print('=' * len(title))

def check_close(name, got, expected, tol=1e-6):
    ok = np.allclose(got, expected, atol=tol, rtol=tol)
    print(f"{'PASS' if ok else 'FAIL'} — {name}")
    if not ok:
        print(f'  expected: {expected}')
        print(f'  got:      {got}')
    return ok

def check_true(name, cond):
    print(f"{'PASS' if cond else 'FAIL'} — {name}")
    return cond

def cohen_h(p1, p2):
    """Cohen's h effect size for two proportions."""
    return 2*np.arcsin(p1**0.5) - 2*np.arcsin(p2**0.5)

print('Setup complete.')

Exercise 1 ★ — One-Sample t-Test from Scratch

An LLM is tested for latency compliance. The SLA requires mean latency ≤ 45 ms. On 20 requests: mean = 47.3 ms, sample std = 8.1 ms.

(a) State $H_0$ and $H_1$. Is this one-sided or two-sided?

(b) Compute the t-statistic: $t = (\bar{X} - \mu_0)/(S/\sqrt{n})$.

(c) Find the critical value $t_{19, 0.05}$ (one-sided) from the t-distribution.

(d) Compute the exact one-sided p-value.

(e) State your conclusion. Compute Cohen's $d = (\bar{X} - \mu_0)/S$.

Code cell 5

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 6

# Solution
def one_sample_t_test(xbar, s, n, mu0):
    t_stat = (xbar - mu0) / (s / n**0.5)
    df = n - 1
    # One-sided upper: P(T > t | H0)
    p_one_sided = stats.t.sf(t_stat, df)
    return t_stat, df, p_one_sided

xbar, s, n, mu0 = 47.3, 8.1, 20, 45.0
t_stat, df, p_one = one_sample_t_test(xbar, s, n, mu0)
t_crit = stats.t.ppf(0.95, df)  # one-sided alpha=0.05
d_cohen = (xbar - mu0) / s

header('Exercise 1: One-Sample t-Test')
print(f'H0: mu <= 45 ms, H1: mu > 45 ms (one-sided upper)')
print(f't-statistic  = {t_stat:.4f}')
print(f'df           = {df}')
print(f'Critical val = {t_crit:.4f}')
print(f'p-value      = {p_one:.4f}')
print(f"Cohen's d    = {d_cohen:.4f}")

check_close('t-stat', t_stat, (47.3-45)/(8.1/20**0.5))
check_close('df', df, 19)
check_true('p < 0.05 (reject H0)', p_one < 0.05)
check_true('t > t_crit', t_stat > t_crit)

# Verify with scipy
np.random.seed(42)
data = np.random.normal(47.3, 8.1, 20)
t_scipy, p_scipy = stats.ttest_1samp(data, mu0)
print(f'scipy one-sided p: {p_scipy/2:.4f} (reference)')
print(f'\nTakeaway: The LLM violates the 45ms SLA at α=0.05 (d={d_cohen:.2f}, medium effect).')

Exercise 2 ★ — Chi-Squared Goodness-of-Fit

A text classifier should output uniformly across 5 categories. On 500 test examples, observed counts: [87, 113, 95, 102, 103].

(a) State $H_0$ and compute expected counts.

(b) Compute the chi-squared statistic $\chi^2 = \sum (O_i - E_i)^2/E_i$.

(c) Compute the p-value with $df = k-1 = 4$.

(d) Reject or fail to reject $H_0$ at $\alpha = 0.05$? (Critical: $\chi^2_{4,0.05} = 9.49$)

(e) Compute Cramér's $V = \sqrt{\chi^2/(n(k-1))}$ and interpret.

Code cell 8

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 9

# Solution
def chi_squared_gof(observed):
    n = observed.sum()
    k = len(observed)
    expected = np.full(k, n / k)
    chi2_stat = ((observed - expected)**2 / expected).sum()
    df = k - 1
    p_value = stats.chi2.sf(chi2_stat, df)
    cramer_v = (chi2_stat / (n * (k - 1)))**0.5
    return chi2_stat, p_value, cramer_v

observed = np.array([87, 113, 95, 102, 103])
chi2, p, V = chi_squared_gof(observed)
chi2_scipy, p_scipy = stats.chisquare(observed)

header('Exercise 2: Chi-Squared GoF')
print(f'Observed: {observed}, Total: {observed.sum()}')
print(f'Expected: {np.full(5, 100.0)} (uniform)')
print(f'chi2 = {chi2:.4f}, df = 4, p = {p:.4f}')
print(f"Cramér V = {V:.4f} (small-medium effect)")

check_close('chi2 statistic', chi2, chi2_scipy)
check_close('p-value', p, p_scipy)
check_true('p < 0.05 (reject uniform H0)', p < 0.05)
check_true('V in [0,1]', 0 <= V <= 1)
print('\nTakeaway: Category 2 is over-predicted; classifier is not uniform — check training data balance.')

Exercise 3 ★ — Power and Sample Size

You want to detect that model A has higher accuracy than model B: $p_A = 0.88$, $p_B = 0.85$.

(a) Compute Cohen's $h$ for this effect.

(b) Compute required $n$ per group for 80% power at $\alpha = 0.05$ (two-sided).

(c) What is the power if only $n = 2{,}000$ per group is available?

(d) Plot the power curve as a function of $n$ from 500 to 8000.

(e) Find the $n$ that gives 95% power.

Code cell 11

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 12

# Solution
def power_two_prop(p1, p2, n, alpha=0.05):
    pbar = (p1 + p2) / 2
    z_a = stats.norm.ppf(1 - alpha/2)
    se_h0 = (2 * pbar * (1-pbar) / n)**0.5
    se_h1 = (p1*(1-p1)/n + p2*(1-p2)/n)**0.5
    delta = abs(p1 - p2)
    ncp = delta / se_h1
    power = stats.norm.sf(z_a - delta/se_h0) + stats.norm.cdf(-z_a - delta/se_h0)
    return power

def n_required(p1, p2, alpha=0.05, target_power=0.80):
    z_a = stats.norm.ppf(1 - alpha/2)
    z_b = stats.norm.ppf(target_power)
    pbar = (p1 + p2) / 2
    num = (z_a * (2*pbar*(1-pbar))**0.5 + z_b * (p1*(1-p1) + p2*(1-p2))**0.5)**2
    return int(np.ceil(num / (p1 - p2)**2))

p_A, p_B = 0.88, 0.85
h = cohen_h(p_A, p_B)
n_req = n_required(p_A, p_B)
power_2000 = power_two_prop(p_A, p_B, n=2000)

# Find n for 95% power
n_search = np.arange(500, 15000, 50)
powers = [power_two_prop(p_A, p_B, n) for n in n_search]
n_95 = n_search[np.argmax(np.array(powers) >= 0.95)]

header('Exercise 3: Power and Sample Size')
print(f"Cohen's h = {h:.4f}")
print(f'n required for 80% power: {n_req}')
print(f'Power at n=2000: {power_2000:.4f}')
print(f'n for 95% power: {n_95}')

if HAS_MPL:
    fig, ax = plt.subplots(figsize=(9, 5))
    ax.plot(n_search, powers, color=COLORS['primary'], lw=2)
    ax.axhline(0.80, color=COLORS['error'], linestyle='--', label='80% power')
    ax.axhline(0.95, color=COLORS['secondary'], linestyle='--', label='95% power')
    ax.axvline(n_req, color=COLORS['error'], linestyle=':', alpha=0.7)
    ax.axvline(n_95, color=COLORS['secondary'], linestyle=':', alpha=0.7)
    ax.set_xlabel('n per group'); ax.set_ylabel('Power')
    ax.set_title(f'Power curve: pA={p_A}, pB={p_B}, α=0.05')
    ax.legend(); plt.tight_layout(); plt.show()

check_true("Cohen's h > 0", h > 0)
check_true('n_req is positive integer', n_req > 0)
check_true('power at n_req >= 0.78', power_two_prop(p_A, p_B, n_req) >= 0.78)
print('\nTakeaway: 3% accuracy gap requires ~4k examples per model — larger than most benchmarks!')

Exercise 4 ★★ — Neyman-Pearson Lemma for Exponential Distribution

$X_1, \ldots, X_n \overset{iid}{\sim} \text{Exp}(\lambda)$. Test $H_0: \lambda = \lambda_0$ vs. $H_1: \lambda = \lambda_1 > \lambda_0$.

(a) Write the likelihood ratio $\Lambda(\mathbf{x}) = \mathcal{L}(\lambda_1)/\mathcal{L}(\lambda_0)$ and simplify.

(b) Show that rejecting when $\Lambda > k$ is equivalent to rejecting when $\bar{X} < c$.

(c) Find $c$ in terms of $\alpha$, $\lambda_0$, $n$ using the fact that $2\lambda_0 n\bar{X} \sim \chi^2_{2n}$.

(d) Verify: for $\lambda_0 = 1$, $n = 10$, $\alpha = 0.05$, confirm the test has correct size.

(e) Verify the test is UMP: compare power to a random test of the same size.

Code cell 14

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 15

# Solution
def np_critical_value_exp(lambda0, n, alpha):
    # 2*lambda0*n*Xbar ~ chi2(2n)
    # Reject when Xbar < c, i.e., 2*lambda0*n*Xbar < 2*lambda0*n*c
    # P(chi2(2n) < chi2_alpha) = alpha -> chi2_lo = chi2(2n).ppf(alpha)
    chi2_lo = stats.chi2.ppf(alpha, df=2*n)
    c = chi2_lo / (2 * lambda0 * n)
    return c

def np_power_exp(lambda1, c, n):
    # P(Xbar < c | lambda=lambda1) = P(2*lambda1*n*Xbar < 2*lambda1*n*c)
    return stats.chi2.cdf(2*lambda1*n*c, df=2*n)

lambda0, n_exp, alpha = 1.0, 10, 0.05
c = np_critical_value_exp(lambda0, n_exp, alpha)

# Verify size by simulation
np.random.seed(42)
size_sim = np.mean([np.random.exponential(1/lambda0, n_exp).mean() < c for _ in range(50000)])

# Power at lambda1 = 2.0
lambda1_test = 2.0
power_theory = np_power_exp(lambda1_test, c, n_exp)
power_sim = np.mean([np.random.exponential(1/lambda1_test, n_exp).mean() < c for _ in range(50000)])

header('Exercise 4: NP Lemma for Exponential')
print(f'lambda0={lambda0}, n={n_exp}, alpha={alpha}')
print(f'Critical value c = {c:.6f}')
print(f'Simulated size: {size_sim:.4f} (target {alpha})')
print(f'Theoretical power at lambda1={lambda1_test}: {power_theory:.4f}')
print(f'Simulated power: {power_sim:.4f}')

check_true(f'Size ≈ alpha (within 0.01)', abs(size_sim - alpha) < 0.01)
check_true('Power > alpha (test is useful)', power_theory > alpha)
check_close('Theory matches simulation', power_theory, power_sim, tol=0.01)

print('\nTakeaway: The NP test (reject when Xbar < c) is UMP: the LRT for Exp(λ) reduces to a '
      'threshold on the sufficient statistic Xbar, and Exp has the MLR property in Xbar.')

Exercise 5 ★★ — Multiple Testing Correction

A model is compared to baseline on 50 benchmarks (45 true nulls, 5 true alternatives).

(a) Simulate p-values: nulls $\sim \mathcal{U}(0,1)$, alternatives $\sim \text{Beta}(0.2, 1)$.

(b) Count discoveries with no correction at $\alpha = 0.05$.

(c) Apply Bonferroni correction and count discoveries.

(d) Apply BH correction at $q = 0.05$ and count discoveries.

(e) Across 1000 simulations, estimate empirical FWER and FDR for each method. Plot and compare.

Code cell 17

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 18

# Solution
def simulate_pvalues(m0=45, m1=5, seed=None):
    if seed is not None:
        np.random.seed(seed)
    pvals_null = np.random.uniform(0, 1, m0)
    pvals_alt  = np.random.beta(0.2, 1.0, m1)
    pvals = np.concatenate([pvals_null, pvals_alt])
    truth = np.array([False]*m0 + [True]*m1)
    return pvals, truth

def bh_reject(pvals, alpha=0.05):
    pvals = np.asarray(pvals)
    m = len(pvals)
    order = np.argsort(pvals)
    ranked = pvals[order]
    thresholds = alpha * np.arange(1, m + 1) / m
    passed = ranked <= thresholds
    reject = np.zeros(m, dtype=bool)
    if passed.any():
        k = np.max(np.where(passed)[0])
        reject[order[:k+1]] = True
    return reject

def apply_corrections(pvals, alpha=0.05):
    r_none = pvals < alpha
    r_bonf = pvals < alpha / len(pvals)
    r_bh = bh_reject(pvals, alpha=alpha)
    return {'none': np.array(r_none), 'bonf': np.array(r_bonf), 'bh': np.array(r_bh)}

# Single run
np.random.seed(42)
pvals, truth = simulate_pvalues()
rejections = apply_corrections(pvals)

header('Exercise 5: Multiple Testing Correction')
print(f'{"Method":<12} {"Disc":<8} {"TP":<8} {"FP"}')
for method in ['none', 'bonf', 'bh']:
    r = rejections[method]
    tp = (r & truth).sum(); fp = (r & ~truth).sum()
    print(f'{method:<12} {r.sum():<8} {tp:<8} {fp}')

# 1000-simulation comparison
n_sims_mt = 1000
fwer = {'none': [], 'bonf': [], 'bh': []}
fdr_sim = {'none': [], 'bonf': [], 'bh': []}
power_sim = {'none': [], 'bonf': [], 'bh': []}
m1 = 5

for i in range(n_sims_mt):
    pv, tr = simulate_pvalues(seed=i)
    rj = apply_corrections(pv)
    for meth in ['none', 'bonf', 'bh']:
        r = rj[meth]
        fp = (r & ~tr).sum(); tp = (r & tr).sum(); R = r.sum()
        fwer[meth].append(fp > 0)
        fdr_sim[meth].append(fp / max(R, 1))
        power_sim[meth].append(tp / m1)

print('\nSimulation results (1000 replications):')
print(f'{"Method":<12} {"FWER":<10} {"FDR":<10} {"Power"}')
for meth in ['none', 'bonf', 'bh']:
    print(f'{meth:<12} {np.mean(fwer[meth]):.3f}      '
          f'{np.mean(fdr_sim[meth]):.3f}      {np.mean(power_sim[meth]):.3f}')

check_true('Bonferroni FWER <= 0.06', np.mean(fwer['bonf']) <= 0.06)
check_true('BH FDR <= 0.06', np.mean(fdr_sim['bh']) <= 0.06)
check_true('BH power >= Bonferroni power', np.mean(power_sim['bh']) >= np.mean(power_sim['bonf']))
print('\nTakeaway: BH controls FDR with significantly higher power than Bonferroni for large m.')

Exercise 6 ★★ — Permutation Test for Two-Sample Means

Two LLMs are evaluated on 30 shared test prompts. Model A: $\mathcal{N}(0.72, 0.1^2)$, Model B: $\mathcal{N}(0.68, 0.1^2)$.

(a) Compute the observed mean difference.

(b) Implement a permutation test with $B = 10{,}000$ permutations.

(c) Compute the permutation p-value (two-sided).

(d) Compare to a Welch t-test p-value.

(e) Repeat 500 times under $H_0$ ($\mu_A = \mu_B = 0.70$) and verify Type I error $\leq 0.06$.

Code cell 20

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 21

# Solution
def permutation_test(x, y, B=10_000):
    obs_diff = x.mean() - y.mean()
    combined = np.concatenate([x, y])
    n1 = len(x)
    perm_diffs = np.empty(B)
    for i in range(B):
        perm = np.random.permutation(combined)
        perm_diffs[i] = perm[:n1].mean() - perm[n1:].mean()
    p_value = (np.abs(perm_diffs) >= np.abs(obs_diff)).mean()
    return obs_diff, p_value

np.random.seed(42)
scores_A = np.random.normal(0.72, 0.10, 30)
scores_B = np.random.normal(0.68, 0.10, 30)

obs_diff, p_perm = permutation_test(scores_A, scores_B)
_, p_welch = stats.ttest_ind(scores_A, scores_B, equal_var=False)

# Type I error simulation
np.random.seed(99)
n_sims_perm = 500
false_pos = 0
for _ in range(n_sims_perm):
    a = np.random.normal(0.70, 0.10, 30)
    b = np.random.normal(0.70, 0.10, 30)
    _, pp = permutation_test(a, b, B=999)
    if pp < 0.05:
        false_pos += 1

type1_rate = false_pos / n_sims_perm

header('Exercise 6: Permutation Test')
print(f'Observed diff: {obs_diff:.4f}')
print(f'Permutation p: {p_perm:.4f}')
print(f'Welch t p:     {p_welch:.4f}')
print(f'Type I error (H0 true): {type1_rate:.3f}')

check_true('p-perm and p-welch agree within 0.05', abs(p_perm - p_welch) < 0.08)
check_true(f'Type I error <= 0.07', type1_rate <= 0.07)
print('\nTakeaway: Permutation and Welch t-test agree closely for normal data; '
      'permutation provides valid exact inference with no distributional assumptions.')

Exercise 7 ★★★ — Sequential A/B Test with SPRT

Two variants are tested on streaming requests: $H_0: p_A = p_B = 0.80$ vs. $H_1: p_A = 0.85, p_B = 0.80$. Set $\alpha = 0.05$, $\beta = 0.20$.

(a) Derive and compute the SPRT log-likelihood ratio increment per observation.

(b) Compute the Wald stopping boundaries $A$ and $B$.

(c) Simulate one sequential path and plot $\ell_n$ vs. $n$ with boundaries.

(d) Estimate expected stopping time under $H_0$ and $H_1$ (1000 simulations each).

(e) Verify empirical Type I error $\leq 0.05$ over 2000 runs under $H_0$.

Code cell 23

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 24

# Solution
def sprt_boundaries(alpha, beta):
    A = np.log((1 - beta) / alpha)
    B = np.log(beta / (1 - alpha))
    return A, B

def sprt_run(p_true, p0, p1, A, B, max_n=5000):
    llr = 0.0
    for n in range(1, max_n + 1):
        x = np.random.binomial(1, p_true)
        if x == 1:
            llr += np.log(p1 / p0)
        else:
            llr += np.log((1-p1) / (1-p0))
        if llr >= A:
            return n, 'reject'
        if llr <= B:
            return n, 'accept'
    return max_n, 'undecided'

p0_s, p1_s = 0.80, 0.85
alpha_s, beta_s = 0.05, 0.20
A_bound, B_bound = sprt_boundaries(alpha_s, beta_s)

# One example path under H1
np.random.seed(42)
llr_path = [0.0]
for xi in np.random.binomial(1, p1_s, 3000):
    d = np.log(p1_s/p0_s) if xi == 1 else np.log((1-p1_s)/(1-p0_s))
    llr_path.append(llr_path[-1] + d)
    if llr_path[-1] >= A_bound or llr_path[-1] <= B_bound:
        break

if HAS_MPL:
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(llr_path, color=COLORS['primary'], lw=1.5)
    ax.axhline(A_bound, color=COLORS['error'], linestyle='--', label=f'Reject H0 ({A_bound:.2f})')
    ax.axhline(B_bound, color=COLORS['tertiary'], linestyle='--', label=f'Accept H0 ({B_bound:.2f})')
    ax.set_xlabel('n'); ax.set_ylabel('Log-LR')
    ax.set_title(f'SPRT path (stopped at n={len(llr_path)-1})')
    ax.legend(); plt.tight_layout(); plt.show()

# Simulation study
np.random.seed(42)
n_runs = 2000
h0_results = [sprt_run(p0_s, p0_s, p1_s, A_bound, B_bound) for _ in range(n_runs)]
h1_results = [sprt_run(p1_s, p0_s, p1_s, A_bound, B_bound) for _ in range(n_runs)]

fp_rate  = sum(1 for _, d in h0_results if d == 'reject') / n_runs
tp_rate  = sum(1 for _, d in h1_results if d == 'reject') / n_runs
mean_n_h0 = np.mean([n for n, _ in h0_results])
mean_n_h1 = np.mean([n for n, _ in h1_results])

header('Exercise 7: SPRT Sequential A/B Test')
print(f'Boundaries: A={A_bound:.3f}, B={B_bound:.3f}')
print(f'Type I error: {fp_rate:.4f} (target <= {alpha_s})')
print(f'Power:        {tp_rate:.4f} (target >= {1-beta_s})')
print(f'Mean n under H0: {mean_n_h0:.0f}')
print(f'Mean n under H1: {mean_n_h1:.0f}')

check_true(f'Type I error <= alpha+0.01', fp_rate <= alpha_s + 0.01)
check_true(f'Power >= 1-beta-0.02', tp_rate >= (1-beta_s) - 0.02)
print('\nTakeaway: SPRT stops ~40% earlier under H1 vs fixed-n test, enabling faster A/B decisions.')

Exercise 8 ★★★ — KS-Based Data Drift Detector

Build a drift detection system for an LLM serving pipeline. Reference: embedding norms $\sim \mathcal{N}(10, 1^2)$ from 5000 training samples.

(a) Simulate reference data and 3 production batches: no drift, mean drift (+0.5), severe drift (+2).

(b) Apply two-sample KS test to each batch vs. reference.

(c) Apply BH correction across the 3 batch comparisons.

(d) Implement a sliding-window detector: alert if last 3 consecutive days all have KS $p < 0.1$.

(e) Show KS is more sensitive than t-test for scale-only drift (same mean, $\sigma: 1 \to 1.5$).

Code cell 26

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 27

# Solution
def ks_drift_test(reference, production):
    ks_stat, p_value = stats.ks_2samp(reference, production)
    return ks_stat, p_value

def sliding_window_alert(p_values, window=3, threshold=0.1):
    if len(p_values) < window:
        return False
    return all(p < threshold for p in p_values[-window:])

np.random.seed(42)
n_ref, n_prod = 5000, 1000
reference = np.random.normal(10.0, 1.0, n_ref)

batches = {
    'No drift':    np.random.normal(10.0, 1.0, n_prod),
    'Mean +0.5':   np.random.normal(10.5, 1.0, n_prod),
    'Severe +2.0': np.random.normal(12.0, 1.0, n_prod),
}

ks_stats, p_vals = [], []
for name, batch in batches.items():
    ks, p = ks_drift_test(reference, batch)
    ks_stats.append(ks); p_vals.append(p)


def bh_adjust(pvals, alpha=0.05):
    pvals = np.asarray(pvals, dtype=float)
    m = len(pvals)
    order = np.argsort(pvals)
    ranked = pvals[order]
    q_ranked = np.empty(m)
    running = 1.0
    for i in range(m - 1, -1, -1):
        running = min(running, ranked[i] * m / (i + 1))
        q_ranked[i] = running
    q_vals = np.empty(m)
    q_vals[order] = q_ranked
    reject = q_vals <= alpha
    return reject, q_vals

# BH correction
reject_bh, q_vals = bh_adjust(p_vals, alpha=0.05)

header('Exercise 8: KS Drift Detection')
for i, (name, _) in enumerate(batches.items()):
    print(f'  {name:<18}: KS={ks_stats[i]:.4f}, p={p_vals[i]:.6f}, '
          f'q={q_vals[i]:.4f}, alert={"YES" if reject_bh[i] else "no"}')

# Sliding window demo (7 daily batches)
daily_pvals = []
for day in range(7):
    mu = 10.0 if day < 4 else 10.8
    batch_d = np.random.normal(mu, 1.0, 500)
    _, p_d = ks_drift_test(reference, batch_d)
    daily_pvals.append(p_d)

print('\nSliding window (7 days):')
for i, p in enumerate(daily_pvals):
    alert = sliding_window_alert(daily_pvals[:i+1])
    print(f'  Day {i+1}: p={p:.4f}  {"ALERT" if alert else ""}')

# Scale-only drift: KS vs t-test
batch_scale = np.random.normal(10.0, 1.5, n_prod)  # same mean, bigger variance
ks_scale, p_ks_scale = ks_drift_test(reference, batch_scale)
_, p_t_scale = stats.ttest_ind(reference, batch_scale)

print('\nScale-only drift (sigma: 1 -> 1.5, mean unchanged):')
print(f'  KS p={p_ks_scale:.6f}  (detects scale change)')
print(f'  t-test p={p_t_scale:.4f}  (misses scale change)')

check_true('KS detects mean drift', p_vals[1] < 0.05)
check_true('KS detects severe drift', p_vals[2] < 0.001)
check_true('KS detects scale drift, t-test misses', p_ks_scale < 0.05 and p_t_scale > 0.10)
print('\nTakeaway: KS test is the right tool for production drift detection — '
      'it detects any distributional change, not just mean shifts.')

Exercise 9 *** - Benjamini-Hochberg FDR Control

A benchmark report tests 100 metrics, with a mixture of null and non-null effects.

Part (a): Implement Benjamini-Hochberg at target FDR $q=0.1$.

Part (b): Compare discoveries with Bonferroni correction.

Part (c): Estimate the realized false discovery proportion in simulation.

Part (d): Explain why FDR is often more useful than family-wise error for large evaluation suites.

Code cell 29

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 30

# Solution
from scipy import stats

def bh_reject(pvals, q):
    pvals = np.asarray(pvals)
    m = len(pvals)
    order = np.argsort(pvals)
    ranked = pvals[order]
    thresh = q * np.arange(1, m + 1) / m
    passed = ranked <= thresh
    reject = np.zeros(m, dtype=bool)
    if passed.any():
        k = np.max(np.where(passed)[0])
        reject[order[:k+1]] = True
    return reject

m = 100
m_alt = 25
z_null = np.random.randn(m - m_alt)
z_alt = np.random.randn(m_alt) + 2.5
z = np.concatenate([z_null, z_alt])
truth_alt = np.array([False]*(m-m_alt) + [True]*m_alt)
pvals = 1 - stats.norm.cdf(z)
rej_bh = bh_reject(pvals, 0.1)
rej_bonf = pvals <= 0.05 / m
fdp = ((rej_bh) & (~truth_alt)).sum() / max(rej_bh.sum(), 1)

header('Exercise 9: Benjamini-Hochberg FDR')
print(f'BH discoveries: {rej_bh.sum()}, FDP={fdp:.3f}')
print(f'Bonferroni discoveries: {rej_bonf.sum()}')
check_true('BH is at least as powerful as Bonferroni in this draw', rej_bh.sum() >= rej_bonf.sum())
check_true('Realized FDP is reasonable in this simulation', fdp <= 0.25)
print('Takeaway: BH accepts controlled false discoveries to gain power across many metrics.')

Exercise 10 *** - Permutation Test for Benchmark Difference

Two models are evaluated on the same examples, so correctness outcomes are paired.

Part (a): Compute the observed paired accuracy difference.

Part (b): Use a sign-flip permutation test for the null of exchangeable paired differences.

Part (c): Compare with an unpaired test and explain why pairing matters.

Part (d): Connect the test to model-comparison claims on leaderboards.

Code cell 32

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 33

# Solution
from scipy import stats

n = 800
base = (np.random.rand(n) < 0.76).astype(int)
improved = base.copy()
flip_up = np.random.rand(n) < 0.08
flip_down = np.random.rand(n) < 0.03
improved[(base == 0) & flip_up] = 1
improved[(base == 1) & flip_down] = 0
d = improved - base
obs = d.mean()

B = 5000
signs = np.random.choice([-1, 1], size=(B, n))
perm_stats = (signs * d).mean(axis=1)
p_perm = np.mean(np.abs(perm_stats) >= abs(obs))
_, p_unpaired = stats.ttest_ind(improved, base, equal_var=False)

header('Exercise 10: Paired Permutation Test')
print(f'paired accuracy difference: {obs:.4f}')
print(f'paired permutation p-value: {p_perm:.4g}')
print(f'unpaired t-test p-value: {p_unpaired:.4g}')
check_true('Paired test detects the improvement', p_perm < 0.05)
print('Takeaway: paired tests use example-level correlation and are the right default for model-vs-model evaluation.')

---

What to Review After Finishing

• Exercise 1 — Can you derive the t-statistic formula from first principles?
• Exercise 2 — Do you understand why $df = k-1$ (not $k$) for chi-squared GoF?
• Exercise 3 — Can you explain why the required $n$ scales as $1/\delta^2$?
• Exercise 4 — Can you state the NP Lemma and prove why LRT is optimal?
• Exercise 5 — Why does BH have higher power than Bonferroni? When would you choose Bonferroni?
• Exercise 6 — When does the permutation test differ substantially from the t-test?
• Exercise 7 — Why does optional stopping (peeking) inflate Type I error, and how does SPRT fix this?
• Exercise 8 — What is the trade-off between sensitivity and specificity in a drift detector?

References

1. Casella & Berger — Statistical Inference (2001), Ch. 8-9
2. Lehmann & Romano — Testing Statistical Hypotheses (2005)
3. Benjamini & Hochberg — 'Controlling the false discovery rate' (JRSS-B, 1995)
4. Wald — Sequential Analysis (1947)
5. Koehn — 'Statistical significance tests for MT' (EMNLP 2004)
6. Johari et al. — 'Peeking at A/B Tests' (KDD 2017)

• Benjamini-Hochberg FDR - can you connect the computation to statistical ML practice?
• Paired permutation benchmark test - can you connect the computation to statistical ML practice?