1#!/usr/bin/env python32"""Aggregate raw Monte Carlo CSVs into results/processed/summary.csv and draw the paper figures (PDF, vector).3Every number in the paper's Monte Carlo tables/figures is produced here from results/raw/*.csv (no manual entry).4Usage: python3 experiments/process_results.py [--figures]5"""6import argparse, glob, json, pathlib, re, sys7import numpy as np, pandas as pd8import matplotlib; matplotlib.use("Agg")9import matplotlib.pyplot as plt10from scipy import stats1112ROOT = pathlib.Path(__file__).resolve().parents[1]13RAW, PROC, FIG = ROOT / "results" / "raw", ROOT / "results" / "processed", ROOT / "results" / "figures"14PAL = {"blue": "#2a78d6", "orange": "#eb6834", "aqua": "#1baf7a", "yellow": "#eda100", "magenta": "#e87ba4", "green": "#008300", "violet": "#6250d6", "red": "#e34948"}15INK, INK2, GRID = "#0b0b0b", "#52514e", "#e6e5e1"16STYLE = { # fixed assignment estimator -> (color, marker, linestyle, label)17 "r1_q0": ("blue", "o", "-", "raw average (r=1, q=0)"), "r1_q1": ("blue", "o", "--", "raw + 1 Newton (1,1)"), "r1_q2": ("blue", "o", ":", "raw + 2 Newton (1,2)"),18 "r2_q0": ("orange", "s", "-", "half-block jackknife (2,0)"), "r2_q1": ("orange", "s", "--", "half-block + 1 Newton (2,1)"), "r2_q2": ("orange", "s", ":", "half-block + 2 Newton (2,2)"),19 "r3_q0": ("aqua", "^", "-", "quarter-block (3,0)"), "r3_q1": ("aqua", "^", "--", "quarter-block + 1 Newton (3,1)"), "r3_q2": ("aqua", "^", ":", "quarter-block + 2 Newton (3,2)"),20 "cs_q0": ("violet", "D", "-", "analytic Cox–Snell (q=0)"), "cs_q1": ("violet", "D", "--", "analytic Cox–Snell + 1 Newton"),21 "csl1": ("green", "v", "-", "CSL, 1 round"), "csl2": ("green", "v", "--", "CSL, 2 rounds"), "aee": ("yellow", "P", "-", "AEE (Lin–Xi)"),22 "stale_r1_q1": ("magenta", "x", "--", "raw + 1 stale-Hessian step"), "stale_r2_q1": ("magenta", "x", ":", "half-block + 1 stale-Hessian step"),23}24ALPHA = {"r1_q0": 0.5, "r1_q1": 0.75, "r1_q2": 0.875, "r2_q0": 0.75, "r2_q1": 0.875, "r2_q2": 0.9375, "r3_q0": 5 / 6, "r3_q1": 11 / 12, "r3_q2": 23 / 24}252627def _setup():28 plt.rcParams.update({"font.family": "serif", "font.size": 9, "axes.edgecolor": INK2, "axes.labelcolor": INK, "xtick.color": INK2, "ytick.color": INK2,29 "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6, "axes.spines.top": False, "axes.spines.right": False,30 "legend.frameon": False, "legend.fontsize": 8, "figure.dpi": 120, "savefig.bbox": "tight", "lines.linewidth": 1.6, "lines.markersize": 4.5,31 "axes.facecolor": "#fcfcfb", "figure.facecolor": "white"})323334def estimators(df):35 names = set()36 for c in df.columns:37 m = re.match(r"^(r\d_q\d|stale_r\d_q\d|cs_q\d|csl\d|aee)_0$", c)38 if m:39 names.add(m.group(1))40 return sorted(names)414243def summarize_cell(df):44 p = int(df["p"].iloc[0]); N = int(df["N"].iloc[0]); sN = np.sqrt(N)45 th0 = df[[f"theta0_{j}" for j in range(p)]].to_numpy()46 cen = df[[f"central_{j}" for j in range(p)]].to_numpy()47 rows = []48 B = len(df)49 for e in estimators(df):50 th = df[[f"{e}_{j}" for j in range(p)]].to_numpy()51 ok = np.all(np.isfinite(th), axis=1)52 gap = sN * (th - cen); err = sN * (th - th0)53 gap_n = np.linalg.norm(gap, axis=1); err_n = np.linalg.norm(err, axis=1)54 row = {"estimator": e, "reps": B, "finite_frac": ok.mean(),55 "mean_gap0": np.nanmean(gap[ok, 0]), "se_gap0": np.nanstd(gap[ok, 0], ddof=1) / np.sqrt(ok.sum()),56 "rms_gap": np.sqrt(np.nanmean(gap_n[ok] ** 2)), "median_gap": np.nanmedian(gap_n[ok]), "p95_gap": np.nanquantile(gap_n[ok], 0.95),57 "rmse": np.sqrt(np.nanmean(err_n[ok] ** 2)), "rmse_se": np.sqrt(np.nanvar(err_n[ok] ** 2, ddof=1) / ok.sum()) / (2 * max(np.sqrt(np.nanmean(err_n[ok] ** 2)), 1e-12)),58 "bias0": np.nanmean(err[ok, 0]), "bias0_se": np.nanstd(err[ok, 0], ddof=1) / np.sqrt(ok.sum()), "sd0": np.nanstd(err[ok, 0], ddof=1),59 "max_abs_bias": np.max(np.abs(np.nanmean(err[ok], axis=0)))}60 if f"se_{e}_0" in df.columns:61 se = df[[f"se_{e}_{j}" for j in range(p)]].to_numpy()62 for lv in (0.90, 0.95, 0.99):63 z = stats.norm.ppf(0.5 + lv / 2)64 hit = np.abs(th - th0) <= z * se65 row[f"cov{int(lv*100)}"] = np.nanmean(hit[ok]); row[f"len{int(lv*100)}"] = np.nanmean(2 * z * se[ok]) * sN66 if f"se_sand_central_0" in df.columns and e == "central":67 ses = df[[f"se_sand_central_{j}" for j in range(p)]].to_numpy()68 z = stats.norm.ppf(0.975); row["cov95_sandwich"] = np.nanmean((np.abs(th - th0) <= z * ses)[ok])69 rows.append(row)70 # centralized reference row71 err = sN * (cen - th0); rows.append({"estimator": "central", "reps": B, "finite_frac": 1.0, "rmse": np.sqrt(np.mean(np.sum(err ** 2, axis=1))), "bias0": err[:, 0].mean(),72 "sd0": err[:, 0].std(ddof=1), "mean_gap0": 0.0, "rms_gap": 0.0})73 if "se_central_0" in df.columns:74 se = df[[f"se_central_{j}" for j in range(p)]].to_numpy()75 for lv in (0.90, 0.95, 0.99):76 z = stats.norm.ppf(0.5 + lv / 2); rows[-1][f"cov{int(lv*100)}"] = np.mean(np.abs(cen - th0) <= z * se)77 if "se_sand_central_0" in df.columns:78 ses = df[[f"se_sand_central_{j}" for j in range(p)]].to_numpy(); rows[-1]["cov95_sandwich"] = np.mean(np.abs(cen - th0) <= stats.norm.ppf(0.975) * ses)79 out = pd.DataFrame(rows)80 for c in ("label", "model", "p", "N", "K", "m"):81 out[c] = df[c].iloc[0]82 for r in (1, 2, 3):83 if f"valid_r{r}" in df.columns:84 out[f"valid_frac_r{r}"] = df[f"valid_r{r}"].mean() / df["K"].iloc[0]85 for s in (1, 2, 4):86 if f"local_fail_s{s}" in df.columns:87 out[f"local_fail_frac_s{s}"] = df[f"local_fail_s{s}"].mean() / (df["K"].iloc[0] * s)88 out["central_fail"] = (df["central_status"] != 0).mean()89 out["t_estimate_s"] = df["t_estimate_s"].mean()90 for c in df.columns:91 if c.startswith("cond_r2_q1"):92 out["cond_r2_q1"] = df[c].mean()93 return out949596def build_summary():97 frames = []98 for f in sorted(glob.glob(str(RAW / "*.csv"))):99 if "boundary" in f or "ablation" in f or "skew" in f:100 continue101 df = pd.read_csv(f)102 if len(df) < 10:103 continue104 frames.append(summarize_cell(df))105 s = pd.concat(frames, ignore_index=True)106 PROC.mkdir(parents=True, exist_ok=True); s.to_csv(PROC / "summary.csv", index=False)107 return s108109110def fig_rmse_vs_m(s, label, p, N, metric="rms_gap", ests=("r1_q0", "r1_q1", "r1_q2", "r2_q0", "r2_q1", "r3_q0", "r3_q1", "cs_q1", "csl1"), fname=None, ylabel=None):111 d = s[(s.label == label) & (s.p == p) & (s.N == N)]112 if d.empty:113 return114 fig, ax = plt.subplots(figsize=(4.6, 3.4))115 for e in ests:116 dd = d[d.estimator == e].sort_values("m")117 if dd.empty or e not in STYLE:118 continue119 col, mk, ls, lab = STYLE[e]120 ax.plot(dd["m"], dd[metric].clip(lower=1e-12), marker=mk, ls=ls, color=PAL[col], label=lab, markerfacecolor="white", markeredgewidth=1.2)121 ax.set_xscale("log", base=2); ax.set_yscale("log"); ax.set_xlabel("local sample size m (N fixed)"); ax.set_ylabel(ylabel or metric)122 ax.set_ylim(bottom=max(ax.get_ylim()[0], 1e-8))123 ax.set_title(f"{label}, p = {p}, N = {N:,}", loc="left", fontsize=9, color=INK)124 ax.legend(ncol=2, loc="upper center", bbox_to_anchor=(0.5, -0.22), fontsize=7); FIG.mkdir(parents=True, exist_ok=True)125 fig.savefig(FIG / (fname or f"fig_{metric}_vs_m_{label}_p{p}_N{N}.pdf")); plt.close(fig)126127128def fig_gap_collapse(s, label, p, ests=("r1_q1", "r2_q1", "r1_q2", "r3_q1"), fname=None):129 """|mean scaled gap| against K / N^alpha(r,q): curves for different N should collapse (Theorem 1 scaling)."""130 d = s[(s.label == label) & (s.p == p)]131 if d.empty:132 return133 fig, axes = plt.subplots(1, len(ests), figsize=(2.4 * len(ests) + 0.6, 2.8), sharey=False)134 Ns = sorted(d.N.unique()); shades = ["#86b6ef", "#5598e7", "#2a78d6", "#1c5cab", "#0d366b"]135 for ax, e in zip(np.atleast_1d(axes), ests):136 for i, N in enumerate(Ns):137 dd = d[(d.estimator == e) & (d.N == N)].sort_values("K")138 if dd.empty:139 continue140 x = dd["K"] / N ** ALPHA[e]141 ax.plot(x, dd["rms_gap"].clip(lower=1e-12), marker="o", ms=3.5, color=shades[i % len(shades)], label=f"N = 2^{int(np.log2(N))}")142 ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel(f"K / N^{{{ALPHA[e]:.4g}}}"); ax.set_title(STYLE[e][3], fontsize=8, loc="left")143 np.atleast_1d(axes)[0].set_ylabel("RMS of √N (θ̃ − θ̂_N)"); np.atleast_1d(axes)[-1].legend(fontsize=7)144 fig.savefig(FIG / (fname or f"fig_gap_collapse_{label}_p{p}.pdf")); plt.close(fig)145146147def fig_coverage(s, label, p, N, ests=("r1_q0", "r1_q1", "r2_q0", "r2_q1", "r3_q1", "cs_q1"), level=95, fname=None):148 d = s[(s.label == label) & (s.p == p) & (s.N == N)]149 if d.empty or f"cov{level}" not in d.columns:150 return151 fig, ax = plt.subplots(figsize=(4.6, 3.0))152 for e in ests:153 dd = d[d.estimator == e].sort_values("m")154 if dd.empty or dd[f"cov{level}"].isna().all():155 continue156 col, mk, ls, lab = STYLE[e]157 ax.plot(dd["m"], dd[f"cov{level}"], marker=mk, ls=ls, color=PAL[col], label=lab, markerfacecolor="white", markeredgewidth=1.2)158 ax.axhline(level / 100, color=INK2, lw=0.8, ls="-"); ax.set_xscale("log", base=2); ax.set_xlabel("local sample size m"); ax.set_ylabel(f"empirical coverage, nominal {level}%")159 ax.set_title(f"{label}, p = {p}, N = {N:,}", loc="left", fontsize=9); ax.legend(); fig.savefig(FIG / (fname or f"fig_coverage{level}_{label}_p{p}_N{N}.pdf")); plt.close(fig)160161162def fig_boundary(fname="fig_boundary_exponential.pdf"):163 """Ratio of the Monte Carlo mean scaled gap to the asymptotic limit, against m, per (r,q) and c; hollow markers = finite-m prediction."""164 fp = PROC / "boundary_summary_with_prediction.csv"; f = PROC / "boundary_summary.csv"165 if not (fp.exists() or f.exists()):166 return167 b = pd.read_csv(fp if fp.exists() else f)168 cases = [(1, 1), (2, 0), (2, 1), (1, 2)]169 fig, axes = plt.subplots(1, len(cases), figsize=(2.45 * len(cases) + 0.6, 2.9), sharey=True)170 cols = ["blue", "orange", "aqua"]; mks = ["o", "s", "^"]171 for ax, (r, q) in zip(axes, cases):172 d = b[(b.r == r) & (b.q == q)]173 for i, c in enumerate(sorted(d.c.unique())):174 dd = d[d.c == c].sort_values("m")175 ratio = dd["mean_gap"] / dd["limit_c_eff"]; err = 1.96 * dd["se"] / dd["limit_c_eff"].abs()176 ax.errorbar(dd["m"], ratio, yerr=err, marker=mks[i], ms=4, color=PAL[cols[i]], label=f"c = {c:g}", capsize=2, lw=1.2, ls="-")177 if "pred_finite_m" in dd:178 ax.plot(dd["m"], dd["pred_finite_m"] / dd["limit_c_eff"], marker=mks[i], ms=6, markerfacecolor="white", color=PAL[cols[i]], lw=0, markeredgewidth=1.1)179 ax.axhline(1.0, color=INK2, lw=0.8); ax.set_xscale("log", base=2); ax.set_xlabel("local sample size m")180 ax.set_title(f"(r, q) = ({r}, {q}), α = {1-1/(2**(q+1)*r):.4g}", fontsize=8, loc="left")181 axes[0].set_ylabel("mean gap / asymptotic limit"); axes[0].legend(fontsize=7, loc="upper right")182 fig.savefig(FIG / fname); plt.close(fig)183184185def main():186 ap = argparse.ArgumentParser(); ap.add_argument("--figures", action="store_true"); a = ap.parse_args()187 _setup(); s = build_summary(); print("cells:", s.groupby(["label", "p", "N", "K"]).ngroups, "rows:", len(s))188 if a.figures:189 for (label, p), d in s.groupby(["label", "p"]):190 Nmax = d.N.max()191 ests = ("r1_q0", "r1_q1", "r1_q2", "r2_q0", "r2_q1", "r3_q0", "r3_q1", "cs_q1") + (() if label.startswith("exp") else ("csl1",))192 fig_rmse_vs_m(s, label, p, Nmax, "rms_gap", ests=ests, ylabel="RMS of √N (θ̃ − θ̂_N)")193 fig_rmse_vs_m(s, label, p, Nmax, "rmse", ests=("r1_q0", "r2_q0", "r3_q0", "r1_q1", "r2_q1", "cs_q0"), ylabel="RMSE of √N (θ̃ − θ0)", fname=f"fig_rmse_theta0_vs_m_{label}_p{p}_N{Nmax}.pdf")194 fig_gap_collapse(s, label, p); fig_coverage(s, label, p, Nmax)195 fig_boundary()196 print("figures:", len(list(FIG.glob("*.pdf"))))197198199if __name__ == "__main__":200 main()201