13.3 从实验到决策

# file: sim/n13_decision.py
# tier: S
# verify: python3 sim/n13_decision.py --config config/ --out out/
"""从实验到决策（13.3）：把一次实验结果压成一个可执行的决策。

对一次实验的观测结果，算：
- 提升点估计、两比例 z 检验的 p 值、差值的 95% 置信区间（不合并方差）。
- 当前样本量下的最小可检测效应 MDE。
- 「让置信区间下界越过业务门槛」所需的每组样本量。
- 护栏指标是否越界。
- 决策判定：回滚 / 继续观察 / 全量。

产出 out/n13_decision.json + out/n13_decision.csv + out/n13_decision.svg。
"""

from __future__ import annotations

import argparse
import csv
import json
import math
import os
import sys

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from n13_stats import norm_cdf, norm_ppf, two_prop_z  # noqa: E402
from svg_line import render  # noqa: E402


def mde_permille(p: float, n: int, alpha: float, power: float) -> float:
    """等样本量下、当前 n 能检出的最小比例差（千分比）。"""
    za = norm_ppf(1.0 - alpha / 2.0)
    zb = norm_ppf(power)
    return (za + zb) * math.sqrt(2.0 * p * (1.0 - p) / n) * 1000.0


def n_for_ci_lower(diff: float, p: float, threshold: float, alpha: float) -> int:
    """让差值置信区间下界 >= threshold 所需的每组样本量（不合并方差近似）。"""
    z = norm_ppf(1.0 - alpha / 2.0)
    slack = diff - threshold
    if slack <= 0:
        return -1  # 点估计都没过门槛，加样本也没用
    return math.ceil(2.0 * p * (1.0 - p) * z * z / (slack * slack))


def main() -> int:
    here = os.path.dirname(os.path.abspath(__file__))
    ap = argparse.ArgumentParser()
    ap.add_argument("--config", default=os.path.join(here, "..", "config"))
    ap.add_argument("--out", default=os.path.join(here, "..", "out"))
    args = ap.parse_args()
    os.makedirs(args.out, exist_ok=True)

    with open(os.path.join(args.config, "n13_experiment.json"), encoding="utf-8") as fh:
        cfg = json.load(fh)
    ex = cfg["experiment"]
    dc = cfg["decision"]
    alpha = ex["alphaPermille"] / 1000.0
    power = ex["powerPermille"] / 1000.0

    o = dc["observed"]
    c1, n1, c2, n2 = o["controlConv"], o["controlN"], o["treatConv"], o["treatN"]
    p1, p2 = c1 / n1, c2 / n2
    diff = p2 - p1
    diff_permille = round(diff * 1000, 2)
    rel_lift_permille = round(diff / p1 * 1000, 1)

    z, pval = two_prop_z(c1, n1, c2, n2)
    z = float(z)
    pval = float(pval)

    # 差值的 95% CI（不合并方差）
    se = math.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    zc = norm_ppf(1.0 - alpha / 2.0)
    ci_lo = (diff - zc * se) * 1000.0
    ci_hi = (diff + zc * se) * 1000.0

    p_pool = (c1 + c2) / (n1 + n2)
    mde = mde_permille(p_pool, min(n1, n2), alpha, power)
    threshold = dc["minLiftPermille"] / 1000.0
    n_needed = n_for_ci_lower(diff, p_pool, threshold, alpha)

    # 护栏
    guardrails = []
    gr_ok = True
    for g in dc["guardrails"]:
        delta = g["treatPermille"] - g["controlPermille"]
        breach = abs(delta) > g["tolerancePermille"]
        gr_ok = gr_ok and not breach
        guardrails.append({"id": g["id"], "deltaPermille": delta,
                           "tolerancePermille": g["tolerancePermille"], "breach": breach})

    significant = pval < alpha
    if not significant:
        verdict = "hold"           # 不显著：继续观察 / 延长实验
    elif diff < 0:
        verdict = "rollback"       # 显著为负：回滚
    elif not gr_ok:
        verdict = "rollback"       # 护栏越界：回滚
    elif ci_lo / 1000.0 >= threshold:
        verdict = "ship"           # 显著为正、护栏正常、CI 下界过门槛：全量
    else:
        verdict = "hold"           # 显著为正，但还不能确信超过业务门槛

    # 样本量 vs MDE 曲线（用于 SVG 与 CSV）
    curve_n = [1000, 2000, 5000, 10000, 20000, 50000]
    curve = [{"nPerArm": n, "mdePermille": round(mde_permille(p_pool, n, alpha, power), 2)} for n in curve_n]

    with open(os.path.join(args.out, "n13_decision.csv"), "w", encoding="utf-8", newline="") as fh:
        w = csv.writer(fh)
        w.writerow(["nPerArm", "mdePermille"])
        for r in curve:
            w.writerow([r["nPerArm"], r["mdePermille"]])

    render([r["mdePermille"] for r in curve], os.path.join(args.out, "n13_decision.svg"),
           "MDE permille vs n per arm")

    report = {
        "version": 1,
        "controlRatePermille": round(p1 * 1000, 2), "treatRatePermille": round(p2 * 1000, 2),
        "diffPermille": diff_permille, "relativeLiftPermille": rel_lift_permille,
        "z": round(z, 4), "pValue": round(pval, 6), "significant": significant,
        "ciLowPermille": round(ci_lo, 2), "ciHighPermille": round(ci_hi, 2),
        "mdePermille": round(mde, 2), "minLiftPermille": dc["minLiftPermille"],
        "nPerArm": min(n1, n2), "nPerArmForThreshold": n_needed,
        "guardrails": guardrails, "guardrailsOk": gr_ok,
        "verdict": verdict, "curve": curve,
    }
    with open(os.path.join(args.out, "n13_decision.json"), "w", encoding="utf-8") as fh:
        json.dump(report, fh, ensure_ascii=False, indent=2)
    print(json.dumps(report, ensure_ascii=False))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())