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 if any(k in str(df["label"].iloc[0]) for k in ("misspec", "contam")):48 # pseudo-true value: the mean of the centralized pseudo-MLE over replications (its O(1/N) bias is negligible at the N used)49 th0 = np.broadcast_to(cen.mean(axis=0), th0.shape)50 rows = []51 B = len(df)52 for e in estimators(df):53 th = df[[f"{e}_{j}" for j in range(p)]].to_numpy()54 ok = np.all(np.isfinite(th), axis=1)55 gap = sN * (th - cen); err = sN * (th - th0)56 gap_n = np.linalg.norm(gap, axis=1); err_n = np.linalg.norm(err, axis=1)57 row = {"estimator": e, "reps": B, "finite_frac": ok.mean(),58 "mean_gap0": np.nanmean(gap[ok, 0]), "se_gap0": np.nanstd(gap[ok, 0], ddof=1) / np.sqrt(ok.sum()),59 "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),60 "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)),61 "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),62 "max_abs_bias": np.max(np.abs(np.nanmean(err[ok], axis=0)))}63 if f"se_{e}_0" in df.columns:64 se = df[[f"se_{e}_{j}" for j in range(p)]].to_numpy()65 for lv in (0.90, 0.95, 0.99):66 z = stats.norm.ppf(0.5 + lv / 2)67 hit = np.abs(th - th0) <= z * se68 row[f"cov{int(lv*100)}"] = np.nanmean(hit[ok]); row[f"len{int(lv*100)}"] = np.nanmean(2 * z * se[ok]) * sN69 if f"se_sand_central_0" in df.columns and e == "central":70 ses = df[[f"se_sand_central_{j}" for j in range(p)]].to_numpy()71 z = stats.norm.ppf(0.975); row["cov95_sandwich"] = np.nanmean((np.abs(th - th0) <= z * ses)[ok])72 rows.append(row)73 # centralized reference row74 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(),75 "sd0": err[:, 0].std(ddof=1), "mean_gap0": 0.0, "rms_gap": 0.0})76 if "se_central_0" in df.columns:77 se = df[[f"se_central_{j}" for j in range(p)]].to_numpy()78 for lv in (0.90, 0.95, 0.99):79 z = stats.norm.ppf(0.5 + lv / 2); rows[-1][f"cov{int(lv*100)}"] = np.mean(np.abs(cen - th0) <= z * se)80 if "se_sand_central_0" in df.columns:81 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)82 out = pd.DataFrame(rows)83 for c in ("label", "model", "p", "N", "K", "m"):84 out[c] = df[c].iloc[0]85 for r in (1, 2, 3):86 if f"valid_r{r}" in df.columns:87 out[f"valid_frac_r{r}"] = df[f"valid_r{r}"].mean() / df["K"].iloc[0]88 for s in (1, 2, 4):89 if f"local_fail_s{s}" in df.columns:90 out[f"local_fail_frac_s{s}"] = df[f"local_fail_s{s}"].mean() / (df["K"].iloc[0] * s)91 out["central_fail"] = (df["central_status"] != 0).mean()92 out["t_estimate_s"] = df["t_estimate_s"].mean()93 for c in df.columns:94 if c.startswith("cond_r2_q1"):95 out["cond_r2_q1"] = df[c].mean()96 return out979899def build_summary():100 frames = []101 for f in sorted(glob.glob(str(RAW / "*.csv"))):102 if "boundary" in f or "ablation" in f or "skew" in f:103 continue104 df = pd.read_csv(f)105 if len(df) < 10:106 continue107 frames.append(summarize_cell(df))108 s = pd.concat(frames, ignore_index=True)109 PROC.mkdir(parents=True, exist_ok=True); s.to_csv(PROC / "summary.csv", index=False)110 return s111112113def 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):114 d = s[(s.label == label) & (s.p == p) & (s.N == N)]115 if d.empty:116 return117 fig, ax = plt.subplots(figsize=(4.6, 3.4))118 for e in ests:119 dd = d[d.estimator == e].sort_values("m")120 if dd.empty or e not in STYLE:121 continue122 col, mk, ls, lab = STYLE[e]123 ax.plot(dd["m"], dd[metric].clip(lower=1e-12), marker=mk, ls=ls, color=PAL[col], label=lab, markerfacecolor="white", markeredgewidth=1.2)124 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)125 ax.set_ylim(bottom=max(ax.get_ylim()[0], 1e-8))126 ax.set_title(f"{label}, p = {p}, N = {N:,}", loc="left", fontsize=9, color=INK)127 ax.legend(ncol=2, loc="upper center", bbox_to_anchor=(0.5, -0.22), fontsize=7); FIG.mkdir(parents=True, exist_ok=True)128 out = FIG / (fname or f"fig_{metric}_vs_m_{label}_p{p}_N{N}.pdf"); fig.savefig(out)129 fig.savefig(FIG / out.name.replace(f"_N{N}", "")) # canonical copy (largest N) referenced by the paper130 plt.close(fig)131132133def fig_gap_collapse(s, label, p, ests=("r1_q1", "r2_q1", "r1_q2", "r3_q1"), fname=None):134 """|mean scaled gap| against K / N^alpha(r,q): curves for different N should collapse (Theorem 1 scaling)."""135 d = s[(s.label == label) & (s.p == p)]136 if d.empty:137 return138 fig, axes = plt.subplots(1, len(ests), figsize=(2.4 * len(ests) + 0.6, 2.8), sharey=False)139 Ns = sorted(d.N.unique()); shades = ["#86b6ef", "#5598e7", "#2a78d6", "#1c5cab", "#0d366b"]140 for ax, e in zip(np.atleast_1d(axes), ests):141 for i, N in enumerate(Ns):142 dd = d[(d.estimator == e) & (d.N == N)].sort_values("K")143 if dd.empty:144 continue145 x = dd["K"] / N ** ALPHA[e]146 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))}")147 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")148 np.atleast_1d(axes)[0].set_ylabel("RMS of √N (θ̃ − θ̂_N)"); np.atleast_1d(axes)[-1].legend(fontsize=7)149 fig.savefig(FIG / (fname or f"fig_gap_collapse_{label}_p{p}.pdf")); plt.close(fig)150151152def fig_coverage(s, label, p, N, ests=("r1_q0", "r1_q1", "r2_q0", "r2_q1", "r3_q1", "cs_q1"), level=95, fname=None):153 d = s[(s.label == label) & (s.p == p) & (s.N == N)]154 if d.empty or f"cov{level}" not in d.columns:155 return156 fig, ax = plt.subplots(figsize=(4.6, 3.0))157 for e in ests:158 dd = d[d.estimator == e].sort_values("m")159 if dd.empty or dd[f"cov{level}"].isna().all():160 continue161 col, mk, ls, lab = STYLE[e]162 ax.plot(dd["m"], dd[f"cov{level}"], marker=mk, ls=ls, color=PAL[col], label=lab, markerfacecolor="white", markeredgewidth=1.2)163 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}%")164 ax.set_title(f"{label}, p = {p}, N = {N:,}", loc="left", fontsize=9); ax.legend(fontsize=7, ncol=2, loc="lower right")165 out = FIG / (fname or f"fig_coverage{level}_{label}_p{p}_N{N}.pdf"); fig.savefig(out); fig.savefig(FIG / out.name.replace(f"_N{N}", "")); plt.close(fig)166167168def fig_boundary(fname="fig_boundary_exponential.pdf"):169 """Ratio of the Monte Carlo mean scaled gap to the asymptotic limit, against m, per (r,q) and c; hollow markers = finite-m prediction."""170 fp = PROC / "boundary_summary_with_prediction.csv"; f = PROC / "boundary_summary.csv"171 if not (fp.exists() or f.exists()):172 return173 b = pd.read_csv(fp if fp.exists() else f)174 cases = [(1, 1), (2, 0), (2, 1), (1, 2)]175 fig, axes = plt.subplots(1, len(cases), figsize=(2.45 * len(cases) + 0.6, 2.9), sharey=True)176 cols = ["blue", "orange", "aqua"]; mks = ["o", "s", "^"]177 for ax, (r, q) in zip(axes, cases):178 d = b[(b.r == r) & (b.q == q)]179 for i, c in enumerate(sorted(d.c.unique())):180 dd = d[d.c == c].sort_values("m")181 ratio = dd["mean_gap"] / dd["limit_c_eff"]; err = 1.96 * dd["se"] / dd["limit_c_eff"].abs()182 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="-")183 if "pred_finite_m" in dd:184 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)185 ax.axhline(1.0, color=INK2, lw=0.8); ax.set_xscale("log", base=2); ax.set_xlabel("local sample size m")186 ax.set_title(f"(r, q) = ({r}, {q}), α = {1-1/(2**(q+1)*r):.4g}", fontsize=8, loc="left")187 axes[0].set_ylabel("mean gap / asymptotic limit"); axes[0].legend(fontsize=7, loc="upper right")188 fig.savefig(FIG / fname); plt.close(fig)189190191def main():192 ap = argparse.ArgumentParser(); ap.add_argument("--figures", action="store_true"); a = ap.parse_args()193 _setup(); s = build_summary(); print("cells:", s.groupby(["label", "p", "N", "K"]).ngroups, "rows:", len(s))194 if a.figures:195 for (label, p), d in s.groupby(["label", "p"]):196 Nmax = d.N.max()197 ests = ("r1_q0", "r1_q1", "r1_q2", "r2_q0", "r2_q1", "r3_q0", "r3_q1", "cs_q1") + (() if label.startswith("exp") else ("csl1",))198 fig_rmse_vs_m(s, label, p, Nmax, "rms_gap", ests=ests, ylabel="RMS of √N (θ̃ − θ̂_N)")199 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")200 fig_gap_collapse(s, label, p); fig_coverage(s, label, p, Nmax)201 fig_boundary()202 print("figures:", len(list(FIG.glob("*.pdf"))))203204205if __name__ == "__main__":206 main()207