1# Verification checks for P4 (reads out/results.json, out/table1.csv, data/design.csv)2import json, numpy as np, pandas as pd3r = json.load(open('out/results.json')); t = pd.read_csv('out/table1.csv'); d = pd.read_csv('data/design.csv')4ok = lambda c, m: print(("PASS " if c else "FAIL ") + m)5ok(len(t) == 27 and set(t.estimator) == {'cr0','cr1','wcb'}, "table1 has 27 rows, 3 estimators")6ok(len(r['coverage']) == 27 and len(r['mc_se_coverage']) == 27 and len(r['rel_bias_var']) == 18, "key counts 27/27/18")7ok(t[t.estimator=='wcb'].rel_bias_var.isna().all(), "rel_bias_var empty for wcb")8ok(all(0 <= v <= 1 for v in r['coverage'].values()), "coverages in [0,1]")9ok(all(v > 0 for v in r['mc_se_coverage'].values()), "MC SEs positive")10ok(all(abs(r['mc_se_coverage'][k] - np.sqrt(p*(1-p)/r['R'])) < 1e-12 for k, p in r['coverage'].items()), "MC SE = sqrt(p(1-p)/R)")11ok(all(abs(v*r['R'] - round(v*r['R'])) < 1e-9 for v in r['coverage'].values()), "coverage*R integer")12f = sorted(k for k, v in r['coverage'].items() if v < 0.9)13ok(f == sorted(r['failure_cells']) and r['n_failure_cells'] == len(f), "failure_cells consistent")14nm = sorted(k for k, v in r['coverage'].items() if abs(v-0.95) <= 2*r['mc_se_coverage'][k])15ok(nm == sorted(r['nominal_cells']), "nominal_cells consistent")16ok(r['n_analysed'] == r['n_cells']*r['R'] == 18000 and r['R']==2000 and r['B']==399, "n_analysed=9*2000, R,B unchanged")17ok(r['random_seed'] == int(d.seed[d.cell_id=='C1'].iloc[0]), "random_seed = C1 seed")18ex = pd.read_csv('out/exclusions.csv'); ok(18000 - ex.n_removed.sum() == r['n_analysed'], "exclusion ledger sums")19# CR1/CR0 ratio must be G/(G-1)20ok(all(abs(r['rel_bias_var']['cr1'+k[3:]]/r['rel_bias_var'][k] - (int(k.split('_')[1][1:])/(int(k.split('_')[1][1:])-1))) < 1e-921 for k in r['rel_bias_var'] if k.startswith('cr0')), "rel_bias cr1/cr0 = G/(G-1)")22# analytic E[V_CR0]/Var(ybar) vs MC estimate23print("Analytic vs MC relative bias of CR0 (z = (MC-exact)/MCSE):")24for _, row in d.iterrows():25 n = np.array([int(x) for x in row.cluster_sizes.split('|')], float); N = n.sum(); icc = row.icc26 V = (icc*np.sum(n**2) + (1-icc)*N)/N**227 VT = icc*n**2 + (1-icc)*n28 EV0 = np.sum(VT*(1-2*n/N) + n**2*V)/N**229 k = f"cr0_G{row.G}_icc{icc:.1f}"30 z = (r['rel_bias_var'][k] - EV0/V)/r['mc_se_rel_bias_var'][k]31 ok(abs(z) < 3, f" {k}: exact {EV0/V:.4f} MC {r['rel_bias_var'][k]:.4f} z={z:.2f}")32 dg = r['cell_diagnostics'][row.cell_id]33 ok(abs(dg['var_exact'] - V) < 1e-15, f" {row.cell_id} var_exact matches formula")34# monotonicity sanity: cr1 halfwidth > cr0 halfwidth35w = t.pivot(index='cell_id', columns='estimator', values='mean_ci_halfwidth')36ok((w.cr1 > w.cr0).all(), "CR1 intervals wider than CR0 in every cell")37