code/notebooks/_m1_src/0209-variance-components-equivalence.py

align / code /notebooks /_m1_src /0209-variance-components-equivalence.py

Add files using upload-large-folder tool

63f41a6 verified 17 days ago

17.7 kB

# %% [markdown]
# # 0209 - Variance components and equivalence tests for evaluation
# Companion to `notes/0209 - Variance Components and Equivalence Tests.pdf`.  CPU only, numpy/scipy, 1 thread.
# Running toy: a benchmark of `N = 244` problems, `S` independently trained models ("seeds") of one recipe, `m = 16` samples per (seed, problem).
# Per-(seed, problem) solve probability `p_sx = expit(mu_x + a_s + b_sx)`, `mu_x ~ N(0.3, 1.5^2)` (problem difficulty), `a_s ~ N(0, 0.25^2)` (seed), `b_sx ~ N(0, 0.5^2)` (seed x problem).
# %%
import os
os.environ["OMP_NUM_THREADS"] = "1"; os.environ["OPENBLAS_NUM_THREADS"] = "1"; os.environ["MKL_NUM_THREADS"] = "1"
import math, itertools
from fractions import Fraction as Fr
from math import comb
import numpy as np
from scipy import stats
from scipy.special import expit
CHECKS = []
def check(name, err, tol, note=""):
    CHECKS.append((name, float(err), float(tol)))
    assert err <= tol, (name, float(err), tol)
    print(f"[ok] {name:100s} err={float(err):.3e} (tol {tol:.1e}) {note}")
N, m = 244, 16
M0, SX, SA, SB = 0.3, 1.5, 0.25, 0.5
# %% [markdown]
# ## 1. Variance components: problems (A), seeds (B), seed x problem (C), generation (e)
# Exchangeable array `p_sx = mu + A_x + B_s + C_sx` (Hoeffding/ANOVA decomposition; components orthogonal).  With outcomes `c_sx ~ Bin(m, p_sx)` and the score of one run `Y_s = (1/(N m)) sum_x c_sx`:
# `Var(Y_s) = sA2/N + sB2 + sC2/N + se2/(N m)`,   `se2 = E p(1-p)`.   Only `sB2` survives `N -> inf` and `m -> inf`.
# Population components by Gauss-Hermite integration of the logit-normal model (no simulation).
# %%
gh_x, gh_w = np.polynomial.hermite_e.hermegauss(60); gh_w = gh_w / gh_w.sum()     # standard normal nodes/weights
MU = M0 + SX * gh_x; A_ = SA * gh_x; B_ = SB * gh_x
P3 = expit(MU[:, None, None] + A_[None, :, None] + B_[None, None, :])                # [x, s, b]
W3 = gh_w[:, None, None] * gh_w[None, :, None] * gh_w[None, None, :]
Ep = (P3 * W3).sum(); Ep2 = (P3 ** 2 * W3).sum()
Ax = (P3 * gh_w[None, :, None] * gh_w[None, None, :]).sum((1, 2))                    # E_{s,b} p | x  (weights over s,b sum to 1)
Bs = (P3 * gh_w[:, None, None] * gh_w[None, None, :]).sum((0, 2))                    # E_{x,b} p | s
sA2 = (gh_w * (Ax - Ep) ** 2).sum(); sB2 = (gh_w * (Bs - Ep) ** 2).sum()
sT2 = Ep2 - Ep ** 2; sC2 = sT2 - sA2 - sB2; se2 = Ep - Ep2
print(f"  E p = {Ep:.4f}; sigma_A^2 = {sA2:.5f}, sigma_B^2 = {sB2:.6f} (sd {math.sqrt(sB2):.4f}), sigma_C^2 = {sC2:.5f}, sigma_e^2 = E p(1-p) = {se2:.5f}")
def var_Y(S, N=N, m=m): return sA2 / N + sB2 / S + sC2 / (N * S) + se2 / (N * S * m)      # variance of a grand mean over S seeds
print("  per-term contribution to Var(Y_s), S=1:  problems %.3e  seed %.3e  seed*problem %.3e  generation %.3e" % (sA2 / N, sB2, sC2 / N, se2 / (N * m)))
check("1. components are non-negative and sum to Var(p)", abs(sA2 + sB2 + sC2 - sT2) + max(0, -min(sA2, sB2, sC2)), 1e-15)
# simulate replications: each replication redraws problems, seeds, samples
rng = np.random.default_rng(0)
def simulate(R, S, N=N, m=m):
    mu = M0 + SX * rng.standard_normal((R, 1, N)); a = SA * rng.standard_normal((R, S, 1)); b = SB * rng.standard_normal((R, S, N))
    p = expit(mu + a + b)
    c = rng.binomial(m, p)
    return p, c
R1 = 6000
p, c = simulate(R1, 1)
Y = c[:, 0, :].sum(1) / (N * m)
se_rel = math.sqrt(2 / (R1 - 1))
check("1. Var(Y_s) over 6000 replications vs sA2/N + sB2 + sC2/N + se2/(Nm)", abs(Y.var(ddof=1) / var_Y(1) - 1), 4 * se_rel, f"emp {Y.var(ddof=1):.4e} vs {var_Y(1):.4e}")
# MoM / ANOVA estimators with S=4 seeds
S = 4; R2 = 4000
p, c = simulate(R2, S)
ybar = c / m                                                                         # cell means [R,S,N]
gm = ybar.mean((1, 2), keepdims=True)
MSx = S * ((ybar.mean(1) - gm[:, :, 0]) ** 2).sum(1) / (N - 1)
MSs = N * ((ybar.mean(2) - gm[:, :, 0]) ** 2).sum(1) / (S - 1)
res = ybar - ybar.mean(1, keepdims=True) - ybar.mean(2, keepdims=True) + gm
MSr = (res ** 2).sum((1, 2)) / ((S - 1) * (N - 1))
pq = (c * (m - c) / (m * (m - 1))).mean((1, 2))                                      # unbiased for E p(1-p)
est = {"sigma_e^2": (pq, se2), "sigma_C^2": (MSr - pq / m, sC2), "sigma_A^2": ((MSx - MSr) / S, sA2), "sigma_B^2": ((MSs - MSr) / N, sB2)}
for k, (v, t) in est.items():
    z = (v.mean() - t) / (v.std(ddof=1) / math.sqrt(R2))
    check(f"1. ANOVA estimator of {k} is unbiased (4000 reps, S=4)", abs(z), 4.0, f"mean {v.mean():.6f} vs true {t:.6f}  (z={z:+.2f})")
