"""Shared-baseline coupling: h = 100 - m and Delta = t - m share the estimated control score m, so control-score noise e_m enters both, h_obs = h - e_m and Delta_obs = Delta - e_m + e_t, and induces a positive covariance between h and Delta that a fit treating h as fixed reads as slope. Frontier report. Local re-tabulation only: no model API, VM, evaluation or raw log. Run (from the paper directory): python3 -I evidence/coupling.py Writes coupling.json and coupling.txt beside this file. Inputs (read only) LEDGER forensics/ledger.csv in the paper directory, the 195-row comparison record: est, se, ctrl, ours, mx, headroom_pts, n In the released ledger.csv these free-text cells are empty, and every numeric column and every categorical column except team is as produced: note on L009 L036 L039 L042 L043 L052 L073 L079 L083 L086 L091 L099 L100 L101 L102 L104 L105 L107 L108 L109 L110 L112 L115 L121 L122 L125 L129 L130 L138 L139 L140 L148 L150 L165 L166 L168 L176 L178 L193; settings_parity on L036 L037 L039 L052 L086 L102 L108 L110 L139 L140 L142 L143 L148 L193; instrument on L169 L170 L171 L173; contrast on L169 L170 L171 L172 L173. Agent and team names are replaced by neutral IDs or plain study names in the team column and in the instrument, contrast, settings_parity and note cells of L030 L034 L038 L044 L045 L052 L082 L103 L107-L110 L133 L134 L151 L152 L184-L191 L193. Every source cell says in plain words where the comparison was recorded (internal document names, paths and commit ids removed), as do the note of L068 and the instrument of L146. No number this script prints reads those cells. N_EFF tasks and repeats per row where the ledger lacks n, from the paper's protocol notes (below) HFJ headroom_forecast.json (this folder): the 21 rows and the headline values, asserted equal Model code (mean, nll, fit, STARTS) as in the synthesis paper's re-test of the law, used to validate the fast fitter. Sections 1 control-score SE per row, se_m = sqrt(m (100 - m) / n_eff), n_eff = tasks x repeats; mean se_m, Var(h) and the expected spurious slope Var(e_m) / Var(h_obs) for the 15 Opus 5.5 and the 6 weaker rows 2 errors-in-variables correction of k (method of moments for the weighted through-origin fit, and SIMEX on the paper's own likelihood fit), and the corrected weak-to-Opus forecast (Bayes factors, own tau; MAE) 3 parametric bootstrap of coupling-only nulls (true lift constant c per population, one common constant, or 0; control noise in both h and Delta; with and without the null model's fitted tau), plus the same nulls without coupling (se_m = 0) to separate coupling from the rest of the null; the weak-fitted forecast scored on each simulated Opus set; null percentiles against the observed values Coupling: the worst case consistent with the row's own SE s. If se_m <= s, Delta noise = -e_m + u with Var(u) = s^2 - se_m^2 (arms independent: Cov(h noise, Delta noise) = se_m^2). If se_m > s (a paired design whose arms are positively correlated), Delta noise = -(s / se_m) e_m: total Delta variance s^2, correlation 1 with the h noise, covariance s se_m. In paired designs the true coupling is smaller than this. """ import io import json import math import multiprocessing as mp import os import sys import time for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"): os.environ.setdefault(_v, "1") import numpy as np import pandas as pd from scipy import optimize, stats HERE = os.path.dirname(os.path.abspath(__file__)) OUT = os.path.join(HERE, "coupling.json") TXT = os.path.join(HERE, "coupling.txt") HFJ = os.path.join(HERE, "headroom_forecast.json") LEDGER = os.path.join(os.path.dirname(HERE), "forensics", "ledger.csv") SEED = 20261030 B_NULL = 4000 # simulated datasets per null configuration B_SIMEX = 500 # SIMEX replicates per lambda RES = {} class _Tee: def __init__(self, *streams): self.streams = streams def write(self, s): for st in self.streams: st.write(s) def flush(self): for st in self.streams: st.flush() _BUF = io.StringIO() sys.stdout = _Tee(sys.__stdout__, _BUF) T0 = time.time() OPUS = ["L003", "L004", "L010", "L011", "L012", "L016", "L057", "L099", "L124", "L126", "L136", "L147", "L163", "L166", "L174"] WEAK = ["L001", "L002", "L006", "L007", "L054", "L058"] # n_eff = tasks x repeats behind the control score m (the paper's protocol notes) N_EFF = { "L003": (87, 1, "87 eligible tasks, one branch per arm"), "L004": (87, 1, "87 eligible tasks, one branch per arm"), "L010": (30, 1, "30 probe items per hour; m = 100, so se_m = 0 for any n"), "L011": (16, 1, "16 Terminal-Bench tasks per probe hour (ledger n)"), "L012": (15, 2, "15 tau2 tasks x 2 trials = 30 items per probe hour"), "L016": (60, 1, "60 tasks, k = 1 per cell"), "L057": (66, 1, "66 tasks per pair"), "L099": (120, 1, "n = 120 = B u B', each episode once"), "L124": (112, 3, "112 post-drift tasks x 3 repeats"), "L126": (173, 3, "173 pooled post-drift tasks x 3 repeats"), "L136": (224, 3, "224 pooled tasks x 3 repeats"), "L147": (61, 3, "61 post-drift tasks x 3 repeats"), "L163": (417, 1, "417 tasks, one run per arm"), "L166": (212, 1, "212 tasks, one run per arm"), "L174": (417, 1, "417 tasks, one run per arm; the row's control 84.9 = 354/417 is the whole mixed-model stream, " "so n is the whole stream (on the Opus 5.5 half both arms are 195/209)"), "L001": (417, 2, "417 test_challenge tasks; control = mean of 2 bare runs"), "L002": (168, 2, "168 test_normal tasks x 2 bare runs"), "L006": (65, 1, "65 eligible tasks"), "L007": (58, 1, "58 eligible tasks"), "L054": (66, 1, "66 tasks per pair"), "L058": (84, 1, "84 tasks per arm, the second half (the ledger's 168 is the whole stream)"), } # ------------------------------------------------------------------------------------------------------------ # rows L = pd.read_csv(LEDGER, keep_default_na=False) for c in ["est", "lo", "hi", "p", "n", "se", "z", "ours", "ctrl", "mx", "headroom_pts"]: L[c] = pd.to_numeric(L[c], errors="coerce") def rows(ids, scale_n="tasks_x_repeats"): d = L.set_index("id").loc[ids].reset_index() d["y"] = d.est * 100 / d.mx d["s"] = d.se * 100 / d.mx d["h"] = d.headroom_pts d["m"] = d.ctrl * 100 / d.mx d["t"] = d.ours * 100 / d.mx d["tasks"] = [N_EFF[i][0] for i in d.id] d["repeats"] = [N_EFF[i][1] for i in d.id] d["n_eff"] = d.tasks * (d.repeats if scale_n == "tasks_x_repeats" else 1) d["se_m"] = np.sqrt(np.clip(d.m * (100 - d.m), 0, None) / d.n_eff) d["cpl"] = np.minimum(d.se_m ** 2, d.s * d.se_m) # worst-case Cov(h noise, Delta noise) return d A = rows(OPUS) W = rows(WEAK) A_cons = rows(OPUS, "tasks") W_cons = rows(WEAK, "tasks") hf = json.load(open(HFJ)) _hf_rows = {r["id"]: r for r in hf["D5_rows"]["opus15"] + hf["D5_rows"]["weaker6"]} for d in (A, W): for _, r in d.iterrows(): q = _hf_rows[r.id] assert abs(q["delta"] - round(r.y, 4)) < 1e-9 and abs(q["s"] - round(r.s, 4)) < 1e-9 and \ abs(q["h"] - round(r.h, 4)) < 1e-9, r.id assert abs(q["control_score"] - round(r.m, 4)) < 1e-9, r.id assert abs((100 - r.m) - r.h) < 0.051, (r.id, r.m, r.h) # h = 100 - m up to the ledger's rounding # ------------------------------------------------------------------------------------------------------------ # model code, verbatim from the synthesis paper's re-test code (used to validate the fast fitter) def mean(kind, th, h): if kind == "M0": return np.zeros_like(h) if kind == "M1": return np.full_like(h, th[0]) if kind == "M2": return th[0] * h if kind == "M2b": return th[0] * h - th[1] if kind == "M3": return np.minimum(th[0], h) if kind == "M5": c = math.exp(th[0]) return c * (1 - np.exp(-h / c)) raise ValueError(kind) NPAR = {"M0": 0, "M1": 1, "M2": 1, "M2b": 2, "M3": 1, "M5": 1} STARTS = {"M0": [[]], "M1": [[0.0], [5.0]], "M2": [[0.3], [0.6]], "M2b": [[0.3, 0.0], [0.6, 1.0]], "M3": [[3.0], [8.0], [0.0]], "M5": [[math.log(3)], [math.log(10)], [math.log(30)]]} def nll(p, d, kind, tau_fixed=None): th, tau = (p, tau_fixed) if tau_fixed is not None else (p[:-1], math.exp(p[-1])) sd = np.sqrt(d.s.values ** 2 + tau ** 2) return -stats.norm.logpdf(d.y.values, mean(kind, th, d.h.values), sd).sum() def fit(d, kind, tau_fixed=None): best = None for st in STARTS[kind]: for lt in ([0.0, 1.5, 2.5] if tau_fixed is None else [None]): x0 = list(st) + ([lt] if lt is not None else []) if not x0: continue r = optimize.minimize(nll, x0, args=(d, kind, tau_fixed), method="Nelder-Mead", options=dict(maxiter=40000, xatol=1e-9, fatol=1e-11)) if best is None or r.fun < best.fun: best = r if best is None: # M0 with fixed tau: no free parameter return [], tau_fixed, nll([], d, kind, tau_fixed) th = list(best.x[:-1]) if tau_fixed is None else list(best.x) tau = math.exp(best.x[-1]) if tau_fixed is None else tau_fixed return th, tau, float(best.fun) # ------------------------------------------------------------------------------------------------------------ # fast fitter (same likelihood, same MLE; exact inner solutions, profiled log tau) for the simulations LOG2PI = math.log(2 * math.pi) LT_GRID = np.linspace(-7.0, 4.5, 24) def _nll(y, mu, var): return 0.5 * float(np.sum(LOG2PI + np.log(var) + (y - mu) ** 2 / var)) def _inner(kind, y, h, var): w = 1.0 / var if kind == "M0": return _nll(y, 0.0, var), [] if kind == "M1": c = float(np.sum(w * y) / np.sum(w)) return _nll(y, c, var), [c] if kind == "M2": k = float(np.sum(w * h * y) / np.sum(w * h * h)) return _nll(y, k * h, var), [k] if kind == "M3": # exact: c in each interval between sorted h, rows with h > c are capped at c hs = np.sort(h) edges = np.concatenate([[-np.inf], hs, [np.inf]]) best = (np.inf, None) for j in range(len(edges) - 1): lo, hi = edges[j], edges[j + 1] S = h > lo if np.isfinite(lo) else np.ones_like(h, dtype=bool) if not S.any(): c = hs[-1] else: c = float(np.sum(w[S] * y[S]) / np.sum(w[S])) c = min(max(c, lo), hi) if np.isfinite(hi) else max(c, lo) f = _nll(y, np.minimum(c, h), var) if f < best[0] - 1e-12: best = (f, [c]) return best if kind == "M5": f = lambda lc: _nll(y, math.exp(lc) * (1 - np.exp(-h / math.exp(lc))), var) g = np.linspace(math.log(0.05), math.log(5000), 61) v = np.array([f(x) for x in g]) i = int(np.argmin(v)) lo, hi = g[max(i - 1, 0)], g[min(i + 1, len(g) - 1)] r = optimize.minimize_scalar(f, bounds=(lo, hi), method="bounded", options=dict(xatol=1e-7)) return (float(r.fun), [float(r.x)]) if r.fun <= v[i] else (float(v[i]), [float(g[i])]) raise ValueError(kind) def fit_fast(y, s, h, kind, k_fixed=None): s2 = s ** 2 if k_fixed is not None: prof = lambda lt: (_nll(y, k_fixed * h, s2 + math.exp(2 * lt)), [k_fixed]) else: prof = lambda lt: _inner(kind, y, h, s2 + math.exp(2 * lt)) v = [prof(x)[0] for x in LT_GRID] i = int(np.argmin(v)) lo, hi = LT_GRID[max(i - 1, 0)], LT_GRID[min(i + 1, len(LT_GRID) - 1)] if i == 0: lo = -30.0 r = optimize.minimize_scalar(lambda x: prof(x)[0], bounds=(lo, hi), method="bounded", options=dict(xatol=1e-7)) lt = float(r.x) if r.fun <= v[i] else float(LT_GRID[i]) f, th = prof(lt) return th, math.exp(lt), f def lpd_sum(kind, th, tau, y, s, h): return float(stats.norm.logpdf(y, mean(kind, th, h), np.sqrt(s ** 2 + tau ** 2)).sum()) KINDS = ["M0", "M1", "M2", "M3", "M5"] print("=" * 100) print("0. FAST FITTER VALIDATION against the verbatim Nelder-Mead fit (nll, parameter, tau)") val = [] for lab, d in (("opus15", A), ("weaker6", W), ("all21", pd.concat([A, W]))): for kind in KINDS: th0, tau0, f0 = fit(d, kind) th1, tau1, f1 = fit_fast(d.y.values, d.s.values, d.h.values, kind) dpar = max([abs(a - b) for a, b in zip(th0, th1)] or [0.0]) val.append(dict(rows=lab, model=kind, nll_nm=round(f0, 5), nll_fast=round(f1, 5), d_nll=round(f1 - f0, 6), param_nm=[round(x, 4) for x in th0], param_fast=[round(x, 4) for x in th1], tau_nm=round(tau0, 4), tau_fast=round(tau1, 4))) print(f" {lab:8s} {kind}: nll NM {f0:9.4f} fast {f1:9.4f} (d {f1 - f0:+.1e}) params {np.round(th0, 4)} " f"vs {np.round(th1, 4)} tau {tau0:.3f} vs {tau1:.3f}") assert max(v["d_nll"] for v in val) < 1e-3, "fast fitter worse than Nelder-Mead" RES["fast_fitter_validation"] = val # observed headline values recomputed with the fast fitter def forecast_stats(yw, sw, hw, yo, so, ho): fw = {k: fit_fast(yw, sw, hw, k) for k in KINDS} lp = {k: lpd_sum(k, fw[k][0], fw[k][1], yo, so, ho) for k in KINDS} mae = {k: float(np.mean(np.abs(yo - mean(k, fw[k][0], ho)))) for k in KINDS} return fw, lp, mae fw_obs, lp_obs, mae_obs = forecast_stats(W.y.values, W.s.values, W.h.values, A.y.values, A.s.values, A.h.values) OBS = dict(logbf_M2_M1=lp_obs["M2"] - lp_obs["M1"], logbf_M2_M3=lp_obs["M2"] - lp_obs["M3"], logbf_M2_M0=lp_obs["M2"] - lp_obs["M0"], logbf_M2_M5=lp_obs["M2"] - lp_obs["M5"], mae_gap_M1_M2=mae_obs["M1"] - mae_obs["M2"], mae_gap_M3_M2=mae_obs["M3"] - mae_obs["M2"], mae_gap_M0_M2=mae_obs["M0"] - mae_obs["M2"], k_weak=fw_obs["M2"][0][0]) _fo1, _fo2 = fit_fast(A.y.values, A.s.values, A.h.values, "M1"), fit_fast(A.y.values, A.s.values, A.h.values, "M2") OBS["aic_gap_opus_M1_M2"] = 2 * (_fo1[2] - _fo2[2]) OBS["k_opus"] = _fo2[0][0] print("\nObserved (fast fitter): BF law:constant %.2f, law:capped %.2f, law:zero %.2f, law:smooth cap %.2f; " "MAE gap constant-law %.3f; Opus-15 AIC gap constant-law %.2f; k weak %.4f, k Opus %.4f" % ( math.exp(OBS["logbf_M2_M1"]), math.exp(OBS["logbf_M2_M3"]), math.exp(OBS["logbf_M2_M0"]), math.exp(OBS["logbf_M2_M5"]), OBS["mae_gap_M1_M2"], OBS["aic_gap_opus_M1_M2"], OBS["k_weak"], OBS["k_opus"])) hb = hf["headline"] assert abs(math.exp(OBS["logbf_M2_M0"]) - hb["forecast_bf_own_tau"]["vs_zero"]) < 0.05 assert abs(math.exp(OBS["logbf_M2_M1"]) - hb["forecast_bf_own_tau"]["vs_constant"]) < 0.01 assert abs(math.exp(OBS["logbf_M2_M3"]) - hb["forecast_bf_own_tau"]["vs_capped"]) < 0.01 RES["observed"] = {k: round(v, 4) for k, v in OBS.items()} RES["observed"].update(bf_M2_M1=round(math.exp(OBS["logbf_M2_M1"]), 2), bf_M2_M3=round(math.exp(OBS["logbf_M2_M3"]), 2), bf_M2_M0=round(math.exp(OBS["logbf_M2_M0"]), 2), bf_M2_M5=round(math.exp(OBS["logbf_M2_M5"]), 2)) # ------------------------------------------------------------------------------------------------------------ print("\n" + "=" * 100) print("1. CONTROL-SCORE STANDARD ERROR PER ROW: se_m = sqrt(m (100 - m) / n_eff)") def coupling_summary(d, lab): tau2 = fit_fast(d.y.values, d.s.values, d.h.values, "M2")[1] ** 2 w = 1 / (d.s.values ** 2 + tau2) var_h = float(np.var(d.h.values, ddof=1)) out = dict(rows=lab, n=int(len(d)), mean_se_m=float(d.se_m.mean()), rms_se_m=float(np.sqrt((d.se_m ** 2).mean())), var_e_m=float((d.se_m ** 2).mean()), var_h_obs=var_h, sd_h_obs=math.sqrt(var_h), spurious_slope=float((d.se_m ** 2).mean() / var_h), spurious_slope_worstcase_cov=float(d.cpl.mean() / var_h), spurious_k_through_origin_weighted=float(np.sum(w * d.cpl) / np.sum(w * d.h.values ** 2)), rows_with_se_m_above_s=d.id[d.se_m > d.s].tolist()) return out S1 = {} for lab, d, dc in (("opus15", A, A_cons), ("weaker6", W, W_cons)): print(f"\n {lab}: id, m, tasks x repeats = n_eff, se_m, s (row SE), h, worst-case Cov(h noise, Delta noise)") for (_, r), (_, rc) in zip(d.iterrows(), dc.iterrows()): print(f" {r.id} m {r.m:6.2f} {r.tasks:4d} x {r.repeats} = {r.n_eff:5d} se_m {r.se_m:5.2f} " f"(tasks only {rc.se_m:5.2f}) s {r.s:5.2f} h {r.h:5.1f} cov {r.cpl:6.2f}") S1[lab] = dict(primary=coupling_summary(d, lab), tasks_only=coupling_summary(dc, lab), per_row=[dict(id=r.id, m=round(r.m, 3), tasks=int(r.tasks), repeats=int(r.repeats), n_eff=int(r.n_eff), se_m=round(r.se_m, 4), se_m_tasks_only=round(rc.se_m, 4), s=round(r.s, 4), h=round(r.h, 3), cov_worstcase=round(r.cpl, 4), n_source=N_EFF[r.id][2]) for (_, r), (_, rc) in zip(d.iterrows(), dc.iterrows())]) for v in ("primary", "tasks_only"): o = S1[lab][v] print(f" {lab} [{v}]: mean se_m {o['mean_se_m']:.3f} (rms {o['rms_se_m']:.3f}); Var(e_m) {o['var_e_m']:.3f}; " f"Var(h_obs) {o['var_h_obs']:.2f} (sd {o['sd_h_obs']:.2f}); expected spurious slope Var(e_m)/Var(h_obs) " f"{o['spurious_slope']:.4f} (worst-case cov {o['spurious_slope_worstcase_cov']:.4f}); through-origin " f"weighted analogue {o['spurious_k_through_origin_weighted']:.4f}; se_m > s on {o['rows_with_se_m_above_s']}") RES["S1_control_se"] = S1 # ------------------------------------------------------------------------------------------------------------ print("\n" + "=" * 100) print("2. ERRORS-IN-VARIABLES CORRECTION OF k") print(" model: h_obs = h - e_m, Delta_obs = k h - e_m + e_t + eta; Var(e_m) = se_m^2, Cov(h noise, Delta noise) = cov") def mom_k(d): """Method of moments for the weighted through-origin fit: E[h_obs D_obs] = k h^2 + cov, E[h_obs^2] = h^2 + se^2.""" y, s, h, se, cv = d.y.values, d.s.values, d.h.values, d.se_m.values, d.cpl.values th, tau, _ = fit_fast(y, s, h, "M2") k = th[0] for _ in range(200): w = 1 / (s ** 2 + tau ** 2) k_new = float(np.sum(w * (h * y - cv)) / np.sum(w * (h * h - se * se))) tau = fit_fast(y, s, h, "M2", k_fixed=k_new)[1] if abs(k_new - k) < 1e-10: k = k_new break k = k_new return k, tau def simex_k(d, rng, lambdas=(0.5, 1.0, 1.5, 2.0)): y, s, h, se, cv = d.y.values, d.s.values, d.h.values, d.se_m.values, d.cpl.values n = len(d) a = np.where(se > 0, cv / np.where(se > 0, se, 1), 0.0) # Delta-noise loading on the h-noise draw resid_sd = np.sqrt(np.clip(s ** 2 - a ** 2, 0, None)) pts = {0.0: fit_fast(y, s, h, "M2")[0][0]} for lam in lambdas: ks = [] for _ in range(B_SIMEX): z, v = rng.standard_normal(n), rng.standard_normal(n) hs = h - math.sqrt(lam) * se * z ys = y - math.sqrt(lam) * a * z + math.sqrt(lam) * resid_sd * v ks.append(fit_fast(ys, s * math.sqrt(1 + lam), hs, "M2")[0][0]) pts[lam] = float(np.mean(ks)) lam = np.array(sorted(pts)) est = np.array([pts[x] for x in lam]) quad = np.polyfit(lam, est, 2) lin = np.polyfit(lam, est, 1) return float(np.polyval(quad, -1.0)), float(np.polyval(lin, -1.0)), {str(k): round(v, 5) for k, v in pts.items()} rng = np.random.default_rng(SEED) S2 = {} for lab, d, dc in (("opus15", A, A_cons), ("weaker6", W, W_cons)): S2[lab] = {} for v, dd in (("primary", d), ("tasks_only", dc)): k_naive = fit_fast(dd.y.values, dd.s.values, dd.h.values, "M2")[0][0] k_mom, tau_mom = mom_k(dd) k_sq, k_sl, pts = simex_k(dd, rng) S2[lab][v] = dict(k_naive=round(k_naive, 4), k_mom=round(k_mom, 4), tau_at_k_mom=round(tau_mom, 3), k_simex_quadratic=round(k_sq, 4), k_simex_linear=round(k_sl, 4), simex_points=pts) print(f" {lab:8s} [{v:10s}] k naive {k_naive:.4f} -> method of moments {k_mom:.4f} (tau {tau_mom:.2f}); " f"SIMEX quadratic {k_sq:.4f}, linear {k_sl:.4f}; SIMEX means by lambda {pts}") # corrected forecast: weak k replaced by its corrected value, weak tau re-profiled at that k, Opus scored (own tau) print("\n Corrected weak-to-Opus forecast (law with corrected weak k; tau re-profiled on the weak rows at that k;" " other models as fitted)") FC = {} for v in ("primary", "tasks_only"): for meth in ("k_mom", "k_simex_quadratic"): kc = S2["weaker6"][v][meth] tau_c = fit_fast(W.y.values, W.s.values, W.h.values, "M2", k_fixed=kc)[1] lp_c = lpd_sum("M2", [kc], tau_c, A.y.values, A.s.values, A.h.values) mae_c = float(np.mean(np.abs(A.y.values - kc * A.h.values))) rmse_c = float(np.sqrt(np.mean((A.y.values - kc * A.h.values) ** 2))) o = dict(k=kc, tau_weak=round(tau_c, 3), bf_vs_zero=round(math.exp(lp_c - lp_obs["M0"]), 2), bf_vs_constant=round(math.exp(lp_c - lp_obs["M1"]), 2), bf_vs_capped=round(math.exp(lp_c - lp_obs["M3"]), 2), bf_vs_smooth_cap=round(math.exp(lp_c - lp_obs["M5"]), 2), mae=round(mae_c, 3), rmse=round(rmse_c, 3)) FC[f"{v}/{meth}"] = o print(f" [{v:10s} {meth:17s}] k {kc:.4f} tau_w {tau_c:.2f}: BF vs zero {o['bf_vs_zero']}, constant " f"{o['bf_vs_constant']}, capped {o['bf_vs_capped']}, smooth cap {o['bf_vs_smooth_cap']}; " f"MAE {o['mae']} RMSE {o['rmse']}") print(f" [uncorrected ] k {OBS['k_weak']:.4f}: BF vs zero {RES['observed']['bf_M2_M0']}, constant " f"{RES['observed']['bf_M2_M1']}, capped {RES['observed']['bf_M2_M3']}; MAE {mae_obs['M2']:.3f}") S2["corrected_forecast"] = FC RES["S2_errors_in_variables"] = S2 # ------------------------------------------------------------------------------------------------------------ print("\n" + "=" * 100) print(f"3. PARAMETRIC BOOTSTRAP OF COUPLING-ONLY NULLS ({B_NULL} datasets each)") print(" each population simulated with true h = observed h, control noise e_m ~ N(0, se_m^2) in both h and Delta") print(" (h kept >= 0), remaining Delta noise so Var(Delta noise) = s^2 (worst-case coupling, see header).") print(" true lift: 'constant' = the constant c fitted on that population (weak 8.12, Opus 5.71: the two populations") print(" differ); 'common' = one constant for both (the all-21 constant fit); 'zero' = 0. 'tau0' = no between-row") print(" spread (the primary specification); 'taufit' = the null model's own fitted tau (standard parametric") print(" bootstrap). 'nocoupling' = the same null with se_m = 0, to separate coupling from the rest of the null.") print(" Bayes factors and MAE: models fitted on the simulated weak set, scored on the simulated Opus set (own tau).") def simulate(rg, d, lift, tau_het, coupled=True): h_t, s = d.h.values, d.s.values se = d.se_m.values if coupled else np.zeros(len(d)) n = len(d) e = np.minimum(rg.standard_normal(n) * se, h_t) h_sim = h_t - e ind = se <= s u = rg.standard_normal(n) * np.sqrt(np.clip(s ** 2 - se ** 2, 0, None)) dn = np.where(ind, -e + u, -(s / np.where(se > 0, se, 1.0)) * e) return h_sim, lift + rg.standard_normal(n) * tau_het + dn fitW = {k: fit_fast(W.y.values, W.s.values, W.h.values, k) for k in ("M0", "M1")} fitA = {k: fit_fast(A.y.values, A.s.values, A.h.values, k) for k in ("M0", "M1")} fit21 = fit_fast(np.concatenate([W.y.values, A.y.values]), np.concatenate([W.s.values, A.s.values]), np.concatenate([W.h.values, A.h.values]), "M1") # name: (true lift, between-row spread, n_eff variant, coupled) CONFIGS = { "constant_tau0": ("constant", False, "primary", True), "constant_taufit": ("constant", True, "primary", True), "common_tau0": ("common", False, "primary", True), "common_taufit": ("common", True, "primary", True), "zero_tau0": ("zero", False, "primary", True), "zero_taufit": ("zero", True, "primary", True), "constant_tau0_nocoupling": ("constant", False, "primary", False), "constant_taufit_nocoupling": ("constant", True, "primary", False), "common_tau0_nocoupling": ("common", False, "primary", False), "common_taufit_nocoupling": ("common", True, "primary", False), "zero_tau0_nocoupling": ("zero", False, "primary", False), "zero_taufit_nocoupling": ("zero", True, "primary", False), "constant_taufit_tasks_only": ("constant", True, "tasks_only", True), "common_taufit_tasks_only": ("common", True, "tasks_only", True), "zero_taufit_tasks_only": ("zero", True, "tasks_only", True), } METRICS = [("lbf_M1", "logbf_M2_M1", "BF law:constant"), ("lbf_M3", "logbf_M2_M3", "BF law:capped"), ("lbf_M0", "logbf_M2_M0", "BF law:zero"), ("lbf_M5", "logbf_M2_M5", "BF law:smooth cap")] GAPS = [("gap_M1", "mae_gap_M1_M2", "MAE constant - law"), ("gap_M3", "mae_gap_M3_M2", "MAE capped - law"), ("gap_M0", "mae_gap_M0_M2", "MAE zero - law")] def qtiles(x): return {f"p{q}": round(float(np.percentile(x, q)), 4) for q in (50, 95, 99)} def run_config(name): null, het, nv, coupled = CONFIGS[name] dW, dA = (W, A) if nv == "primary" else (W_cons, A_cons) kind0 = "M0" if null == "zero" else "M1" if null == "constant": liftW, liftA = fitW["M1"][0][0], fitA["M1"][0][0] elif null == "common": liftW = liftA = fit21[0][0] else: liftW = liftA = 0.0 tauW = fitW[kind0][1] if het else 0.0 tauA = fitA[kind0][1] if het else 0.0 rg = np.random.default_rng([SEED, sum(map(ord, name)), len(name)]) rec = {k: [] for k in ("lbf_M1", "lbf_M3", "lbf_M0", "lbf_M5", "gap_M1", "gap_M3", "gap_M0", "k_weak", "k_opus", "aic_gap_opus")} for _ in range(B_NULL): hw, yw = simulate(rg, dW, liftW, tauW, coupled) ho, yo = simulate(rg, dA, liftA, tauA, coupled) fw, lp, mae = forecast_stats(yw, dW.s.values, hw, yo, dA.s.values, ho) for o in ("M1", "M3", "M0", "M5"): rec[f"lbf_{o}"].append(lp["M2"] - lp[o]) for o in ("M1", "M3", "M0"): rec[f"gap_{o}"].append(mae[o] - mae["M2"]) rec["k_weak"].append(fw["M2"][0][0]) f1, f2 = fit_fast(yo, dA.s.values, ho, "M1"), fit_fast(yo, dA.s.values, ho, "M2") rec["aic_gap_opus"].append(2 * (f1[2] - f2[2])) rec["k_opus"].append(f2[0][0]) R = {k: np.array(v) for k, v in rec.items()} out = dict(null=null, coupled=coupled, true_lift_weak=round(liftW, 4), true_lift_opus=round(liftA, 4), tau_het_weak=round(tauW, 3), tau_het_opus=round(tauA, 3), n_eff=nv, datasets=B_NULL) for key, obs_key, lab in METRICS: lbf = R[key] out[lab] = dict(observed=round(math.exp(OBS[obs_key]), 2), observed_log10=round(OBS[obs_key] / math.log(10), 3), null_log10={q: round(v / math.log(10), 3) for q, v in qtiles(lbf).items()}, frac_null_ge_observed=round(float(np.mean(lbf >= OBS[obs_key])), 5), count_ge_observed=int(np.sum(lbf >= OBS[obs_key])), frac_null_bf_gt_1=round(float(np.mean(lbf > 0)), 4)) for key, obs_key, lab in GAPS: out[lab] = dict(observed=round(OBS[obs_key], 3), null=qtiles(R[key]), frac_null_ge_observed=round(float(np.mean(R[key] >= OBS[obs_key])), 5), count_ge_observed=int(np.sum(R[key] >= OBS[obs_key]))) g = R["aic_gap_opus"] out["Opus-15 AIC constant - law"] = dict(observed=round(OBS["aic_gap_opus_M1_M2"], 2), null=qtiles(g), frac_null_ge_observed=round(float(np.mean(g >= OBS["aic_gap_opus_M1_M2"])), 5), count_ge_observed=int(np.sum(g >= OBS["aic_gap_opus_M1_M2"]))) out["k weak (fitted on null data)"] = dict(observed=round(OBS["k_weak"], 4), null_mean=round(float(R["k_weak"].mean()), 4), null=qtiles(R["k_weak"])) out["k Opus-15 (fitted on null data)"] = dict(observed=round(OBS["k_opus"], 4), null_mean=round(float(R["k_opus"].mean()), 4), null=qtiles(R["k_opus"]), count_ge_observed=int(np.sum(R["k_opus"] >= OBS["k_opus"]))) return name, out t1 = time.time() with mp.get_context("fork").Pool(processes=min(len(CONFIGS), max(1, (os.cpu_count() or 2) // 4))) as pool: S3 = dict(pool.map(run_config, list(CONFIGS))) print(f" ({len(CONFIGS)} configurations in {time.time() - t1:.0f} s)") for name in CONFIGS: out = S3[name] print(f"\n [{name}] true lift weak {out['true_lift_weak']:.3f} / Opus {out['true_lift_opus']:.3f}; spread tau " f"weak {out['tau_het_weak']:.2f} / Opus {out['tau_het_opus']:.2f}; n_eff {out['n_eff']}; coupled " f"{out['coupled']}") for _, _, lab in METRICS: o = out[lab] print(f" {lab:22s} observed {o['observed']:>8} (log10 {o['observed_log10']:6.2f}) null log10 median " f"{o['null_log10']['p50']:7.2f} 95th {o['null_log10']['p95']:7.2f} 99th {o['null_log10']['p99']:7.2f} " f"null >= observed {o['count_ge_observed']}/{B_NULL} (null BF > 1: {o['frac_null_bf_gt_1']:.2f})") for _, _, lab in GAPS + [(None, None, "Opus-15 AIC constant - law")]: o = out[lab] print(f" {lab:26s} observed {o['observed']:>7} null median {o['null']['p50']:8.3f} 95th " f"{o['null']['p95']:8.3f} 99th {o['null']['p99']:8.3f} null >= observed {o['count_ge_observed']}/{B_NULL}") for lab in ("k weak (fitted on null data)", "k Opus-15 (fitted on null data)"): o = out[lab] print(f" {lab:32s} observed {o['observed']} null mean {o['null_mean']} 95th {o['null']['p95']} " f"99th {o['null']['p99']}" + (f" null >= observed {o['count_ge_observed']}/{B_NULL}" if "count_ge_observed" in o else "")) RES["S3_bootstrap_nulls"] = S3 print("\n" + "=" * 100) print("4. WHAT COUPLING ADDS: each coupled null against the same null with se_m = 0 (share of null datasets at or") print(" beyond the observed value; k = mean Opus-15 k fitted on null data)") SUMM = [] for base in ("constant_tau0", "constant_taufit", "common_tau0", "common_taufit", "zero_tau0", "zero_taufit"): c, u = S3[base], S3[base + "_nocoupling"] row = dict(null=base) for lab in ("BF law:constant", "BF law:capped", "BF law:zero", "MAE constant - law", "Opus-15 AIC constant - law"): row[lab] = dict(coupled=c[lab]["frac_null_ge_observed"], uncoupled=u[lab]["frac_null_ge_observed"]) row["k Opus-15 null mean"] = dict(coupled=c["k Opus-15 (fitted on null data)"]["null_mean"], uncoupled=u["k Opus-15 (fitted on null data)"]["null_mean"]) row["k weak null mean"] = dict(coupled=c["k weak (fitted on null data)"]["null_mean"], uncoupled=u["k weak (fitted on null data)"]["null_mean"]) SUMM.append(row) print(f" {base:16s} " + "; ".join(f"{k} {v['coupled']:.4f} vs {v['uncoupled']:.4f}" for k, v in row.items() if k != "null")) RES["S4_coupling_contribution"] = SUMM RES["meta"] = dict(seed=SEED, B_null=B_NULL, B_simex=B_SIMEX, n_eff=N_EFF, runtime_s=round(time.time() - T0, 1), coupling_rule="worst case: Cov(h noise, Delta noise) = min(se_m^2, s se_m), Var(Delta noise) = s^2") with open(OUT, "w") as fh: json.dump(RES, fh, indent=1, default=float) print("\nwrote", os.path.basename(OUT), f"({time.time() - T0:.0f} s)") sys.stdout = sys.__stdout__ with open(TXT, "w") as fh: fh.write(_BUF.getvalue()) print("wrote", os.path.basename(TXT))