#!/usr/bin/env python3 """Aggregate raw Monte Carlo CSVs into results/processed/summary.csv and draw the paper figures (PDF, vector). Every number in the paper's Monte Carlo tables/figures is produced here from results/raw/*.csv (no manual entry). Usage: python3 experiments/process_results.py [--figures] """ import argparse, glob, json, pathlib, re, sys import numpy as np, pandas as pd import matplotlib; matplotlib.use("Agg") import matplotlib.pyplot as plt from scipy import stats ROOT = pathlib.Path(__file__).resolve().parents[1] RAW, PROC, FIG = ROOT / "results" / "raw", ROOT / "results" / "processed", ROOT / "results" / "figures" PAL = {"blue": "#2a78d6", "orange": "#eb6834", "aqua": "#1baf7a", "yellow": "#eda100", "magenta": "#e87ba4", "green": "#008300", "violet": "#6250d6", "red": "#e34948"} INK, INK2, GRID = "#0b0b0b", "#52514e", "#e6e5e1" STYLE = { # fixed assignment estimator -> (color, marker, linestyle, label) "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)"), "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)"), "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)"), "cs_q0": ("violet", "D", "-", "analytic Cox–Snell (q=0)"), "cs_q1": ("violet", "D", "--", "analytic Cox–Snell + 1 Newton"), "csl1": ("green", "v", "-", "CSL, 1 round"), "csl2": ("green", "v", "--", "CSL, 2 rounds"), "aee": ("yellow", "P", "-", "AEE (Lin–Xi)"), "stale_r1_q1": ("magenta", "x", "--", "raw + 1 stale-Hessian step"), "stale_r2_q1": ("magenta", "x", ":", "half-block + 1 stale-Hessian step"), } ALPHA = {"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} def _setup(): plt.rcParams.update({"font.family": "serif", "font.size": 9, "axes.edgecolor": INK2, "axes.labelcolor": INK, "xtick.color": INK2, "ytick.color": INK2, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6, "axes.spines.top": False, "axes.spines.right": False, "legend.frameon": False, "legend.fontsize": 8, "figure.dpi": 120, "savefig.bbox": "tight", "lines.linewidth": 1.6, "lines.markersize": 4.5, "axes.facecolor": "#fcfcfb", "figure.facecolor": "white"}) def estimators(df): names = set() for c in df.columns: m = re.match(r"^(r\d_q\d|stale_r\d_q\d|cs_q\d|csl\d|aee)_0$", c) if m: names.add(m.group(1)) return sorted(names) def summarize_cell(df): p = int(df["p"].iloc[0]); N = int(df["N"].iloc[0]); sN = np.sqrt(N) th0 = df[[f"theta0_{j}" for j in range(p)]].to_numpy() cen = df[[f"central_{j}" for j in range(p)]].to_numpy() rows = [] B = len(df) for e in estimators(df): th = df[[f"{e}_{j}" for j in range(p)]].to_numpy() ok = np.all(np.isfinite(th), axis=1) gap = sN * (th - cen); err = sN * (th - th0) gap_n = np.linalg.norm(gap, axis=1); err_n = np.linalg.norm(err, axis=1) row = {"estimator": e, "reps": B, "finite_frac": ok.mean(), "mean_gap0": np.nanmean(gap[ok, 0]), "se_gap0": np.nanstd(gap[ok, 0], ddof=1) / np.sqrt(ok.sum()), "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), "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)), "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), "max_abs_bias": np.max(np.abs(np.nanmean(err[ok], axis=0)))} if f"se_{e}_0" in df.columns: se = df[[f"se_{e}_{j}" for j in range(p)]].to_numpy() for lv in (0.90, 0.95, 0.99): z = stats.norm.ppf(0.5 + lv / 2) hit = np.abs(th - th0) <= z * se row[f"cov{int(lv*100)}"] = np.nanmean(hit[ok]); row[f"len{int(lv*100)}"] = np.nanmean(2 * z * se[ok]) * sN if f"se_sand_central_0" in df.columns and e == "central": ses = df[[f"se_sand_central_{j}" for j in range(p)]].to_numpy() z = stats.norm.ppf(0.975); row["cov95_sandwich"] = np.nanmean((np.abs(th - th0) <= z * ses)[ok]) rows.append(row) # centralized reference row 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(), "sd0": err[:, 0].std(ddof=1), "mean_gap0": 0.0, "rms_gap": 0.0}) if "se_central_0" in df.columns: se = df[[f"se_central_{j}" for j in range(p)]].to_numpy() for lv in (0.90, 0.95, 0.99): z = stats.norm.ppf(0.5 + lv / 2); rows[-1][f"cov{int(lv*100)}"] = np.mean(np.abs(cen - th0) <= z * se) if "se_sand_central_0" in df.columns: 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) out = pd.DataFrame(rows) for c in ("label", "model", "p", "N", "K", "m"): out[c] = df[c].iloc[0] for r in (1, 2, 3): if f"valid_r{r}" in df.columns: out[f"valid_frac_r{r}"] = df[f"valid_r{r}"].mean() / df["K"].iloc[0] for s in (1, 2, 4): if f"local_fail_s{s}" in df.columns: out[f"local_fail_frac_s{s}"] = df[f"local_fail_s{s}"].mean() / (df["K"].iloc[0] * s) out["central_fail"] = (df["central_status"] != 0).mean() out["t_estimate_s"] = df["t_estimate_s"].mean() for c in df.columns: if c.startswith("cond_r2_q1"): out["cond_r2_q1"] = df[c].mean() return out def build_summary(): frames = [] for f in sorted(glob.glob(str(RAW / "*.csv"))): if "boundary" in f or "ablation" in f or "skew" in f: continue df = pd.read_csv(f) if len(df) < 10: continue frames.append(summarize_cell(df)) s = pd.concat(frames, ignore_index=True) PROC.mkdir(parents=True, exist_ok=True); s.to_csv(PROC / "summary.csv", index=False) return s def 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): d = s[(s.label == label) & (s.p == p) & (s.N == N)] if d.empty: return fig, ax = plt.subplots(figsize=(4.6, 3.4)) for e in ests: dd = d[d.estimator == e].sort_values("m") if dd.empty or e not in STYLE: continue col, mk, ls, lab = STYLE[e] ax.plot(dd["m"], dd[metric].clip(lower=1e-12), marker=mk, ls=ls, color=PAL[col], label=lab, markerfacecolor="white", markeredgewidth=1.2) 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) ax.set_title(f"{label}, p = {p}, N = {N:,}", loc="left", fontsize=9, color=INK) ax.legend(ncol=1, loc="best"); FIG.mkdir(parents=True, exist_ok=True) fig.savefig(FIG / (fname or f"fig_{metric}_vs_m_{label}_p{p}_N{N}.pdf")); plt.close(fig) def fig_gap_collapse(s, label, p, ests=("r1_q1", "r2_q1", "r1_q2", "r3_q1"), fname=None): """|mean scaled gap| against K / N^alpha(r,q): curves for different N should collapse (Theorem 1 scaling).""" d = s[(s.label == label) & (s.p == p)] if d.empty: return fig, axes = plt.subplots(1, len(ests), figsize=(2.4 * len(ests) + 0.6, 2.8), sharey=False) Ns = sorted(d.N.unique()); shades = ["#86b6ef", "#5598e7", "#2a78d6", "#1c5cab", "#0d366b"] for ax, e in zip(np.atleast_1d(axes), ests): for i, N in enumerate(Ns): dd = d[(d.estimator == e) & (d.N == N)].sort_values("K") if dd.empty: continue x = dd["K"] / N ** ALPHA[e] 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))}") 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") np.atleast_1d(axes)[0].set_ylabel("RMS of √N (θ̃ − θ̂_N)"); np.atleast_1d(axes)[-1].legend(fontsize=7) fig.savefig(FIG / (fname or f"fig_gap_collapse_{label}_p{p}.pdf")); plt.close(fig) def fig_coverage(s, label, p, N, ests=("r1_q0", "r1_q1", "r2_q0", "r2_q1", "r3_q1", "cs_q1"), level=95, fname=None): d = s[(s.label == label) & (s.p == p) & (s.N == N)] if d.empty or f"cov{level}" not in d.columns: return fig, ax = plt.subplots(figsize=(4.6, 3.0)) for e in ests: dd = d[d.estimator == e].sort_values("m") if dd.empty or dd[f"cov{level}"].isna().all(): continue col, mk, ls, lab = STYLE[e] ax.plot(dd["m"], dd[f"cov{level}"], marker=mk, ls=ls, color=PAL[col], label=lab, markerfacecolor="white", markeredgewidth=1.2) 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}%") 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) def fig_boundary(fname="fig_boundary_exponential.pdf"): f = PROC / "boundary_summary.csv" if not f.exists(): return b = pd.read_csv(f) cases = [(2, 1), (1, 2), (3, 1), (2, 0)] fig, axes = plt.subplots(1, len(cases), figsize=(2.5 * len(cases) + 0.5, 2.9)) cols = ["blue", "orange", "aqua"] for ax, (r, q) in zip(axes, cases): d = b[(b.r == r) & (b.q == q)] for i, c in enumerate(sorted(d.c.unique())): dd = d[d.c == c].sort_values("N") ax.errorbar(dd["N"], dd["mean_gap"], yerr=1.96 * dd["se"], marker="o", ms=3.5, color=PAL[cols[i]], label=f"c = {c:g}", capsize=2, lw=1.2) ax.plot(dd["N"], dd["limit_c_eff"], ls="--", color=PAL[cols[i]], lw=1.0) ax.set_xscale("log", base=2); ax.set_xlabel("N"); ax.set_title(f"(r, q) = ({r}, {q}), α = {1-1/(2**(q+1)*r):.4g}", fontsize=8, loc="left") axes[0].set_ylabel("mean of √N (θ̃ − θ̂_N)"); axes[0].legend(fontsize=7) fig.savefig(FIG / fname); plt.close(fig) def main(): ap = argparse.ArgumentParser(); ap.add_argument("--figures", action="store_true"); a = ap.parse_args() _setup(); s = build_summary(); print("cells:", s.groupby(["label", "p", "N", "K"]).ngroups, "rows:", len(s)) if a.figures: for (label, p), d in s.groupby(["label", "p"]): Nmax = d.N.max() fig_rmse_vs_m(s, label, p, Nmax, "rms_gap", ylabel="RMS of √N (θ̃ − θ̂_N)") 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") fig_gap_collapse(s, label, p); fig_coverage(s, label, p, Nmax) fig_boundary() print("figures:", len(list(FIG.glob("*.pdf")))) if __name__ == "__main__": main()