# %% [markdown]
# ## 2. Which term do you pay down with what?  Two recipes, S seeds each, same problems
# `delta_hat = Ybar_A - Ybar_B`.  Under H0 (same recipe) the shared problem effect A_x cancels: `Var(delta_hat) = (2/S) [ sB2 + sC2/N + se2/(N m) ]`.
# The seed term does not shrink with N or m: a floor on the minimum detectable effect.
# %%
def var_diff(S, N=N, m=m): return (2 / S) * (sB2 + sC2 / N + se2 / (N * m))
R3 = 6000
for S in (1, 3):
    mu = M0 + SX * rng.standard_normal((R3, 1, N))
    d = []
    for _ in range(2):
        a = SA * rng.standard_normal((R3, S, 1)); b = SB * rng.standard_normal((R3, S, N))
        cc = rng.binomial(m, expit(mu + a + b)); d.append(cc.sum((1, 2)) / (S * N * m))
    dd = d[0] - d[1]
    check(f"2. Var(delta_hat) for S={S} seeds per recipe vs (2/S)[sB2 + sC2/N + se2/(Nm)]", abs(dd.var(ddof=1) / var_diff(S) - 1), 4 * math.sqrt(2 / (R3 - 1)), f"emp {dd.var(ddof=1):.3e} vs {var_diff(S):.3e}")
z = stats.norm.ppf(0.975) + stats.norm.ppf(0.8)
for S in (1, 3, 5):
    print(f"  S={S}: MDE(80% power, 5% two-sided) at N=244: {z*math.sqrt(var_diff(S))*100:.2f} pp;  N->inf, m->inf floor: {z*math.sqrt(2*sB2/S)*100:.2f} pp")
check("2. floor is strictly positive: MDE(N=1e9, m=1e9, S=3) equals z sqrt(2 sB2 / 3)", abs(z * math.sqrt(var_diff(3, N=10**9, m=10**9)) - z * math.sqrt(2 * sB2 / 3)), 1e-6)
# %% [markdown]
# ## 3. Unbiased pass@k from n samples, and its exact variance
# `e(n,c,k) = 1 - C(n-c,k)/C(n,k)`; `c ~ Bin(n,p)`  =>  `E e = 1 - (1-p)^k`.  Exact variance by summing over `c = 0..n`.
# NAIVE plug-in `1 - (1 - c/n)^k` is biased.  Special cases: `k = 1` gives `p(1-p)/n`;  `k = n` gives `q(1-q)`, `q = 1-(1-p)^n`.
# %%
def e_nck(n, c, k): return 1 - Fr(comb(n - c, k), comb(n, k))
def pk_moments(n, k, p):
    p = Fr(p); pm = [Fr(comb(n, c)) * p ** c * (1 - p) ** (n - c) for c in range(n + 1)]
    E1 = sum(w * e_nck(n, c, k) for c, w in enumerate(pm)); E2 = sum(w * e_nck(n, c, k) ** 2 for c, w in enumerate(pm))
    return E1, E2 - E1 ** 2, sum(w * (1 - (1 - Fr(c, n)) ** k) for c, w in enumerate(pm))
worst = 0; worst_k1 = 0; worst_kn = 0
for n, k, p in [(64, 8, Fr(1, 10)), (64, 32, Fr(1, 50)), (16, 4, Fr(3, 10)), (64, 1, Fr(1, 7)), (64, 64, Fr(1, 100)), (10, 10, Fr(1, 5))]:
    E1, V, Enaive = pk_moments(n, k, p)
    worst = max(worst, abs(float(E1 - (1 - (1 - p) ** k))))
    if k == 1: worst_k1 = max(worst_k1, abs(float(V - p * (1 - p) / n)))
    if k == n: q = 1 - (1 - p) ** n; worst_kn = max(worst_kn, abs(float(V - q * (1 - q))))
check("3. E[e(n,c,k)] = 1-(1-p)^k exactly (6 (n,k,p) triples, rational arithmetic)", worst, 1e-15)
check("3. k=1: Var = p(1-p)/n exactly", worst_k1, 1e-15)
check("3. k=n: Var = q(1-q) with q = 1-(1-p)^n exactly", worst_kn, 1e-15)
print("  n=64:        p     k   true pass@k   E[unbiased]   E[naive plug-in] (bias)   sd(unbiased)   sd/true")
for p in (0.02, 0.1, 0.5):
    for k in (1, 8, 32):
        E1, V, En = pk_moments(64, k, Fr(p).limit_denominator(1000)); g = 1 - (1 - p) ** k
        print(f"           {p:5.2f} {k:4d}   {g:.5f}      {float(E1):.5f}      {float(En):.5f} ({float(En)-g:+.5f})     {math.sqrt(float(V)):.5f}      {math.sqrt(float(V))/g:.3f}")
