1"""T2-simulation: finite-sample coverage of CR0, CR1 and wild-cluster-bootstrap-t intervals for a mean.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/design.csv; writes out/results.json, out/table1.csv, out/figure1.png.4RNG: per cell numpy default_rng(seed from design.csv); draw order per cell:5 u (R x G) ~ N(0, icc), e (R x N) ~ N(0, 1-icc), then Rademacher weights (B x G per replication)6 drawn chunk-wise (chunks of 200 replications, in replication order).7"""8import json, os9import numpy as np, pandas as pd10from scipy import stats11import matplotlib, matplotlib.ticker12matplotlib.use("Agg")13import matplotlib.pyplot as plt1415ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))16OUT = os.path.join(ROOT, "out")17design = pd.read_csv(os.path.join(ROOT, "data", "design.csv"))18CHUNK = 20019Z = 1.959963984540054 # 0.975 normal quantile ("1.96")202122def simulate_cell(row):23 G, icc, N, R, B, mu = int(row.G), float(row.icc), int(row.N), int(row.R), int(row.B), float(row.mu)24 n = np.array([int(x) for x in row.cluster_sizes.split("|")])25 assert len(n) == G and n.sum() == N26 gid = np.repeat(np.arange(G), n)27 var_exact = (icc * np.sum(n.astype(float) ** 2) + (1 - icc) * N) / N ** 228 rng = np.random.default_rng(int(row.seed))29 u = rng.normal(0.0, np.sqrt(icc), size=(R, G))30 e = rng.normal(0.0, np.sqrt(1 - icc), size=(R, N))31 y = mu + u[:, gid] + e # R x N32 ybar = y.mean(axis=1) # R33 # cluster sums of residuals S_g = sum_i (y_ig - ybar)34 csum = np.zeros((R, G))35 for g in range(G):36 csum[:, g] = y[:, gid == g].sum(axis=1)37 S = csum - n[None, :] * ybar[:, None] # R x G38 v0 = (S ** 2).sum(axis=1) / N ** 239 c = G / (G - 1)40 v1 = c * v041 # naive iid baseline (ignores clustering): s^2/N, normal critical value42 s2 = y.var(axis=1, ddof=1)43 vnaive = s2 / N44 # wild cluster bootstrap-t (unrestricted, Rademacher)45 k = int(np.ceil(0.95 * B)) # order statistic: smallest x with ECDF(x) >= 0.95 (=380th of 399)46 qwcb = np.empty(R)47 min_vstar = np.inf48 for start in range(0, R, CHUNK):49 stop = min(start + CHUNK, R)50 W = rng.integers(0, 2, size=(stop - start, B, G)).astype(np.float64) * 2 - 151 A = W * S[start:stop, None, :] # w_g * S_g : bootstrap cluster sums of (y* - ybar)52 D = A.sum(axis=2) / N # ybar* - ybar53 Sstar = A - n[None, None, :] * D[:, :, None] # cluster sums of bootstrap residuals54 vstar = c * (Sstar ** 2).sum(axis=2) / N ** 255 min_vstar = min(min_vstar, float(vstar.min()))56 tabs = np.abs(D) / np.sqrt(vstar)57 qwcb[start:stop] = np.sort(tabs, axis=1)[:, k - 1]58 tcrit = stats.t.ppf(0.975, G - 1)59 hw = {"cr0": Z * np.sqrt(v0), "cr1": tcrit * np.sqrt(v1), "wcb": qwcb * np.sqrt(v1),60 "naive": Z * np.sqrt(vnaive)}61 err = np.abs(ybar - mu)62 res = {}63 for est, h in hw.items():64 p = float(np.mean(err <= h))65 res[est] = dict(coverage=p, mc_se=float(np.sqrt(p * (1 - p) / R)), mean_hw=float(h.mean()))66 res["cr0"]["rel_bias"] = float(v0.mean() / var_exact)67 res["cr1"]["rel_bias"] = float(v1.mean() / var_exact)68 res["cr0"]["rel_bias_mc_se"] = float(v0.std(ddof=1) / np.sqrt(R) / var_exact)69 res["cr1"]["rel_bias_mc_se"] = float(v1.std(ddof=1) / np.sqrt(R) / var_exact)70 # analytic E[V_CR0]: S_g = T_g - n_g*ybar, Var(T_g)=s2g, Cov(T_g,ybar)=s2g/N71 s2g = icc * n.astype(float) ** 2 + (1 - icc) * n72 exact_rb0 = float(np.sum(s2g - 2 * n * s2g / N + n.astype(float) ** 2 * var_exact) / N ** 2 / var_exact)73 res["_diag"] = dict(var_exact=var_exact, emp_var_ybar=float(ybar.var(ddof=1)),74 min_vstar=min_vstar, exact_rel_bias_cr0=exact_rb0, exact_rel_bias_cr1=exact_rb0 * c, mean_q_wcb=float(qwcb.mean()), tcrit=float(tcrit),75 n_distinct_q=int(len(np.unique(np.round(qwcb, 10)))))76 return res777879def key(est, G, icc):80 return f"{est}_G{int(G)}_icc{float(icc):.1f}"818283rows, cov, se, rb, rbse = [], {}, {}, {}, {}84base_cov, base_se, diag = {}, {}, {}85for _, r in design.iterrows():86 res = simulate_cell(r)87 for est in ["cr0", "cr1", "wcb"]:88 kk = key(est, r.G, r.icc)89 cov[kk] = res[est]["coverage"]; se[kk] = res[est]["mc_se"]90 if est != "wcb":91 rb[kk] = res[est]["rel_bias"]; rbse[kk] = res[est]["rel_bias_mc_se"]92 rows.append(dict(estimator=est, cell_id=r.cell_id, G=int(r.G), icc=float(r.icc),93 coverage=res[est]["coverage"], mc_se_coverage=res[est]["mc_se"],94 rel_bias_var=res[est].get("rel_bias", np.nan),95 mean_ci_halfwidth=res[est]["mean_hw"], R=int(r.R)))96 kb = key("naive", r.G, r.icc)97 base_cov[kb] = res["naive"]["coverage"]; base_se[kb] = res["naive"]["mc_se"]98 diag[r.cell_id] = res["_diag"]99 print(r.cell_id, {e: round(res[e]["coverage"], 4) for e in ["naive", "cr0", "cr1", "wcb"]})100101table = pd.DataFrame(rows)102table.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")103104failure = sorted([k for k, p in cov.items() if p < 0.90])105nominal = sorted([k for k, p in cov.items() if abs(p - 0.95) <= 2 * se[k]])106smallG = {est: float(np.mean([cov[key(est, 5, i)] for i in (0.0, 0.2, 0.5)])) for est in ["cr0", "cr1", "wcb"]}107best = max(smallG, key=smallG.get)108R = int(design.R.iloc[0]); B = int(design.B.iloc[0])109results = {110 "task_id": "T2-simulation",111 "R": R, "B": B, "n_cells": int(len(design)),112 "coverage": cov, "mc_se_coverage": se, "rel_bias_var": rb,113 "failure_cells": failure, "n_failure_cells": len(failure),114 "nominal_cells": nominal,115 "best_procedure_small_G": best,116 "feasibility": "complete", "feasibility_reason": None,117 "data_issues": [118 {"issue": "ICC=0 cells: cluster effect u_g has variance 0 (degenerate normal).",119 "handling": "Drawn with numpy normal(scale=0), which returns exact zeros; DGP is then i.i.d. as intended."},120 {"issue": "G=5 with Rademacher weights: only 2^5=32 distinct weight vectors (16 distinct |t*| values), so B=399 draws resample a tiny discrete set and the WCB critical value is coarse.",121 "handling": "Kept as specified (task forbids changing B or weights); documented as a limitation."},122 ],123 "n_input_rows": {"design.csv": int(len(design))},124 "n_analysed": int(len(design) * R),125 "random_seed": int(design.loc[design.cell_id == "C1", "seed"].iloc[0]),126 "rel_bias_var_mc_se": rbse,127 "baseline_naive_iid_coverage": base_cov,128 "baseline_naive_iid_mc_se": base_se,129 "avg_coverage_small_G": smallG,130 "wcb_quantile_rule": "order statistic ceil(0.95*B) of sorted |t*| (380th of 399)",131 "diagnostics": diag,132}133with open(os.path.join(OUT, "results.json"), "w") as f:134 json.dump(results, f, indent=1, sort_keys=False)135136# figure137fig, axes = plt.subplots(1, 3, figsize=(13, 4.2), sharey=True)138style = {"cr0": ("o-", "CR0 + N(0,1)"), "cr1": ("s-", "CR1 + t(G-1)"), "wcb": ("^-", "WCB-t (Rademacher)")}139Gs = [5, 10, 30]140for ax, icc in zip(axes, [0.0, 0.2, 0.5]):141 for est, (m, lab) in style.items():142 p = np.array([cov[key(est, g, icc)] for g in Gs]); s = np.array([se[key(est, g, icc)] for g in Gs])143 ax.errorbar(Gs, p, yerr=1.96 * s, fmt=m, capsize=3, label=lab)144 ax.axhline(0.95, color="k", ls="--", lw=1, label="nominal 0.95")145 ax.axhline(0.90, color="r", ls=":", lw=1, label="failure line 0.90")146 ax.set_xscale("log"); ax.set_xticks(Gs); ax.set_xticklabels(Gs); ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())147 ax.set_title(f"ICC = {icc}"); ax.set_xlabel("number of clusters G")148axes[0].set_ylabel("coverage (bars: ±1.96 MC SE)")149axes[0].legend(fontsize=8, loc="lower right")150fig.suptitle(f"Coverage of nominal 95% intervals for the mean (R={R}, B={B})")151fig.tight_layout()152fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)153print("failure cells:", failure)154print("best procedure at G=5:", best, smallG)155