arikw/intervention-learning-framework

1"""Reproducible milestone-1 experiment: simulate -> detect -> estimate -> gate.2 3Runs the full milestone-1 loop on simulated data and prints:4  * simulator parameter recovery summary,5  * per-detector report (AUC, ECE, cost-optimal threshold, P/R),6  * per-moment effect table (pooled / stratified / IPW with 95% CIs, true ATE,7    confounding warnings),8  * a power/runnability gate for the headline effect.9 10Usage:11    python scripts/run_experiment.py --n-calls 25000 --seed 20260910 [--out results.json]12"""13 14from __future__ import annotations15 16import argparse17import json18from pathlib import Path19 20import numpy as np21from sklearn.metrics import roc_auc_score22 23from intervene.detect import (24    MomentDetector,25    expected_calibration_error,26    select_threshold_by_cost,27)28from intervene.effects import moment_effect_table29from intervene.planning import plan30from intervene.sim import ArmSpec, CallSimulator, DetectionSpec, MomentSpec, SimulatorConfig31 32FEATURES = [f"f{d}" for d in range(8)]33 34 35def build_config(seed: int, n_calls: int) -> SimulatorConfig:36    return SimulatorConfig(37        seed=seed,38        n_calls=n_calls,39        moments=(40            MomentSpec(41                id="pricing_objection",42                name="Pricing objection",43                base_logodds=-1.35,44                rep_coef=0.3,45                quality_coef=0.5,46                arms=(47                    ArmSpec(id="control", is_control=True),48                    ArmSpec(49                        id="discount_offer",50                        effect_logodds=0.0,51                        policy_intercept=-0.35,52                        policy_quality_coef=1.3,53                    ),54                ),55                detector=DetectionSpec(beta_pos=(9.0, 2.0), beta_neg=(1.5, 8.0)),56            ),57            MomentSpec(58                id="value_reframe",59                name="Value reframe",60                base_logodds=-1.9,61                rep_coef=0.2,62                quality_coef=0.3,63                arms=(64                    ArmSpec(id="control", is_control=True),65                    ArmSpec(id="reframe", effect_logodds=0.25),66                ),67                detector=DetectionSpec(beta_pos=(7.0, 2.5), beta_neg=(1.2, 7.0)),68            ),69        ),70    )71 72 73def main() -> None:74    ap = argparse.ArgumentParser(description=__doc__)75    ap.add_argument("--seed", type=int, default=20260910)76    ap.add_argument("--n-calls", type=int, default=25000)77    ap.add_argument("--out", type=Path, default=Path("results.json"))78    args = ap.parse_args()79 80    cfg = build_config(args.seed, args.n_calls)81    data = CallSimulator(cfg).generate()82    gt = data.ground_truth83 84    print("=" * 78)85    print(f"MILESTONE 1 EXPERIMENT  seed={args.seed}  n_calls={args.n_calls}")86    print(f"decisions logged: {data.meta['decided_calls']}  proxy-censored dropped: {data.n_proxy_censored}")87    print("=" * 78)88 89    # --- parameter recovery -------------------------------------------------90    print("\n[1] Simulator parameter recovery (empirical vs analytic quadrature)")91    for spec in cfg.moments:92        emp = data.calls_df["moment_id"].eq(spec.id).mean()93        true = gt.true_marginal_rate(spec.id)94        print(95            f"  {spec.id:22s} base rate: empirical {emp:.4f}  true {true:.4f}  "96            f"diff {emp - true:+.4f}"97        )98 99    # --- detection ----------------------------------------------------------100    print("\n[2] Detection (isotonic-calibrated logistic regression, 30/70 split, cost 3:1)")101    detector_report = {}102    for spec in cfg.moments:103        df = data.detector_frame[data.detector_frame["moment_id"] == spec.id].sample(104            frac=1.0, random_state=args.seed105        )106        n_train = int(0.3 * len(df))107        train, ev = df.iloc[:n_train], df.iloc[n_train:]108        det = MomentDetector(spec.id, calibration="isotonic", seed=args.seed).fit(109            train[FEATURES].to_numpy(), train["label"].to_numpy()110        )111        proba = det.predict_proba(ev[FEATURES].to_numpy())112        y = ev["label"].to_numpy()113        sel = select_threshold_by_cost(y, proba, c_fp=3.0, c_fn=1.0)114        k = int(np.argmin(np.abs(sel.curve.thresholds - sel.threshold)))115        detector_report[spec.id] = {116            "auc": float(roc_auc_score(y, proba)),117            "ece": expected_calibration_error(y, proba),118            "threshold": sel.threshold,119            "expected_cost": sel.expected_cost,120            "precision": float(sel.curve.precision[k]),121            "recall": float(sel.curve.recall[k]),122            "n_eval": int(len(y)),123        }124        r = detector_report[spec.id]125        print(126            f"  {spec.id:22s} AUC {r['auc']:.3f}  ECE {r['ece']:.4f}  "127            f"tau* {sel.threshold:.3f}  P {r['precision']:.3f}  R {r['recall']:.3f}  "128            f"cost {r['expected_cost']:.4f} (n_eval={r['n_eval']})"129        )130 131    # --- effects ------------------------------------------------------------132    print("\n[3] Offline effect estimation (95% Wald intervals; true ATE from quadrature)")133    all_rows = []134    for spec in cfg.moments:135        rows = moment_effect_table(data.decisions_df, spec.id, ground_truth=gt)136        all_rows.extend(rows)137        print(f"  moment: {spec.id}")138        for arm in sorted({r["arm"] for r in rows}):139            for est in ("pooled", "stratified", "ipw"):140                r = next(x for x in rows if x["arm"] == arm and x["estimator"] == est)141                star = " <-- covers true ATE" if r["estimate"].covers(r["true_ate"]) else ""142                print(143                    f"    arm {r['arm']:16s} {est:10s} "144                    f"{r['value']:+.4f} [{r['ci_low']:+.4f}, {r['ci_high']:+.4f}] "145                    f"(n={r['n']}, true ATE {r['true_ate']:+.4f}){star}"146                )147            warn = next(x for x in rows if x["estimator"] == "ipw")["confounding_warning"]148            if warn:149                print(f"    !! {warn}")150 151    # --- power gate ---------------------------------------------------------152    print("\n[4] Planning gate: detecting +3pt proxy lift at p1=0.24")153    gate = plan(154        p1=0.24,155        p2=0.27,156        eligible_rate=0.2,157        calls_per_day=500.0,158        objective_resolution_lag_days=5.0,159    )160    print(161        f"  n/arm {gate['n_per_arm']}  proxy {gate['proxy_days']:.1f} days  "162        f"objective {gate['objective_days']:.1f} days (+ objective resolution lag)"163    )164 165    out = {166        "seed": args.seed,167        "n_calls": args.n_calls,168        "decisions": data.meta["decided_calls"],169        "proxy_censored": data.n_proxy_censored,170        "recovery": {171            s.id: {172                "base_rate_empirical": float(data.calls_df["moment_id"].eq(s.id).mean()),173                "base_rate_true": gt.true_marginal_rate(s.id),174            }175            for s in cfg.moments176        },177        "detectors": detector_report,178        "effects": [{k: v for k, v in r.items() if k != "estimate"} for r in all_rows],179        "planning_gate": gate,180    }181    args.out.write_text(json.dumps(out, indent=2, default=str))182    print(f"\nwrote {args.out}")183 184 185if __name__ == "__main__":186    main()