""" ============================================================================= 第 5 讲 · 代码 1：最大似然的三件套 —— 似然、梯度、海森矩阵，以及不确定性 ============================================================================= 对应讲义：lectures/lecture05_mle.md 第 2、3、6 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 前几讲一直把 fit_bradley_terry 当作黑盒。这个文件把它打开： 演示 1：最简单的情形 —— 只有两个选项，a 赢了 70 次、输了 30 次。似然函数长什么样？ 最大值在哪里（有封闭解）？曲线的"尖锐程度"怎样告诉我们估计有多确定？ 演示 2：梯度 = "观察到的胜场 − 模型预期的胜场"。用数值差分验证解析梯度。 演示 3：海森矩阵 = 负的加权图拉普拉斯矩阵。它的特征值说明了两件事： ① 有一个特征值是 0（所有分数加同一个常数，似然不变）；② 其余都是负的（目标函数是凹的， 只有一个最高点）。比较图的"连通得好不好"决定了最小的那个非零特征值，也就决定了估计有多准。 演示 4：不确定性 —— 用海森矩阵算出的标准误，和"重复做 300 次实验"的实际波动、自助法的结果对比； 95% 置信区间是不是真的有 95% 的时候包含真值？ ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 01_likelihood_gradient.py # 约 3 秒 需要 numpy（第 5 讲起使用）。只依赖本目录的 bt_numpy.py。 """ import math import numpy as np from bt_numpy import (bar, covariance, diff_se, fit_newton, gradient, hessian, log_sigmoid, objective, pad, print_header, random_pairs, sigmoid, simulate_bt) # ============================================================================= # 演示 1：两个选项 # ============================================================================= def demo_two_items(): print_header("演示 1：两个选项，a 赢 70 次、输 30 次 —— 似然函数长什么样？") wins, losses = 70, 30 n = wins + losses print("设 Δ = s_a − s_b。对数似然 = 70·log σ(Δ) + 30·log σ(−Δ)。在一串 Δ 上算出来（画成横条，越长越大）：\n") deltas = [0.0, 0.25, 0.5, 0.75, 0.85, 1.0, 1.25, 1.5, 2.0] lls = [wins * float(log_sigmoid(d)) + losses * float(log_sigmoid(-d)) for d in deltas] top, low = max(lls), min(lls) for d, ll in zip(deltas, lls): print(f" Δ = {d:>4.2f} 对数似然 {ll:>8.3f} {bar((ll - low) / (top - low))}") mle = math.log(wins / losses) p = wins / n se = 1.0 / math.sqrt(n * p * (1 - p)) print(f""" 最大似然估计：令导数 70·(1 − σ(Δ)) − 30·σ(Δ) = 0，得 σ(Δ̂) = 0.7，即 Δ̂ = log(70/30) = {mle:.3f}。 曲率（二阶导数）= −n·p(1−p) = −{n * p * (1 - p):.1f}，标准误 ≈ 1/√{n * p * (1 - p):.1f} = {se:.3f}， 所以 Δ 的 95% 置信区间大约是 {mle:.3f} ± 1.96 × {se:.3f} = [{mle - 1.96 * se:.3f}, {mle + 1.96 * se:.3f}]。 如果数据多 100 倍（7000 胜 3000 负），最大值的位置不变，曲线尖锐 100 倍，标准误缩小 10 倍。""") print(""" 观察： · 最大似然估计有一个漂亮的直觉：让模型预测的胜率 σ(Δ̂) 等于观察到的胜率 70%。 · 曲线在最高点附近越“尖”，数据越能确定 Δ；越“平”，越多的 Δ 都差不多一样好，越不确定。 “尖锐程度”就是二阶导数（曲率），这是第 7 讲 Fisher 信息的直觉来源。 · 注意 p(1 − p) 这个因子：p = 0.5 时最大。实力接近的两个选项，比一次学到的东西最多。""") # ============================================================================= # 演示 2：梯度 = 观察到的 − 预期的 # ============================================================================= def demo_gradient(rng): print_header("演示 2：梯度 = 观察到的胜场 − 模型预期的胜场") M = 5 true_s = np.array([1.0, 0.5, 0.0, -0.5, -1.0]) i, j, y = simulate_bt(true_s, 400, rng) s = np.zeros(M) # 在“所有分数都是 0”这一点上看梯度 p = sigmoid(s[i] - s[j]) winner = np.where(y == 1, i, j) observed = np.bincount(winner, minlength=M) expected = np.bincount(i, weights=p, minlength=M) + np.bincount(j, weights=1 - p, minlength=M) g = gradient(s, i, j, y) * len(y) # 乘以 N，得到"总"梯度（不是平均） print("在 s = 0（模型认为每场都是五五开）这一点：\n") print(f" {pad('选项', 4, 'right')} {pad('真实分数', 9, 'right')} {pad('实际胜场', 9, 'right')} " f"{pad('预期胜场', 9, 'right')} {pad('梯度', 8, 'right')}") for k in range(M): print(f" {k:>4d} {true_s[k]:>9.1f} {observed[k]:>9d} {expected[k]:>9.1f} {g[k]:>8.1f}") # 数值验证：中心差分 s0 = rng.normal(0, 0.5, size=M) h = 1e-6 num = np.array([(objective(s0 + h * np.eye(M)[k], i, j, y) - objective(s0 - h * np.eye(M)[k], i, j, y)) / (2 * h) for k in range(M)]) ana = gradient(s0, i, j, y) print(f"\n 在一个随机的点上，解析梯度和数值差分梯度的最大差别：{np.max(np.abs(num - ana)):.1e}") print(""" 观察： · 梯度的第 k 个分量 = 选项 k 实际赢的场数 − 模型（按当前分数）预期它赢的场数。赢得比预期多，就把分数往上调。 这就是第 2 讲“被选中的次数 − 预期的次数”在两两比较上的样子。 · 梯度为 0 的点（最大似然估计）满足：**每个选项的预期胜场 = 实际胜场**。这叫“矩匹配”，是指数族模型 最大似然的普遍性质（逻辑回归也是这样）。Zermelo 1929 年就是从这个方程出发解出排名的。 · 所有分量加起来是 0：一场比赛有一胜一负，总的实际胜场 = 总的预期胜场。这是“平移不变”在梯度上的体现。 · 数值差分（每个方向挪一点点看函数变了多少）和解析公式吻合到 1e-10 量级：写梯度代码时，永远先做这个检查。""") # ============================================================================= # 演示 3：海森矩阵与比较图 # ============================================================================= def chain_pairs(M, n, rng): """“链式”设计：每个选项只和编号相邻的选项（差 1 或 2）比较。""" i = rng.integers(0, M, size=n) step = rng.choice(np.array([-2, -1, 1, 2]), size=n) j = np.clip(i + step, 0, M - 1) keep = i != j return i[keep], j[keep] def demo_hessian(rng): print_header("演示 3：海森矩阵 = 负的加权图拉普拉斯 —— 凹、平移不变、连通性决定精度") M = 30 true_s = np.linspace(-1.5, 1.5, M) designs = {"随机配对（连通得好）": random_pairs(M, 1500, rng), "链式（只和邻居比）": chain_pairs(M, 1500, rng)} print(f"{M} 个选项，1500 条比较，两种“谁和谁比”的设计。在真实分数处计算海森矩阵 H 的特征值（−H 的特征值）：\n") print(f" {pad('设计', 22)} {pad('最小特征值', 12, 'right')} {pad('第二小（代数连通度）', 22, 'right')} " f"{pad('最大特征值', 12, 'right')} {pad('s_0 − s_29 的标准误', 20, 'right')}") for name, (i, j) in designs.items(): _, _, y = simulate_bt(true_s, len(i), rng, pairs=(i, j)) H = hessian(true_s, i, j, y) * len(y) ev = np.sort(np.linalg.eigvalsh(-H)) se = diff_se(covariance(true_s, i, j, y), 0, M - 1) print(f" {pad(name, 22)} {ev[0]:>12.2e} {ev[1]:>22.3f} {ev[-1]:>12.2f} {se:>20.3f}") print(""" 观察： · 最小的特征值是 0（只差浮点误差），对应的特征向量是全 1 向量：所有分数加同一个常数，似然不变。 这就是“只有差值可识别”在线性代数上的样子。加 L2 正则化（λ > 0）会把它变成 λ，让解唯一。 · 其余特征值都是正的（−H 半正定）：对数似然是**凹**函数，只有一个最高点，梯度上升不会被困在“局部最优”。 这是 BT（以及所有逻辑回归）好优化的根本原因 —— 第 6 讲的神经网络奖励模型就没有这么幸运了。 · 第二小的特征值叫比较图的“代数连通度”（Fiedler 值）。链式设计里，相隔很远的两个选项只能通过一长串 中间人间接比较，误差层层累积：代数连通度小了几十倍，s_0 − s_29 的标准误大了将近三倍。 数据量一样，**谁和谁比**决定了能学多准 —— 这是第 7 讲（实验设计）的出发点。""") # ============================================================================= # 演示 4：不确定性 # ============================================================================= def demo_uncertainty(rng, reps=300): print_header("演示 4：标准误靠不靠得住？—— 公式 vs 重复实验 vs 自助法，以及置信区间的覆盖率") M = 10 true_s = np.linspace(-1.0, 1.0, M) true_s -= true_s.mean() n = 600 a, b = 0, M - 1 estimates, formula_se, covered = [], [], 0 for _ in range(reps): i, j, y = simulate_bt(true_s, n, rng) s, _ = fit_newton(i, j, y, M) se = diff_se(covariance(s, i, j, y), a, b) est = s[a] - s[b] estimates.append(est) formula_se.append(se) truth = true_s[a] - true_s[b] if abs(est - truth) <= 1.96 * se: covered += 1 # 自助法：只用最后一份数据，有放回地重抽 200 次 boots = [] for _ in range(200): idx = rng.integers(0, n, size=n) s, _ = fit_newton(i[idx], j[idx], y[idx], M) boots.append(s[a] - s[b]) print(f"{M} 个选项，每次实验 {n} 条比较，关心 s_0 − s_9（真实值 {true_s[a] - true_s[b]:.2f}）。重复 {reps} 次实验：\n") print(f" 估计值的平均： {np.mean(estimates):.3f}") print(f" 估计值的实际标准差： {np.std(estimates):.3f} （“真正的”不确定性）") print(f" 公式给出的标准误（平均）：{np.mean(formula_se):.3f} （来自负海森矩阵的伪逆）") print(f" 自助法的标准误（一份数据）：{np.std(boots):.3f}") print(f" 95% 置信区间（估计 ± 1.96 × 标准误）包含真值的比例：{covered / reps:.1%}") print(""" 观察： · 三种不确定性的估计基本一致：数据足够多时，最大似然估计近似服从正态分布，协方差 ≈ 负海森矩阵的逆 （这是统计学的经典结论，第 7 讲讲它和 Fisher 信息的关系）。 · 95% 置信区间确实大约 95% 的时候包含真值 —— 这叫“覆盖率”。数据很少、或者某些选项几乎总赢 / 总输时， 正态近似变差，覆盖率会低于 95%，这时自助法或贝叶斯方法（第 6 讲）更可靠。 · 估计值的平均略大于真实值（绝对值偏大）：最大似然在小样本下有轻微的“往外推”的偏差，第 4 讲的附带参数问题 是它的极端情形；L2 正则化会把它往回拉（第 4 节）。""") def main(): rng = np.random.default_rng(1952) demo_two_items() demo_gradient(rng) demo_hessian(rng) demo_uncertainty(rng) if __name__ == "__main__": main()