score_calibration.py

1"""2多业务分数校准：IR / SIR / Bayesian-SIR（第 8 章 §8.2.1 分数校准方法）34对应正文概念5------------6正文指出混排前处理必须做业务分校准，否则不同业务模型的高估/低估会直接导致7流量错配。给出的三级方案为：8  1. IR（Isotonic Regression，保序回归，fielding1974statistical）：分桶统计真实后验9     r_k = Σ_{i∈B_k} y_i / |B_k|，用 PAVA（de2010isotone）合并违反单调的相邻桶；10  2. SIR（Smoothed Isotonic Regression，deng2020calibrating）：在 IR 的分段常数映射上11     做线性插值，消除桶边界跳变；12  3. Bayesian-SIR：对小样本桶做 Beta-Binomial 共轭贝叶斯平滑，抑制高方差。1314关键公式15--------16    单调约束:     f(p_i) <= f(p_j),  ∀ p_i <= p_j17    桶后验:       r_k = Σ_{i∈B_k} y_i / |B_k|18    PAVA 合并:    r_new = (|B_k| r_k + |B_{k+1}| r_{k+1}) / (|B_k| + |B_{k+1}|)19    SIR 插值:     f(p) = r_k + (p - p_k)/(p_{k+1} - p_k) · (r_{k+1} - r_k)20    贝叶斯平滑:   CTR|Data ~ Beta(α + Click, β + PV - Click)21                  E[CTR] = (Click + α) / (PV + α + β)2223本实现做了什么24--------------25在合成的「模型系统性高估」日志上依次拟合 IR / SIR / Bayesian-SIR，26打印校准前后的 COPC（= Σpred / Σlabel，正文称 PCOC）与 ECE（期望校准误差），27并复现正文表 8.1 的小样本分桶案例（PV=50, Click=20, 经验 CTR=0.40）说明28为什么需要贝叶斯平滑。2930注：正文称 SIR 由「阿里妈妈团队在 2018 年提出」，而引用文献 deng2020calibrating31为 ECML-PKDD 2020 的 "Calibrating User Response Predictions in Online Advertising"。32本实现以正文给出的公式为准。33"""3435import math36import random37import statistics383940# ------------------------------------------------------------------ 指标41def copc(preds: list[float], labels: list[float]) -> float:42    """COPC / PCOC = Σ预估 / Σ真实，等于 1 表示无偏。"""43    s = sum(labels)44    return float("inf") if s == 0 else sum(preds) / s454647def ece(preds: list[float], labels: list[float], n_bin: int = 20) -> float:48    """Expected Calibration Error = Σ_b (n_b/N)·|acc_b - conf_b|。"""49    bins: list[list[tuple[float, float]]] = [[] for _ in range(n_bin)]50    for p, y in zip(preds, labels):51        bins[min(n_bin - 1, int(p * n_bin))].append((p, y))52    n = len(preds)53    return sum(len(b) / n * abs(sum(y for _, y in b) / len(b) - sum(p for p, _ in b) / len(b))54               for b in bins if b)555657# ------------------------------------------------------------------ 校准器58class IsotonicCalibrator:59    """保序回归：等频分桶 + PAVA，映射为分段常数函数 f(p) = c_k。6061    prior_strength > 0 时启用分层贝叶斯平滑（Bayesian-SIR）：每个细桶向62    「同区间的粗桶经验 CTR」收缩，即 Beta 先验的均值取局部估计而非全局常数，63    这样既能缩减稀疏桶的方差，又不会给整条校准曲线引入全局偏置。64        r_k = (Click_k + κ·m_k) / (PV_k + κ),  m_k = 所属粗桶的经验 CTR65    """6667    def __init__(self, n_bucket: int = 20, prior_strength: float = 0.0,68                 coarse_ratio: int = 5):69        self.n_bucket = n_bucket70        self.prior_strength = prior_strength71        self.coarse_ratio = coarse_ratio72        self.blocks: list[dict] = []7374    def fit(self, preds: list[float], labels: list[float]) -> "IsotonicCalibrator":75        pairs = sorted(zip(preds, labels))76        size = max(1, len(pairs) // self.n_bucket)77        raw: list[dict] = []78        for s in range(0, len(pairs), size):79            chunk = pairs[s:s + size]80            if not chunk:81                continue82            pv, click = len(chunk), sum(y for _, y in chunk)83            raw.append({"lo": chunk[0][0], "hi": chunk[-1][0],84                        "mid": sum(p for p, _ in chunk) / pv,85                        "pv": pv, "click": click, "emp": click / pv, "r": click / pv})86        # 分层先验：把相邻 coarse_ratio 个细桶聚成一个粗桶，粗桶经验 CTR 作为先验均值87        if self.prior_strength > 0:88            kappa, cr = self.prior_strength, max(1, self.coarse_ratio)89            for s in range(0, len(raw), cr):90                grp = raw[s:s + cr]91                pv = sum(b["pv"] for b in grp)92                m = sum(b["click"] for b in grp) / pv93                for b in grp:94                    b["r"] = (b["click"] + kappa * m) / (b["pv"] + kappa)95        self.n_raw_block = len(raw)96        merged: list[dict] = []97        for b in raw:98            merged.append(b)99            while len(merged) >= 2 and merged[-2]["r"] > merged[-1]["r"]:100                b2, b1 = merged.pop(), merged.pop()101                pv, click = b1["pv"] + b2["pv"], b1["click"] + b2["click"]102                merged.append({103                    "lo": b1["lo"], "hi": b2["hi"],104                    "mid": (b1["mid"] * b1["pv"] + b2["mid"] * b2["pv"]) / pv,105                    "pv": pv, "click": click, "emp": click / pv,106                    "r": (b1["pv"] * b1["r"] + b2["pv"] * b2["r"]) / pv,107                })108        self.blocks = merged109        return self110111    def predict(self, p: float) -> float:112        for b in self.blocks:113            if p <= b["hi"]:114                return b["r"]115        return self.blocks[-1]["r"]116117118class SmoothedIsotonicCalibrator(IsotonicCalibrator):119    """SIR：在 IR 的单调锚点 (p_k, r_k) 之间做线性插值，得到连续分段线性校准函数。"""120121    def predict(self, p: float) -> float:122        bs = self.blocks123        if p <= bs[0]["mid"]:124            return bs[0]["r"]125        if p >= bs[-1]["mid"]:126            return bs[-1]["r"]127        for k in range(len(bs) - 1):128            if bs[k]["mid"] <= p <= bs[k + 1]["mid"]:129                w = (p - bs[k]["mid"]) / max(1e-12, bs[k + 1]["mid"] - bs[k]["mid"])130                return bs[k]["r"] + w * (bs[k + 1]["r"] - bs[k]["r"])131        return bs[-1]["r"]132133134# ------------------------------------------------------------------ 数据135def make_logs(n: int, bias: float):136    """合成日志：真实后验 t ~ Beta(2,5)，模型预估 p = clip(bias·t + 噪声)。"""137    preds, labels, truths = [], [], []138    for _ in range(n):139        t = random.betavariate(2.0, 5.0)140        p = max(1e-4, min(0.9999, bias * t + random.gauss(0.0, 0.03)))141        preds.append(p)142        truths.append(t)143        labels.append(1.0 if random.random() < t else 0.0)144    return preds, labels, truths145146147def main() -> None:148    random.seed(42)149    print("=" * 76)150    print("多业务分数校准：IR / SIR / Bayesian-SIR（§8.2.1）")151    print("=" * 76)152153    # ---- 1) 大样本：IR / SIR 校准效果154    tr_p, tr_y, _ = make_logs(40000, bias=1.45)      # 训练集：广告 CTR 模型高估 45%155    te_p, te_y, _ = make_logs(20000, bias=1.45)156    print(f"\n【1】大样本场景：训练 {len(tr_p)} / 测试 {len(te_p)}，模型高估系数 bias=1.45")157158    ir = IsotonicCalibrator(n_bucket=20).fit(tr_p, tr_y)159    sir = SmoothedIsotonicCalibrator(n_bucket=20).fit(tr_p, tr_y)160161    rows = [("校准前(原始pxtr)", te_p),162            ("IR  保序回归", [ir.predict(p) for p in te_p]),163            ("SIR 平滑保序回归", [sir.predict(p) for p in te_p])]164165    print(f"\n{'方法':<20}{'COPC':>10}{'ECE':>10}{'|COPC-1|':>12}")166    print("-" * 52)167    for name, pr in rows:168        c = copc(pr, te_y)169        print(f"{name:<20}{c:>10.4f}{ece(pr, te_y):>10.5f}{abs(c - 1):>12.4f}")170171    print(f"\nIR 分桶：等频 {ir.n_raw_block} 桶 → PAVA 合并后 {len(ir.blocks)} 桶（单调）")172    print(f"{'区间':<22}{'PV':>7}{'Click':>7}{'经验CTR':>10}{'IR值':>10}")173    for b in ir.blocks[:6]:174        print(f"[{b['lo']:.4f},{b['hi']:.4f}]{'':<4}{b['pv']:>7d}{int(b['click']):>7d}"175              f"{b['emp']:>10.4f}{b['r']:>10.4f}")176177    # ---- 2) 桶边界跳变：IR 阶梯 vs SIR 连续178    print("\n【2】桶边界跳变对比（相邻两个预估分几乎相同，但落在不同桶）：")179    edge = ir.blocks[len(ir.blocks) // 2]["hi"]180    print(f"{'p':>10}{'IR f(p)':>12}{'SIR f(p)':>12}")181    for d in (-0.004, -0.001, 0.001, 0.004):182        p = edge + d183        print(f"{p:>10.4f}{ir.predict(p):>12.4f}{sir.predict(p):>12.4f}")184    jump_ir = abs(ir.predict(edge + 0.001) - ir.predict(edge - 0.001))185    jump_sir = abs(sir.predict(edge + 0.001) - sir.predict(edge - 0.001))186    print(f"边界 ±0.001 跳变幅度：IR={jump_ir:.5f}  SIR={jump_sir:.5f}")187188    # ---- 3) 小样本：Bayesian-SIR 的方差缩减价值189    #     正文强调贝叶斯平滑是针对「分桶样本稀疏 → 后验统计高方差」的方差缩减技术，190    #     因此只有在小样本分桶下它才优于纯 SIR；大样本时先验会被数据淹没。191    print("\n【3】小样本场景：SIR vs Bayesian-SIR（2500 条日志切 50 桶，每桶约 50 条）")192    print("     Bayesian-SIR 采用分层先验：细桶向所属粗桶（5 个细桶一组）的经验 CTR 收缩，κ=60")193    grid = [i / 40 for i in range(1, 40)]194    n_trial = 15195    acc = {"SIR": [0.0, 0.0], "Bayesian-SIR": [0.0, 0.0]}196    curves: dict[str, list[list[float]]] = {"SIR": [], "Bayesian-SIR": []}197    for _ in range(n_trial):198        s_p, s_y, _ = make_logs(2500, bias=1.45)      # 稀疏日志199        q_p, q_y, _ = make_logs(20000, bias=1.45)200        for name, cal in (("SIR", SmoothedIsotonicCalibrator(50)),201                          ("Bayesian-SIR", SmoothedIsotonicCalibrator(50, 60.0, 5))):202            cal.fit(s_p, s_y)203            out = [cal.predict(p) for p in q_p]204            acc[name][0] += abs(copc(out, q_y) - 1.0)205            acc[name][1] += ece(out, q_y)206            curves[name].append([cal.predict(g) for g in grid])207    print(f"  {'方法':<16}{'平均|COPC-1|':>14}{'平均ECE':>12}{'曲线方差':>12}"208          f"  （{n_trial} 次重复实验）")209    for name, (a, b) in acc.items():210        cs = curves[name]211        var = statistics.fmean([statistics.pvariance([c[j] for c in cs])212                                for j in range(len(grid))])213        print(f"  {name:<16}{a / n_trial:>14.4f}{b / n_trial:>12.5f}{var:>12.6f}")214    print("  「曲线方差」= 不同随机日志拟合出的校准函数在同一 p 处的方差均值，越小越稳。")215216    # ---- 4) 正文表 8.1 小样本案例：为什么要贝叶斯平滑217    print("\n【4】正文表 8.1 小样本分桶案例（Beta 先验 α=30, β=170，先验均值 0.15）：")218    alpha, beta = 30.0, 170.0219    print(f"{'pCTR分桶':<14}{'PV':>8}{'Click':>8}{'经验CTR':>10}{'平滑后CTR':>12}{'收缩量':>10}")220    for lo, hi, pv, ck in [(0.1, 0.2, 100000, 12000), (0.2, 0.3, 80000, 18000),221                           (0.3, 0.4, 50, 20)]:222        emp = ck / pv223        sm = (ck + alpha) / (pv + alpha + beta)224        print(f"[{lo},{hi}]{'':<6}{pv:>8d}{ck:>8d}{emp:>10.4f}{sm:>12.4f}{emp - sm:>+10.4f}")225    print("可见大样本桶几乎不变，小样本桶 0.4000 被拉回 0.2000，噪声被有效抑制。")226227    # ---- 5) 多业务：校准让不同业务的分数进入同一可信空间228    print("\n【5】多业务校准效果（各业务 bias 不同，SIR 校准后 COPC 全部回到 1 附近）：")229    print(f"{'业务':<10}{'bias':>7}{'COPC前':>10}{'COPC后':>10}{'ECE前':>10}{'ECE后':>10}")230    for biz, bias in [("广告", 1.45), ("直播", 0.80), ("电商", 1.15), ("短视频", 1.02)]:231        p_tr, y_tr, _ = make_logs(30000, bias)232        p_te, y_te, _ = make_logs(15000, bias)233        c = SmoothedIsotonicCalibrator(20).fit(p_tr, y_tr)234        out = [c.predict(p) for p in p_te]235        print(f"{biz:<10}{bias:>7.2f}{copc(p_te, y_te):>10.4f}{copc(out, y_te):>10.4f}"236              f"{ece(p_te, y_te):>10.5f}{ece(out, y_te):>10.5f}")237238239if __name__ == "__main__":240    main()