E1, V, En = pk_moments(64, 32, Fr(1, 50))
check("3. naive plug-in is biased low at (n=64,k=32,p=.02) by more than 0.07", 0.0 if (1 - (1 - 0.02) ** 32) - float(En) > 0.07 else 1.0, 0.0, f"bias {float(En)-(1-0.98**32):+.4f}")
# MC check of the exact variance (independent implementation via binomial draws) and of the problem-level decomposition
rng = np.random.default_rng(1)
n, k, p0 = 64, 8, 0.1
cs = rng.binomial(n, p0, size=2_000_000)
lg = np.array([float(e_nck(n, c_, k)) for c_ in range(n + 1)])
est = lg[cs]
Vex = float(pk_moments(n, k, Fr(1, 10))[1])
check("3. exact Var vs 2e6 Monte Carlo draws (n=64,k=8,p=.1)", abs(est.var() / Vex - 1), 4 * math.sqrt(2 / 2e6))
# across problems: p_x ~ Beta(1.4, 0.7); Var(mean_x e_x) = (1/N)[Var_x(g(p_x)) + E_x Var(e | p_x)]
a_, b_ = 1.4, 0.7; k = 8; n = 64
xs = rng.beta(a_, b_, size=400_000)                              # for the population quantities
cs = rng.binomial(n, xs)
lgk = np.array([float(e_nck(n, c_, k)) for c_ in range(n + 1)])
es = lgk[cs]; gx = 1 - (1 - xs) ** k
pred = (gx.var() + (es - gx).var()) / N
reps = np.array([lgk[rng.binomial(n, rng.beta(a_, b_, size=N))].mean() for _ in range(6000)])
check("3. Var of the benchmark mean pass@8 (N=244, n=64) vs (1/N)[Var_x g + E Var(e|p)]", abs(reps.var(ddof=1) / pred - 1), 5 * math.sqrt(2 / 5999) + 0.01, f"emp {reps.var(ddof=1):.3e} vs {pred:.3e}")
# %% [markdown]
# ## 4. Paired sign-flip test on per-problem differences
# `d_x = e^A_x - e^B_x`, statistic `T = sum_x d_x`, null: each `d_x` symmetric about 0 independently  =>  signs exchangeable: p = fraction of the `2^N` sign vectors with `|T'| >= |T|`.
# Binary `d_x in {-1,0,1}`: reduces to exact McNemar (note 0208).  Exact enumeration for N<=18; Monte Carlo for N=244; size under the null.
# %%
def signflip_exact(d):
    d = np.asarray(d, float); n = len(d)
    signs = np.array(list(itertools.product([-1.0, 1.0], repeat=n)))
    T = np.abs(signs @ d); return float((T >= abs(d.sum()) - 1e-12).mean())
def mcnemar(b, c):
    nd = b + c
    return 1.0 if nd == 0 else min(1.0, 2 * sum(comb(nd, i) for i in range(min(b, c) + 1)) / 2 ** nd)
worst = 0
for b_ in range(0, 8):
    for c_ in range(0, 8):
        if b_ + c_ == 0: continue
        d = np.array([1] * b_ + [-1] * c_ + [0] * 3)
        worst = max(worst, abs(signflip_exact(d) - mcnemar(b_, c_)))
check("4. sign-flip p on {-1,0,+1} equals exact McNemar for all b,c <= 7 (3 concordant zeros)", worst, 1e-12)
check("4. (5,0): two-sided sign-flip p = 0.0625", abs(signflip_exact([1] * 5) - 0.0625), 1e-15)
rng = np.random.default_rng(2)
def signflip_mc(d, B=20000):
    d = np.asarray(d, float); s = rng.choice([-1.0, 1.0], size=(B, len(d))); T = np.abs(s @ d)
    return (1 + (T >= abs(d.sum()) - 1e-12).sum()) / (B + 1)
d = rng.standard_normal(14) * 0.2 + 0.05
check("4. Monte Carlo p-value (B=20000) vs exact enumeration at N=14", abs(signflip_mc(d) - signflip_exact(d)), 4 * 0.5 / math.sqrt(20000))
# size under the null: A and B are two draws of the same system; per-problem pass@8 estimators from n=64 samples with problem-specific p_x
reject = 0; R = 1500; N_ = 60
for r in range(R):
    px = rng.beta(a_, b_, size=N_)
    d = lgk[rng.binomial(n, px)] - lgk[rng.binomial(n, px)]
    reject += signflip_mc(d, B=2000) <= 0.05
size = reject / R
check("4. size of the sign-flip test under H0 (N=60 problems, shared p_x, 1500 benchmarks)", abs(size - 0.05) , 4 * math.sqrt(0.05 * 0.95 / R) + 0.01, f"size {size:.4f}")
# the null is symmetry, not mean zero: a skewed mean-zero d inflates/deflates size
reject = 0; R = 1500; N_ = 12
for r in range(R):
    d = np.where(rng.random(N_) < 0.1, 9.0, -1.0)           # mean = 0.1*9 - 0.9 = 0, strongly skewed
    reject += signflip_exact(d) <= 0.05 if N_ <= 12 else 0
size_skew = reject / R
print(f"  skewed mean-zero differences (N=12): size = {size_skew:.4f} (nominal 0.05): the test's null is symmetry")
check("4. skewed mean-zero null: size differs from 0.05 by more than 0.02 (the test is not a test of the mean)", 0.0 if abs(size_skew - 0.05) > 0.02 else 1.0, 0.0)
# %% [markdown]
# ## 5. Fresh seed on the same problems removes generation noise, not selection bias
# K = 8 candidate checkpoints with *identical* true population quality theta (null).  On the fixed benchmark each has a quenched problem-specific advantage
# `beta_k ~ N(0, sC2/N)` (seed x problem) plus a fresh-generation error `g_k ~ N(0, se2/(N m))`.  Select `k* = argmax (beta_k + g_k)`.
# (i) report the same estimate: bias `sigma_tot c_K`;  (ii) re-sample generations on the SAME problems: bias `E beta_{k*} = (sigma_b^2/sigma_tot) c_K`;  (iii) fresh problems: bias 0.
# `c_K = E max of K standard normals`.
# %%
sb2 = sC2 / N; sg2 = se2 / (N * m); st2 = sb2 + sg2; K = 8
cK = float(np.mean(np.max(np.random.default_rng(3).standard_normal((2_000_000, K)), 1)))
from scipy import integrate
cK_exact = integrate.quad(lambda x: K * x * stats.norm.pdf(x) * stats.norm.cdf(x) ** (K - 1), -12, 12)[0]
check("5. c_8 = E max of 8 N(0,1): quadrature vs 2e6 Monte Carlo", abs(cK_exact - cK), 5 * 1.0 / math.sqrt(2e6) , f"c_8={cK_exact:.5f}")
rng = np.random.default_rng(4); R = 400_000
beta = math.sqrt(sb2) * rng.standard_normal((R, K)); g1 = math.sqrt(sg2) * rng.standard_normal((R, K)); g2 = math.sqrt(sg2) * rng.standard_normal((R, K))
win = np.argmax(beta + g1, axis=1); ix = np.arange(R)
b_same = (beta + g1)[ix, win]; b_regen = (beta + g2)[ix, win]; b_fresh = (math.sqrt(sb2) * rng.standard_normal(R) + math.sqrt(sg2) * rng.standard_normal(R))
th_same, th_regen = math.sqrt(st2) * cK_exact, (sb2 / math.sqrt(st2)) * cK_exact
sd = math.sqrt(st2)
check("5(i). same-sample winner estimate bias vs sigma_tot c_K", abs(b_same.mean() - th_same), 4 * sd / math.sqrt(R), f"{b_same.mean()*100:.3f} pp vs {th_same*100:.3f} pp")
check("5(ii). fresh generations on the same problems: bias vs (sb2/sigma_tot) c_K", abs(b_regen.mean() - th_regen), 4 * sd / math.sqrt(R), f"{b_regen.mean()*100:.3f} pp vs {th_regen*100:.3f} pp  (retained fraction {th_regen/th_same:.3f} = sb2/sigma_tot^2 = {sb2/st2:.3f})")
check("5(iii). fresh problems: bias 0", abs(b_fresh.mean()), 4 * sd / math.sqrt(R))
check("5. retained fraction of the bias after regenerating = sb2/(sb2+sg2) (with m=16: generation is a minor term)", abs(th_regen / th_same - sb2 / st2), 1e-12)
# %% [markdown]
# ## 6. Equivalence (TOST) versus "no significant difference"; minimum detectable effect
# `delta_hat ~ N(delta, s^2)`.  Equivalence margin `Delta`.  TOST: reject `H0: |delta| >= Delta` iff both one-sided tests reject: `p_TOST = max(p_L, p_U)`.
# Duality: TOST at level alpha rejects iff the (1-2 alpha) CI lies inside `(-Delta, Delta)`.  Power at `delta=0`: `2 Phi(Delta/s - z_{1-alpha}) - 1`.
# %%
al = 0.05; zA = stats.norm.ppf(1 - al)
def tost_p(dh, s, Dl):
    p_upper = stats.norm.cdf((dh - Dl) / s)          # H0: delta >= Delta
    p_lower = stats.norm.sf((dh + Dl) / s)           # H0: delta <= -Delta
    return np.maximum(p_lower, p_upper)
