d2co.py

1"""2D²Co：两成分 GMM 估 w⁺/w⁻，再用校正函数还原兴趣（第23章 §23.4.3）34正文的分解：5    w = p_x^r · w_d⁺ + (1 − p_x^r) · w_d⁻6其中 w_d⁺ 是「感兴趣时的平均观看时长」（含 Duration Bias），7w_d⁻ 是「不感兴趣时的噪声观看」（几秒的误触/滑走，Noisy Watching）。89流程：10    1) 按 Duration 分组，在每组内用两成分 GMM 拟合观看时长分布，11       得到 ŵ_d⁺ 与 ŵ_d⁻；12    2) 用频次加权的双向移动平均抹平相邻 Duration 上的估计噪声；13    3) 线性校正 D²Co(A)：r = (w − w_d⁻) / (w_d⁺ − w_d⁻)；14       指数校正 D²Co(S)：r = (e^{αw} − e^{αw_d⁻}) / (e^{αw_d⁺} − e^{αw_d⁻})。1516本例的特殊之处：合成数据里保留了真实兴趣 r，所以可以拿 GMM 估出来的17w⁺/w⁻ 去和「按 r 分组的真实均值」对比——这是线上做不到的，只在这里18用来判断估计误差有多大、又会怎样传到最终 Label 上。1920工程取舍：GMM 拟合在 log(1+W) 上做，而不是原始秒数。原始 W 右偏且被人为21封顶在 [0, d]，直接套高斯混合会让两个成分都去追那条长尾。log 空间里22「长视频看满」与「划走」才是两个可分的团。回变换用 E[W]=exp(μ+σ²/2)−1，23带方差修正，否则会系统性低估均值。2425运行：python3 ch04_classic_methods/d2co.py26"""2728from __future__ import annotations2930import pathlib31import sys3233import numpy as np3435sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent.parent))3637from common import (  # noqa: E40238    auc, mae, make_watch_data, split, table,39)4041SEED = 2026424344# --------------------------------------------------------------------------45# 一维两成分 GMM（EM，多次随机重启，取对数似然最大的解）46# --------------------------------------------------------------------------47def _log_norm_pdf(x: np.ndarray, mu: float, sd: float) -> np.ndarray:48    return -0.5 * np.log(2 * np.pi) - np.log(sd) - 0.5 * ((x - mu) / sd) ** 2495051def fit_gmm1d(x: np.ndarray, k: int = 2, iters: int = 300,52              restarts: int = 4, seed: int = SEED):53    """54    返回按均值升序排列的 (mu, sd, pi)。5556    第 0 次用分位点初始化（稳定），其余用随机样本初始化（避免所有重启57    都落在同一个局部最优）。EM 的对数似然单调不减，所以「多重启取最大」58    是安全的。59    """60    rng = np.random.default_rng(seed)61    best = None62    for r in range(restarts):63        if r == 0:64            mu = np.sort(np.quantile(x, np.linspace(0.15, 0.85, k)))65        else:66            mu = np.sort(rng.choice(x, size=k, replace=False))67            # 随机重启可能抽到几乎重合的两点，加一点抖动保证可分68            mu = mu + rng.normal(0, 1e-3, size=k)69        sd = np.full(k, max(x.std(), 1e-3))70        pi = np.full(k, 1.0 / k)71        ll_old = -np.inf72        for _ in range(iters):73            logp = np.stack([_log_norm_pdf(x, mu[j], sd[j]) for j in range(k)], axis=1)74            logp = logp + np.log(pi + 1e-300)75            m = logp.max(axis=1, keepdims=True)76            lse = m[:, 0] + np.log(np.exp(logp - m).sum(axis=1))77            resp = np.exp(logp - lse[:, None])78            nk = resp.sum(axis=0) + 1e-30079            mu = (resp * x[:, None]).sum(axis=0) / nk80            sd = np.sqrt((resp * (x[:, None] - mu) ** 2).sum(axis=0) / nk) + 1e-681            pi = nk / len(x)82            ll = float(lse.sum())83            if ll - ll_old < 1e-6:84                break85            ll_old = ll86        if best is None or ll > best[0]:87            best = (ll, mu.copy(), sd.copy(), pi.copy())88    _, mu, sd, pi = best89    order = np.argsort(mu)90    return mu[order], sd[order], pi[order]919293def gmm_means(x: np.ndarray):94    """在 log(1+W) 空间拟合，回变换成秒：E[W] = exp(μ+σ²/2) − 1。"""95    mu, sd, pi = fit_gmm1d(x)96    means = np.expm1(mu + 0.5 * sd ** 2)97    return means, mu, sd, pi9899100# --------------------------------------------------------------------------101# 主流程102# --------------------------------------------------------------------------103def main() -> None:104    print("D²Co：两成分 GMM 估 w⁺/w⁻ + 校正还原兴趣（§23.4.3）\n")105106    data = make_watch_data(n=40000, dim=16, seed=SEED)107    X, w, d, r, p = (data["x"], data["w"].astype(np.float64), data["d"],108                     data["r"], data["p"])109    idx = np.arange(len(X))110    tr, te = split(idx, len(X), 0.8)111112    # ---------------- 1) 分布长什么样 ----------------113    print("① 固定 Duration 区间内的观看时长分布：真的有两个团吗？")114    band = (d >= 40) & (d < 60)115    wb = w[band]116    rows = []117    for lo, hi in ((0, 0.5), (0.5, 2), (2, 5), (5, 10), (10, 20),118                   (20, 40), (40, 60), (60, 1e9)):119        m = (wb >= lo) & (wb < hi)120        rows.append([f"[{lo:g},{hi:g})" if hi < 1e9 else f"[{lo:g},∞)",121                     int(m.sum()), f"{m.mean():.1%}"])122    table(["观看时长区间（秒）", "样本数", "占比"], rows)123    print(f"\n  {int(band.sum()):,} 条 40–60 秒视频：{np.mean(wb == 0):.1%} 是零观看，")124    print("  另有一大团挤在 1–3 秒（划走/误触），剩下的才是不连续的「看点内容」。")125    print("  两条曲线叠在一起，这就是 w⁺/w⁻ 分解在数据上的样子。\n")126127    # ---------------- 2) 按 Duration 分组拟合 GMM ----------------128    print("② 每个 Duration 组内拟合两成分 GMM，与「按真实兴趣分组」的均值对比")129    edges = np.unique(np.quantile(d[tr], np.linspace(0, 1, 13)))130    edges[0], edges[-1] = -np.inf, np.inf131    n_b = len(edges) - 1132    kb = np.searchsorted(edges, d, side="right") - 1133134    est_p = np.full(n_b, np.nan)     # GMM 估的 w⁺135    est_n = np.full(n_b, np.nan)     # GMM 估的 w⁻136    cnt = np.zeros(n_b)137    rows = []138    for b in range(n_b):139        m = (kb == b) & (w > 0) & np.isin(idx, tr)140        cnt[b] = ((kb == b) & np.isin(idx, tr)).sum()141        if m.sum() < 50:142            continue143        means, mu, sd, pi = gmm_means(np.log1p(w[m]))144        est_n[b], est_p[b] = means[0], means[-1]145        # 真实分量：直接用生成时的 r。注意 GMM 是在 w>0 的样本上拟合的，146        # 所以对照的基准也取 w>0，否则「12% 的零观看」会把基准压低，147        # 让误差看起来比实际大。148        mt = (kb == b) & np.isin(idx, tr)149        true_p = w[mt & (r == 1) & (w > 0)].mean()150        true_n = w[mt & (r == 0)].mean()151        lo = edges[b] if np.isfinite(edges[b]) else d.min()152        hi = edges[b + 1] if np.isfinite(edges[b + 1]) else d.max()153        rows.append([f"[{lo:.0f},{hi:.0f})", int(cnt[b]),154                     f"{true_p:.1f}", f"{est_p[b]:.1f}",155                     f"{true_n:.1f}", f"{est_n[b]:.1f}",156                     f"{est_p[b] / true_p - 1:+.1%}"])157    table(["Duration 组", "样本数", "真实 w⁺", "GMM w⁺", "真实 w⁻", "GMM w⁻",158           "w⁺ 相对误差"], rows)159    print("\n  w⁻ 估得很准（几秒的噪声观看本来就好分辨），w⁺ 则是系统性高估，")160    print("  短的组 8%、最长的组 36%。原因是高成分抓到的其实是一个「看完」的")161    print("  尖峰：感兴趣的人里有一半以上看到视频结束，w 就贴在 d 上，")162    print("  log 空间里它是一个又窄又高的峰，回变换后的均值自然高于")163    print("  「把部分观看也平均进来」的 E[W|r=1]。这不是 GMM 收敛不好的抖动，")164    print("  而是「高成分 ≠ 条件均值」这个结构性差异，它会直接进分母。\n")165166    # ---------------- 3) 频次加权移动平均 ----------------167    print("③ 频次加权双向移动平均（T=2）平滑相邻 Duration 的估计")168    T = 2169    sm_p, sm_n = est_p.copy(), est_n.copy()170    for b in range(n_b):171        lo, hi = max(0, b - T), min(n_b, b + T + 1)172        sl = slice(lo, hi)173        valid = np.isfinite(est_p[sl])174        wts = cnt[sl][valid]175        if wts.sum() == 0:176            continue177        sm_p[b] = (est_p[sl][valid] * wts).sum() / wts.sum()178        sm_n[b] = (est_n[sl][valid] * wts).sum() / wts.sum()179    pick = [b for b in range(n_b) if np.isfinite(est_p[b])][:8]180    rows = []181    for b in pick:182        lo = edges[b] if np.isfinite(edges[b]) else d.min()183        hi = edges[b + 1] if np.isfinite(edges[b + 1]) else d.max()184        rows.append([f"[{lo:.0f},{hi:.0f})", f"{est_p[b]:.1f}", f"{sm_p[b]:.1f}",185                     f"{est_n[b]:.2f}", f"{sm_n[b]:.2f}"])186    table(["Duration 组", "w⁺ 原始", "w⁺ 平滑", "w⁻ 原始", "w⁻ 平滑"], rows)187    print("\n  平滑主要作用在样本量小的组：长 Duration 组只有一两千条，")188    print("  GMM 的单次估计抖动明显；借相邻组的样本量加权后曲线更连续。")189    print("  注意它只能降方差，降不了②里那种系统性低估——偏差要靠模型结构解决。\n")190191    # ---------------- 4) 校正与评价 ----------------192    print("④ 三种兴趣 Label 与真实兴趣 r 的一致性（测试集）")193    wte, dte, rte, pte = w[te], d[te], r[te], p[te]194    kte = np.searchsorted(edges, dte, side="right") - 1195    kte = np.clip(kte, 0, n_b - 1)196    wp = sm_p[kte]197    wn = sm_n[kte]198199    pcr = np.clip(wte / np.maximum(dte, 1e-6), 0, 1)200    lin = np.clip((wte - wn) / np.maximum(wp - wn, 1e-6), 0, 1)201202    def expo_label(alpha: float) -> np.ndarray:203        # (e^{αw} − e^{αw⁻})/(e^{αw⁺} − e^{αw⁻})，公因子 e^{αw⁻} 直接约掉，204        # 用 expm1 计算差值，既避免大数溢出的相减抵消，也不会因 α 小而丢精度。205        num = np.expm1(alpha * (wte - wn))206        den = np.expm1(alpha * (wp - wn))207        return np.clip(num / np.maximum(den, 1e-12), 0, 1)208209    # 对照：直接用「真实的」条件均值 E[W|r=1]、E[W|r=0] 做同一套线性校正。210    # 它不是上限——真实条件均值同样不等于 GMM 的高成分，而且它混入了211    # 「感兴趣但只看几秒」的样本，分母反而更像噪声。212    true_p = np.array([w[(kb == b) & (r == 1)].mean() if np.isfinite(est_p[b])213                       else np.nan for b in range(n_b)])214    true_n = np.array([w[(kb == b) & (r == 0)].mean() if np.isfinite(est_p[b])215                       else np.nan for b in range(n_b)])216    tp, tn = true_p[kte], true_n[kte]217    orc = np.clip((wte - tn) / np.maximum(tp - tn, 1e-6), 0, 1)218219    # 桶内 AUC：把 Duration 的影响先扣掉，只看组内排序220    def within_bin_auc(label: np.ndarray) -> float:221        tot, nsum = 0.0, 0222        for b in range(n_b):223            m = kte == b224            if m.sum() > 10 and 0 < rte[m].sum() < m.sum():225                tot += auc(rte[m], label[m]) * m.sum()226                nsum += m.sum()227        return tot / max(nsum, 1)228229    rows = [230        ["原始观看时长 w（不归一）", f"{auc(rte, wte):.3f}", "—",231         f"{within_bin_auc(wte):.3f}"],232        ["PCR：w/d", f"{auc(rte, pcr):.3f}", f"{mae(pte, pcr):.3f}",233         f"{within_bin_auc(pcr):.3f}"],234        ["D²Co(A)：线性校正", f"{auc(rte, lin):.3f}", f"{mae(pte, lin):.3f}",235         f"{within_bin_auc(lin):.3f}"],236        ["D²Co(S)：指数校正 α=0.05", f"{auc(rte, expo_label(0.05)):.3f}",237         f"{mae(pte, expo_label(0.05)):.3f}",238         f"{within_bin_auc(expo_label(0.05)):.3f}"],239        ["对照：真实条件均值做线性校正", f"{auc(rte, orc):.3f}",240         f"{mae(pte, orc):.3f}", f"{within_bin_auc(orc):.3f}"],241    ]242    table(["兴趣 Label", "AUC(r)", "MAE(p)", "桶内 AUC"], rows)243    print("\n  三个数各自说明一件事：")244    print("    · AUC(r)：D²Co(A) 0.862 vs PCR 0.830，+0.032。改善来自")245    print("      「先减掉 w⁻ 再加权」这个结构——它把「光划走」的样本压到 0 附近，")246    print("      让真正看了内容的人和它们分开；")247    print("    · MAE(p)：0.358 → 0.358，几乎没动。Label 更适合排序，")248    print("      但直接把它当兴趣概率用（比如拿去算 pCTR）并不比 PCR 准，")249    print("      因为分母 w⁺ 高估，Label 整体被压缩向 0；")250    print("    · 桶内 AUC：0.829 → 0.860，和整体 AUC 几乎同幅提升。")251    print("      说明增益不是「Duration 信息被扣掉了」的假象，")252    print("      而是同一 Duration 内的排序确实变好了。")253    print("  对照行更值得注意：用「真实条件均值」做同一套线性校正只有 0.859，")254    print("  还略低于 GMM 版本的 0.862。条件均值把「感兴趣但只看三秒」")255    print("  也平均进来，分母被稀释；GMM 的高成分更接近「看完」的尖峰，")256    print("  恰好把这类样本的贡献压了下去。所以「估得越像 E[W|r=1] 越好」")257    print("  在这里并不成立——高成分不是条件均值，而是条件均值里的一支。")258    print("  作为参照，直接用生成用的兴趣概率 p 当分数，AUC(r) 只有 0.817。")259    print("  这不构成「上限」——r 是按 p 抽样出来的，而 w 又是在 r 的基础上")260    print("  再生成的，所以 w 天然比 p 多带了信息。有意义的比较只有同一行内")261    print("  横向看：w / PCR / D²Co 三者用的是同一份观测，D²Co 多出 0.032。\n")262    print("  正文说 w⁺ 和 w⁻ 接近时分母小、误差被放大；本例里 w⁻ 稳定在 2.5 秒，")263    print("  w⁺ 从十几秒到上百秒，分母并不小。真正放大误差的是反方向：")264    print("  高成分高估 8%~36%，且长视频组更严重，于是长视频的 Label 被压得")265    print("  更靠近 0，排序时它们的相对差异被缩小——但方向没变，所以 AUC 仍涨。\n")266267    # ---------------- 5) α 的稳健性 ----------------268    print("⑤ 指数校正的 α：它控制对估计误差的敏感度")269    rows = []270    for alpha in (0.01, 0.02, 0.05, 0.1, 0.2):271        lab = expo_label(alpha)272        rows.append([f"{alpha:.2f}", f"{auc(rte, lab):.3f}", f"{mae(pte, lab):.3f}",273                     f"{lab.std():.3f}", f"{np.mean((lab < 1e-6) | (lab > 1 - 1e-6)):.1%}"])274    table(["α", "AUC(r)", "MAE(p)", "Label 标准差", "落在端点的比例"], rows)275    print("\n  α→0 时 expm1(αx)/α → x，指数校正精确退化成线性校正，")276    print("  所以 0.862 这条基线不会因为换公式而改变。α 变大后，")277    print("  落在两端的样本变多（56% → 64%）、标准差反而下降（0.436 → 0.375）：")278    print("  中间地带被推到 0 或 1，分布从「连续」变成「两端堆积」。")279    print("  代价在 MAE 上：α=0.05 之前 0.356，α=0.2 时劣化到 0.376。")280    print("  也就是说 α 在这里不是提升排序的旋钮，只是改变 Label 的离散程度；")281    print("  论文按数据分布调它，换数据集不能照搬。\n")282283    print("⑥ 一句话：D²Co 把观看时长拆成「兴趣 × 偏差 + (1−兴趣) × 噪声」，")284    print("   再用同 Duration 下的大量样本把两个潜在分量估出来反解兴趣。")285    print("   它的收益主要体现在排序（AUC +0.03）而不是概率校准（MAE 不动），")286    print("   而高成分估的是「看完」那一支、不是条件均值——这一点决定了")287    print("   它的 Label 会比「用真实条件均值」更锐利，也更容易受长尾影响。")288289290if __name__ == "__main__":291    main()