rng = np.random.default_rng(5)
dh = rng.normal(0, 0.03, 100000); s_ = rng.uniform(0.005, 0.03, 100000); Dl = 0.02
rej_p = tost_p(dh, s_, Dl) < al
rej_ci = (dh - zA * s_ > -Dl) & (dh + zA * s_ < Dl)
check("6. TOST p<alpha  <=>  90% CI inside (-Delta, Delta) (100000 random cases)", float((rej_p != rej_ci).sum()), 0.0)
# size at the boundary and power at delta = 0
R = 400_000
for dl in (Dl, -Dl):
    sm = 0.008; x = rng.normal(dl, sm, R); r = (tost_p(x, sm, Dl) < al).mean()
    check(f"6. size at delta={dl:+.2f} (boundary of H0), s=0.008: TOST rejection prob = alpha", abs(r - al), 4 * math.sqrt(al * (1 - al) / R), f"{r:.4f}")
sm = 0.008; x = rng.normal(0, sm, R); r = (tost_p(x, sm, Dl) < al).mean(); pw = 2 * stats.norm.cdf(Dl / sm - zA) - 1
check("6. power at delta=0: 2 Phi(Delta/s - z) - 1 vs MC", abs(r - pw), 4 * math.sqrt(0.25 / R), f"{r:.4f} vs {pw:.4f}")
# difference test vs equivalence at a true effect larger than the margin: N=244, discordance d=0.1
dpair = 0.1; sd_pair = math.sqrt(dpair); s244 = sd_pair / math.sqrt(244)
true_delta = 0.03
x = rng.normal(true_delta, s244, R)
p_diff = 2 * stats.norm.sf(np.abs(x) / s244)
miss = (p_diff >= 0.05).mean(); fake_equiv = (tost_p(x, s244, Dl) < al).mean()
print(f"  N=244, d=0.10 (s = {s244:.4f}), true delta=+0.03 > Delta=0.02:  difference test misses it {miss*100:.1f}% of the time; TOST declares equivalence {fake_equiv*100:.2f}%")
check("6. 'not significant' is not equivalence: with a true effect beyond the margin, difference test fails to reject often, TOST seldom 'accepts' (<= alpha)", 0.0 if (miss > 0.4 and fake_equiv <= al + 4 * math.sqrt(al * 0.95 / R)) else 1.0, 0.0)
# MDE (paired, one run per benchmark) and required N, cross-checked against note 0208's table (normal approximation, d = 0.10, delta = 5/244)
Nn = 244; mde = (stats.norm.ppf(0.975) + stats.norm.ppf(0.8)) * sd_pair / math.sqrt(Nn)
x = rng.normal(mde, s244, R); pw = (2 * stats.norm.sf(np.abs(x) / s244) < 0.05).mean()
check("6. MDE = (z_{.975}+z_{.8}) sigma_d/sqrt(N) gives 80% power (Gaussian paired mean)", abs(pw - 0.8), 4 * math.sqrt(0.16 / R) + 0.0035, f"MDE = {mde*Nn:.2f} problems, power {pw:.4f}  (0208 exact trinomial: 14.3)")
dlt = 5 / 244; d_ = 0.10
N_norm = ((stats.norm.ppf(0.975) * math.sqrt(d_) + stats.norm.ppf(0.8) * math.sqrt(d_ - dlt ** 2)) / dlt) ** 2
check("6. N for 80% power at delta=5/244, d=.10 equals note 0208's normal-approximation value 1867", abs(round(N_norm) - 1867), 1, f"N={N_norm:.1f}")
Req = (( stats.norm.ppf(1 - al) + stats.norm.ppf(0.9)) * sd_pair / Dl) ** 2        # TOST 80% power at delta=0:  Delta/s - z = z_{0.9}
pw_req = 2 * stats.norm.cdf(Dl / (sd_pair / math.sqrt(Req)) - zA) - 1
check("6. N for 80% TOST power at delta=0, Delta=0.02, d=.10: closed form reproduces 80% power", abs(pw_req - 0.8), 1e-9, f"N = {Req:.0f}")
print(f"  (equivalence within 2 points needs N ~ {Req:.0f} problems at 10% discordance; detecting a 5-problem gain needs ~1867)")
# %% [markdown]
# ## Summary
# %%
for n_, e_, t_ in CHECKS: print(f"PASS {n_[:105]}")
print(f"\n{len(CHECKS)} checks, all passed")

	# %% [markdown]
	# # 0209 - Variance components and equivalence tests for evaluation
	# Companion to `notes/0209 - Variance Components and Equivalence Tests.pdf`. CPU only, numpy/scipy, 1 thread.
	# Running toy: a benchmark of `N = 244` problems, `S` independently trained models ("seeds") of one recipe, `m = 16` samples per (seed, problem).
	# Per-(seed, problem) solve probability `p_sx = expit(mu_x + a_s + b_sx)`, `mu_x ~ N(0.3, 1.5^2)` (problem difficulty), `a_s ~ N(0, 0.25^2)` (seed), `b_sx ~ N(0, 0.5^2)` (seed x problem).
	# %%
	import os
	os.environ["OMP_NUM_THREADS"] = "1"; os.environ["OPENBLAS_NUM_THREADS"] = "1"; os.environ["MKL_NUM_THREADS"] = "1"
	import math, itertools
	from fractions import Fraction as Fr
	from math import comb
	import numpy as np
	from scipy import stats
	from scipy.special import expit
	CHECKS = []
	def check(name, err, tol, note=""):
	CHECKS.append((name, float(err), float(tol)))
	assert err <= tol, (name, float(err), tol)
	print(f"[ok] {name:100s} err={float(err):.3e} (tol {tol:.1e}) {note}")
	N, m = 244, 16
	M0, SX, SA, SB = 0.3, 1.5, 0.25, 0.5
	# %% [markdown]
	# ## 1. Variance components: problems (A), seeds (B), seed x problem (C), generation (e)
	# Exchangeable array `p_sx = mu + A_x + B_s + C_sx` (Hoeffding/ANOVA decomposition; components orthogonal). With outcomes `c_sx ~ Bin(m, p_sx)` and the score of one run `Y_s = (1/(N m)) sum_x c_sx`:
	# `Var(Y_s) = sA2/N + sB2 + sC2/N + se2/(N m)`, `se2 = E p(1-p)`. Only `sB2` survives `N -> inf` and `m -> inf`.
	# Population components by Gauss-Hermite integration of the logit-normal model (no simulation).
	# %%
	gh_x, gh_w = np.polynomial.hermite_e.hermegauss(60); gh_w = gh_w / gh_w.sum() # standard normal nodes/weights
	MU = M0 + SX * gh_x; A_ = SA * gh_x; B_ = SB * gh_x
	P3 = expit(MU[:, None, None] + A_[None, :, None] + B_[None, None, :]) # [x, s, b]
	W3 = gh_w[:, None, None] * gh_w[None, :, None] * gh_w[None, None, :]
	Ep = (P3 * W3).sum(); Ep2 = (P3 ** 2 * W3).sum()
	Ax = (P3 * gh_w[None, :, None] * gh_w[None, None, :]).sum((1, 2)) # E_{s,b} p \| x (weights over s,b sum to 1)
	Bs = (P3 * gh_w[:, None, None] * gh_w[None, None, :]).sum((0, 2)) # E_{x,b} p \| s
	sA2 = (gh_w * (Ax - Ep) ** 2).sum(); sB2 = (gh_w * (Bs - Ep) ** 2).sum()
	sT2 = Ep2 - Ep ** 2; sC2 = sT2 - sA2 - sB2; se2 = Ep - Ep2
	print(f" E p = {Ep:.4f}; sigma_A^2 = {sA2:.5f}, sigma_B^2 = {sB2:.6f} (sd {math.sqrt(sB2):.4f}), sigma_C^2 = {sC2:.5f}, sigma_e^2 = E p(1-p) = {se2:.5f}")
	def var_Y(S, N=N, m=m): return sA2 / N + sB2 / S + sC2 / (N * S) + se2 / (N * S * m) # variance of a grand mean over S seeds
	print(" per-term contribution to Var(Y_s), S=1: problems %.3e seed %.3e seedproblem %.3e generation %.3e" % (sA2 / N, sB2, sC2 / N, se2 / (N m)))
	check("1. components are non-negative and sum to Var(p)", abs(sA2 + sB2 + sC2 - sT2) + max(0, -min(sA2, sB2, sC2)), 1e-15)
	# simulate replications: each replication redraws problems, seeds, samples
	rng = np.random.default_rng(0)
	def simulate(R, S, N=N, m=m):
	mu = M0 + SX * rng.standard_normal((R, 1, N)); a = SA * rng.standard_normal((R, S, 1)); b = SB * rng.standard_normal((R, S, N))
	p = expit(mu + a + b)
	c = rng.binomial(m, p)
	return p, c
	R1 = 6000
	p, c = simulate(R1, 1)
	Y = c[:, 0, :].sum(1) / (N * m)
	se_rel = math.sqrt(2 / (R1 - 1))
	check("1. Var(Y_s) over 6000 replications vs sA2/N + sB2 + sC2/N + se2/(Nm)", abs(Y.var(ddof=1) / var_Y(1) - 1), 4 * se_rel, f"emp {Y.var(ddof=1):.4e} vs {var_Y(1):.4e}")
	# MoM / ANOVA estimators with S=4 seeds
	S = 4; R2 = 4000
	p, c = simulate(R2, S)
	ybar = c / m # cell means [R,S,N]
	gm = ybar.mean((1, 2), keepdims=True)
	MSx = S * ((ybar.mean(1) - gm[:, :, 0]) ** 2).sum(1) / (N - 1)
	MSs = N * ((ybar.mean(2) - gm[:, :, 0]) ** 2).sum(1) / (S - 1)
	res = ybar - ybar.mean(1, keepdims=True) - ybar.mean(2, keepdims=True) + gm
	MSr = (res ** 2).sum((1, 2)) / ((S - 1) * (N - 1))
	pq = (c * (m - c) / (m * (m - 1))).mean((1, 2)) # unbiased for E p(1-p)
	est = {"sigma_e^2": (pq, se2), "sigma_C^2": (MSr - pq / m, sC2), "sigma_A^2": ((MSx - MSr) / S, sA2), "sigma_B^2": ((MSs - MSr) / N, sB2)}
	for k, (v, t) in est.items():
	z = (v.mean() - t) / (v.std(ddof=1) / math.sqrt(R2))
	check(f"1. ANOVA estimator of {k} is unbiased (4000 reps, S=4)", abs(z), 4.0, f"mean {v.mean():.6f} vs true {t:.6f} (z={z:+.2f})")
	# %% [markdown]
	# ## 2. Which term do you pay down with what? Two recipes, S seeds each, same problems
	# `delta_hat = Ybar_A - Ybar_B`. Under H0 (same recipe) the shared problem effect A_x cancels: `Var(delta_hat) = (2/S) [ sB2 + sC2/N + se2/(N m) ]`.
	# The seed term does not shrink with N or m: a floor on the minimum detectable effect.
	# %%
	def var_diff(S, N=N, m=m): return (2 / S) * (sB2 + sC2 / N + se2 / (N * m))
	R3 = 6000
	for S in (1, 3):
	mu = M0 + SX * rng.standard_normal((R3, 1, N))
	d = []
	for _ in range(2):
	a = SA * rng.standard_normal((R3, S, 1)); b = SB * rng.standard_normal((R3, S, N))
	cc = rng.binomial(m, expit(mu + a + b)); d.append(cc.sum((1, 2)) / (S * N * m))
	dd = d[0] - d[1]
	check(f"2. Var(delta_hat) for S={S} seeds per recipe vs (2/S)[sB2 + sC2/N + se2/(Nm)]", abs(dd.var(ddof=1) / var_diff(S) - 1), 4 * math.sqrt(2 / (R3 - 1)), f"emp {dd.var(ddof=1):.3e} vs {var_diff(S):.3e}")
	z = stats.norm.ppf(0.975) + stats.norm.ppf(0.8)
	for S in (1, 3, 5):
	print(f" S={S}: MDE(80% power, 5% two-sided) at N=244: {zmath.sqrt(var_diff(S))100:.2f} pp; N->inf, m->inf floor: {zmath.sqrt(2sB2/S)*100:.2f} pp")
	check("2. floor is strictly positive: MDE(N=1e9, m=1e9, S=3) equals z sqrt(2 sB2 / 3)", abs(z * math.sqrt(var_diff(3, N=109, m=109)) - z * math.sqrt(2 * sB2 / 3)), 1e-6)
	# %% [markdown]
	# ## 3. Unbiased pass@k from n samples, and its exact variance
	# `e(n,c,k) = 1 - C(n-c,k)/C(n,k)`; `c ~ Bin(n,p)` => `E e = 1 - (1-p)^k`. Exact variance by summing over `c = 0..n`.
	# NAIVE plug-in `1 - (1 - c/n)^k` is biased. Special cases: `k = 1` gives `p(1-p)/n`; `k = n` gives `q(1-q)`, `q = 1-(1-p)^n`.
	# %%
	def e_nck(n, c, k): return 1 - Fr(comb(n - c, k), comb(n, k))
	def pk_moments(n, k, p):
	p = Fr(p); pm = [Fr(comb(n, c)) * p ** c * (1 - p) ** (n - c) for c in range(n + 1)]
	E1 = sum(w * e_nck(n, c, k) for c, w in enumerate(pm)); E2 = sum(w * e_nck(n, c, k) ** 2 for c, w in enumerate(pm))
	return E1, E2 - E1 ** 2, sum(w * (1 - (1 - Fr(c, n)) ** k) for c, w in enumerate(pm))
	worst = 0; worst_k1 = 0; worst_kn = 0
	for n, k, p in [(64, 8, Fr(1, 10)), (64, 32, Fr(1, 50)), (16, 4, Fr(3, 10)), (64, 1, Fr(1, 7)), (64, 64, Fr(1, 100)), (10, 10, Fr(1, 5))]:
	E1, V, Enaive = pk_moments(n, k, p)
	worst = max(worst, abs(float(E1 - (1 - (1 - p) ** k))))
	if k == 1: worst_k1 = max(worst_k1, abs(float(V - p * (1 - p) / n)))
	if k == n: q = 1 - (1 - p) ** n; worst_kn = max(worst_kn, abs(float(V - q * (1 - q))))
	check("3. E[e(n,c,k)] = 1-(1-p)^k exactly (6 (n,k,p) triples, rational arithmetic)", worst, 1e-15)
	check("3. k=1: Var = p(1-p)/n exactly", worst_k1, 1e-15)
	check("3. k=n: Var = q(1-q) with q = 1-(1-p)^n exactly", worst_kn, 1e-15)
	print(" n=64: p k true pass@k E[unbiased] E[naive plug-in] (bias) sd(unbiased) sd/true")
	for p in (0.02, 0.1, 0.5):
	for k in (1, 8, 32):
	E1, V, En = pk_moments(64, k, Fr(p).limit_denominator(1000)); g = 1 - (1 - p) ** k
	print(f" {p:5.2f} {k:4d} {g:.5f} {float(E1):.5f} {float(En):.5f} ({float(En)-g:+.5f}) {math.sqrt(float(V)):.5f} {math.sqrt(float(V))/g:.3f}")
	E1, V, En = pk_moments(64, 32, Fr(1, 50))
	check("3. naive plug-in is biased low at (n=64,k=32,p=.02) by more than 0.07", 0.0 if (1 - (1 - 0.02) 32) - float(En) > 0.07 else 1.0, 0.0, f"bias {float(En)-(1-0.9832):+.4f}")
	# MC check of the exact variance (independent implementation via binomial draws) and of the problem-level decomposition
	rng = np.random.default_rng(1)
	n, k, p0 = 64, 8, 0.1
	cs = rng.binomial(n, p0, size=2_000_000)
	lg = np.array([float(e_nck(n, c_, k)) for c_ in range(n + 1)])
	est = lg[cs]
	Vex = float(pk_moments(n, k, Fr(1, 10))[1])
	check("3. exact Var vs 2e6 Monte Carlo draws (n=64,k=8,p=.1)", abs(est.var() / Vex - 1), 4 * math.sqrt(2 / 2e6))
	# across problems: p_x ~ Beta(1.4, 0.7); Var(mean_x e_x) = (1/N)[Var_x(g(p_x)) + E_x Var(e \| p_x)]
	a_, b_ = 1.4, 0.7; k = 8; n = 64
	xs = rng.beta(a_, b_, size=400_000) # for the population quantities
	cs = rng.binomial(n, xs)
	lgk = np.array([float(e_nck(n, c_, k)) for c_ in range(n + 1)])
	es = lgk[cs]; gx = 1 - (1 - xs) ** k
	pred = (gx.var() + (es - gx).var()) / N
	reps = np.array([lgk[rng.binomial(n, rng.beta(a_, b_, size=N))].mean() for _ in range(6000)])
	check("3. Var of the benchmark mean pass@8 (N=244, n=64) vs (1/N)[Var_x g + E Var(e\|p)]", abs(reps.var(ddof=1) / pred - 1), 5 * math.sqrt(2 / 5999) + 0.01, f"emp {reps.var(ddof=1):.3e} vs {pred:.3e}")
	# %% [markdown]
	# ## 4. Paired sign-flip test on per-problem differences
	# `d_x = e^A_x - e^B_x`, statistic `T = sum_x d_x`, null: each `d_x` symmetric about 0 independently => signs exchangeable: p = fraction of the `2^N` sign vectors with `\|T'\| >= \|T\|`.
	# Binary `d_x in {-1,0,1}`: reduces to exact McNemar (note 0208). Exact enumeration for N<=18; Monte Carlo for N=244; size under the null.
	# %%
	def signflip_exact(d):
	d = np.asarray(d, float); n = len(d)
	signs = np.array(list(itertools.product([-1.0, 1.0], repeat=n)))
	T = np.abs(signs @ d); return float((T >= abs(d.sum()) - 1e-12).mean())
	def mcnemar(b, c):
	nd = b + c
	return 1.0 if nd == 0 else min(1.0, 2 * sum(comb(nd, i) for i in range(min(b, c) + 1)) / 2 ** nd)
	worst = 0
	for b_ in range(0, 8):
	for c_ in range(0, 8):
	if b_ + c_ == 0: continue
	d = np.array([1] * b_ + [-1] * c_ + [0] * 3)
	worst = max(worst, abs(signflip_exact(d) - mcnemar(b_, c_)))
	check("4. sign-flip p on {-1,0,+1} equals exact McNemar for all b,c <= 7 (3 concordant zeros)", worst, 1e-12)
	check("4. (5,0): two-sided sign-flip p = 0.0625", abs(signflip_exact([1] * 5) - 0.0625), 1e-15)
	rng = np.random.default_rng(2)
	def signflip_mc(d, B=20000):
	d = np.asarray(d, float); s = rng.choice([-1.0, 1.0], size=(B, len(d))); T = np.abs(s @ d)
	return (1 + (T >= abs(d.sum()) - 1e-12).sum()) / (B + 1)
	d = rng.standard_normal(14) * 0.2 + 0.05
	check("4. Monte Carlo p-value (B=20000) vs exact enumeration at N=14", abs(signflip_mc(d) - signflip_exact(d)), 4 * 0.5 / math.sqrt(20000))
	# size under the null: A and B are two draws of the same system; per-problem pass@8 estimators from n=64 samples with problem-specific p_x
	reject = 0; R = 1500; N_ = 60
	for r in range(R):
	px = rng.beta(a_, b_, size=N_)
	d = lgk[rng.binomial(n, px)] - lgk[rng.binomial(n, px)]
	reject += signflip_mc(d, B=2000) <= 0.05
	size = reject / R
	check("4. size of the sign-flip test under H0 (N=60 problems, shared p_x, 1500 benchmarks)", abs(size - 0.05) , 4 * math.sqrt(0.05 * 0.95 / R) + 0.01, f"size {size:.4f}")
	# the null is symmetry, not mean zero: a skewed mean-zero d inflates/deflates size
	reject = 0; R = 1500; N_ = 12
	for r in range(R):
	d = np.where(rng.random(N_) < 0.1, 9.0, -1.0) # mean = 0.1*9 - 0.9 = 0, strongly skewed
	reject += signflip_exact(d) <= 0.05 if N_ <= 12 else 0
	size_skew = reject / R
	print(f" skewed mean-zero differences (N=12): size = {size_skew:.4f} (nominal 0.05): the test's null is symmetry")
	check("4. skewed mean-zero null: size differs from 0.05 by more than 0.02 (the test is not a test of the mean)", 0.0 if abs(size_skew - 0.05) > 0.02 else 1.0, 0.0)
	# %% [markdown]
	# ## 5. Fresh seed on the same problems removes generation noise, not selection bias
	# K = 8 candidate checkpoints with identical true population quality theta (null). On the fixed benchmark each has a quenched problem-specific advantage
	# `beta_k ~ N(0, sC2/N)` (seed x problem) plus a fresh-generation error `g_k ~ N(0, se2/(N m))`. Select `k* = argmax (beta_k + g_k)`.
	# (i) report the same estimate: bias `sigma_tot c_K`; (ii) re-sample generations on the SAME problems: bias `E beta_{k*} = (sigma_b^2/sigma_tot) c_K`; (iii) fresh problems: bias 0.
	# `c_K = E max of K standard normals`.
	# %%
	sb2 = sC2 / N; sg2 = se2 / (N * m); st2 = sb2 + sg2; K = 8
	cK = float(np.mean(np.max(np.random.default_rng(3).standard_normal((2_000_000, K)), 1)))
	from scipy import integrate
	cK_exact = integrate.quad(lambda x: K * x * stats.norm.pdf(x) * stats.norm.cdf(x) ** (K - 1), -12, 12)[0]
	check("5. c_8 = E max of 8 N(0,1): quadrature vs 2e6 Monte Carlo", abs(cK_exact - cK), 5 * 1.0 / math.sqrt(2e6) , f"c_8={cK_exact:.5f}")
	rng = np.random.default_rng(4); R = 400_000
	beta = math.sqrt(sb2) * rng.standard_normal((R, K)); g1 = math.sqrt(sg2) * rng.standard_normal((R, K)); g2 = math.sqrt(sg2) * rng.standard_normal((R, K))
	win = np.argmax(beta + g1, axis=1); ix = np.arange(R)
	b_same = (beta + g1)[ix, win]; b_regen = (beta + g2)[ix, win]; b_fresh = (math.sqrt(sb2) * rng.standard_normal(R) + math.sqrt(sg2) * rng.standard_normal(R))
	th_same, th_regen = math.sqrt(st2) * cK_exact, (sb2 / math.sqrt(st2)) * cK_exact
	sd = math.sqrt(st2)
	check("5(i). same-sample winner estimate bias vs sigma_tot c_K", abs(b_same.mean() - th_same), 4 * sd / math.sqrt(R), f"{b_same.mean()100:.3f} pp vs {th_same100:.3f} pp")
	check("5(ii). fresh generations on the same problems: bias vs (sb2/sigma_tot) c_K", abs(b_regen.mean() - th_regen), 4 * sd / math.sqrt(R), f"{b_regen.mean()100:.3f} pp vs {th_regen100:.3f} pp (retained fraction {th_regen/th_same:.3f} = sb2/sigma_tot^2 = {sb2/st2:.3f})")
	check("5(iii). fresh problems: bias 0", abs(b_fresh.mean()), 4 * sd / math.sqrt(R))
	check("5. retained fraction of the bias after regenerating = sb2/(sb2+sg2) (with m=16: generation is a minor term)", abs(th_regen / th_same - sb2 / st2), 1e-12)
	# %% [markdown]
	# ## 6. Equivalence (TOST) versus "no significant difference"; minimum detectable effect
	# `delta_hat ~ N(delta, s^2)`. Equivalence margin `Delta`. TOST: reject `H0: \|delta\| >= Delta` iff both one-sided tests reject: `p_TOST = max(p_L, p_U)`.
	# Duality: TOST at level alpha rejects iff the (1-2 alpha) CI lies inside `(-Delta, Delta)`. Power at `delta=0`: `2 Phi(Delta/s - z_{1-alpha}) - 1`.
	# %%
	al = 0.05; zA = stats.norm.ppf(1 - al)
	def tost_p(dh, s, Dl):
	p_upper = stats.norm.cdf((dh - Dl) / s) # H0: delta >= Delta
	p_lower = stats.norm.sf((dh + Dl) / s) # H0: delta <= -Delta
	return np.maximum(p_lower, p_upper)
	rng = np.random.default_rng(5)
	dh = rng.normal(0, 0.03, 100000); s_ = rng.uniform(0.005, 0.03, 100000); Dl = 0.02
	rej_p = tost_p(dh, s_, Dl) < al
	rej_ci = (dh - zA * s_ > -Dl) & (dh + zA * s_ < Dl)
	check("6. TOST p<alpha <=> 90% CI inside (-Delta, Delta) (100000 random cases)", float((rej_p != rej_ci).sum()), 0.0)
	# size at the boundary and power at delta = 0
	R = 400_000
	for dl in (Dl, -Dl):
	sm = 0.008; x = rng.normal(dl, sm, R); r = (tost_p(x, sm, Dl) < al).mean()
	check(f"6. size at delta={dl:+.2f} (boundary of H0), s=0.008: TOST rejection prob = alpha", abs(r - al), 4 * math.sqrt(al * (1 - al) / R), f"{r:.4f}")
	sm = 0.008; x = rng.normal(0, sm, R); r = (tost_p(x, sm, Dl) < al).mean(); pw = 2 * stats.norm.cdf(Dl / sm - zA) - 1
	check("6. power at delta=0: 2 Phi(Delta/s - z) - 1 vs MC", abs(r - pw), 4 * math.sqrt(0.25 / R), f"{r:.4f} vs {pw:.4f}")
	# difference test vs equivalence at a true effect larger than the margin: N=244, discordance d=0.1
	dpair = 0.1; sd_pair = math.sqrt(dpair); s244 = sd_pair / math.sqrt(244)
	true_delta = 0.03
	x = rng.normal(true_delta, s244, R)
	p_diff = 2 * stats.norm.sf(np.abs(x) / s244)
	miss = (p_diff >= 0.05).mean(); fake_equiv = (tost_p(x, s244, Dl) < al).mean()
	print(f" N=244, d=0.10 (s = {s244:.4f}), true delta=+0.03 > Delta=0.02: difference test misses it {miss100:.1f}% of the time; TOST declares equivalence {fake_equiv100:.2f}%")
	check("6. 'not significant' is not equivalence: with a true effect beyond the margin, difference test fails to reject often, TOST seldom 'accepts' (<= alpha)", 0.0 if (miss > 0.4 and fake_equiv <= al + 4 * math.sqrt(al * 0.95 / R)) else 1.0, 0.0)
	# MDE (paired, one run per benchmark) and required N, cross-checked against note 0208's table (normal approximation, d = 0.10, delta = 5/244)
	Nn = 244; mde = (stats.norm.ppf(0.975) + stats.norm.ppf(0.8)) * sd_pair / math.sqrt(Nn)
	x = rng.normal(mde, s244, R); pw = (2 * stats.norm.sf(np.abs(x) / s244) < 0.05).mean()
	check("6. MDE = (z_{.975}+z_{.8}) sigma_d/sqrt(N) gives 80% power (Gaussian paired mean)", abs(pw - 0.8), 4 * math.sqrt(0.16 / R) + 0.0035, f"MDE = {mde*Nn:.2f} problems, power {pw:.4f} (0208 exact trinomial: 14.3)")
	dlt = 5 / 244; d_ = 0.10
	N_norm = ((stats.norm.ppf(0.975) * math.sqrt(d_) + stats.norm.ppf(0.8) * math.sqrt(d_ - dlt 2)) / dlt) 2
	check("6. N for 80% power at delta=5/244, d=.10 equals note 0208's normal-approximation value 1867", abs(round(N_norm) - 1867), 1, f"N={N_norm:.1f}")
	Req = (( stats.norm.ppf(1 - al) + stats.norm.ppf(0.9)) * sd_pair / Dl) ** 2 # TOST 80% power at delta=0: Delta/s - z = z_{0.9}
	pw_req = 2 * stats.norm.cdf(Dl / (sd_pair / math.sqrt(Req)) - zA) - 1
	check("6. N for 80% TOST power at delta=0, Delta=0.02, d=.10: closed form reproduces 80% power", abs(pw_req - 0.8), 1e-9, f"N = {Req:.0f}")
	print(f" (equivalence within 2 points needs N ~ {Req:.0f} problems at 10% discordance; detecting a 5-problem gain needs ~1867)")
	# %% [markdown]
	# ## Summary
	# %%
	for n_, e_, t_ in CHECKS: print(f"PASS {n_[:105]}")
	print(f"\n{len(CHECKS)} checks, all passed")