""" ============================================================================= 第 12 讲 · 代码 1：贝叶斯推断——Beta–Bernoulli 共轭先验 ============================================================================= 对应讲义：lectures/lecture12_thompson_sampling.md 第 1–2 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 第 11 讲的 UCB 用"置信区间"表示不确定性。贝叶斯方法换了一种说法： 把未知的成功率 μ 看成一个**随机变量**，用一个概率分布描述"我相信 μ 是多少"。 - 看数据之前的信念 = **先验** p(μ) - 看到数据之后的信念 = **后验** p(μ | 数据) ∝ p(数据 | μ) · p(μ) （贝叶斯公式） 对伯努利奖励（成功/失败），如果先验是 Beta 分布，后验**还是** Beta 分布，而且更新只需要加法： 先验 Beta(α, β) + 观察到 s 次成功、f 次失败 → 后验 Beta(α + s, β + f) 这叫**共轭先验**。α、β 可以理解为"虚拟的成功次数和失败次数"。 本文件： 1. 一步一步看后验如何随数据变化（用字符画出密度形状） 2. 验证"共轭"：用最笨的数值方法（网格上逐点做贝叶斯公式）算后验，与 Beta 公式对照 3. 先验的影响：先验相当于多少次虚拟观察？数据多了之后先验还重要吗？ 4. 可信区间是否"可信"：先验正确时，90% 可信区间恰好 90% 的时候包含真值；先验错误时呢？ ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- cd code/lecture12 python 01_beta_bernoulli.py （约 1 秒） """ import numpy as np from bayes_bandits import beta_interval, beta_mean_std, beta_pdf, pad, sparkline np.set_printoptions(precision=4, suppress=True) def grid_posterior(prior_pdf, successes, failures, n_grid=2001): """ 用最直接的数值方法计算后验：在 [0,1] 上取 n_grid 个点，逐点计算"先验 × 似然"，再归一化。 【参数】 prior_pdf : 函数，输入 μ 的数组，返回先验密度 successes, failures : 成功、失败次数 【返回】(网格点 x, 后验密度) 【为什么要这样做？】这是贝叶斯公式最字面的实现，不依赖任何"共轭"技巧，适用于任何先验。 它的缺点是只能用于一两个未知数——未知数多了，网格点的数量会指数爆炸。 共轭先验的价值就在于：对特定的"先验 + 似然"组合，后验有现成的公式，不需要数值积分。 【边界】似然 μ^s (1-μ)^f 在 s、f 很大时会下溢为 0，所以在对数空间里算。 """ x = np.linspace(1e-6, 1 - 1e-6, n_grid) log_post = np.log(prior_pdf(x)) + successes * np.log(x) + failures * np.log(1 - x) post = np.exp(log_post - log_post.max()) # 减去最大值再取 exp，防止下溢/上溢 post /= np.trapezoid(post, x) if hasattr(np, "trapezoid") else np.trapz(post, x) # 【Python 小知识】np.trapezoid 用梯形法则做数值积分（numpy 2.0 之前叫 np.trapz） return x, post def main(): rng = np.random.default_rng(0) # --------------------------------------------------------------------- print("=" * 78) print("演示 1：一步一步的贝叶斯更新（真实成功率 μ = 0.7，先验 Beta(1, 1) = 均匀分布）") print("=" * 78) data = (rng.random(200) < 0.7).astype(int) print(f" 前 20 个观察：{data[:20]}（1 = 成功，0 = 失败）") print(" 看过的数据 后验 均值 标准差 90% 可信区间 密度形状（横轴 0 → 1）") for n in [0, 1, 2, 5, 10, 50, 200]: s = int(data[:n].sum()) a, b = 1 + s, 1 + n - s mean, std = beta_mean_std(a, b) lo, hi = beta_interval(a, b, rng=rng) print(f" n = {n:3d} Beta({a:3d},{b:3d}) {mean:.3f} {std:.3f} [{lo:.3f}, {hi:.3f}] |{sparkline(a, b)}|") print(" 观察：") print(" - 没有数据时，后验就是先验（均匀分布：μ 取 0～1 之间任何值都一样可能）。") print(" - 第一个观察是\"成功\"：密度变成一条斜线（μ 越大越可能），但仍然很宽。") print(" - 数据越多，后验越窄。标准差大约按 1/sqrt(n) 缩小——") print(" 与第 11 讲霍夫丁置信区间的宽度 ∝ 1/sqrt(n) 一致：两种语言描述的是同一件事。") print(" - 注意 n = 200 时，90% 可信区间 [0.577, 0.688] **没有**包含真实值 0.7：这 200 个观察中恰好只有 63% 是成功。") print(" 这并不矛盾——\"90% 可信区间\"本来就有大约 10% 的时候不包含真值（演示 4 会验证这一点）。") print(" 后验描述的是\"根据这些数据应该相信什么\"，数据本身不走运时，后验也会跟着偏。") print(" - 更新只需要加法：成功一次 α + 1，失败一次 β + 1。这就是 Beta 分布作为\"共轭先验\"的便利之处。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 2：验证共轭——网格数值计算的后验 vs Beta 公式") print("=" * 78) for (pa, pb), (s, f) in [((1, 1), (7, 3)), ((2, 8), (30, 70)), ((5, 5), (0, 4))]: x, post = grid_posterior(lambda m: beta_pdf(m, pa, pb), s, f) formula = beta_pdf(x, pa + s, pb + f) err = np.max(np.abs(post - formula)) / np.max(formula) print(f" 先验 Beta({pa},{pb})，观察 {s} 次成功 {f} 次失败 → 公式给出 Beta({pa + s},{pb + f})；" f"与网格计算的最大相对差异 = {err:.1e}") seq = rng.random(30) < 0.4 a, b = 1.0, 1.0 for obs in seq: # 一个一个地更新 a, b = a + obs, b + (1 - obs) print(f" 顺序无关：30 个观察逐个更新得到 Beta({a:.0f},{b:.0f})；一次性更新得到 " f"Beta({1 + seq.sum()},{1 + (~seq).sum()})——完全相同。") print(" → 公式与最笨的数值方法完全一致（差异在 10⁻¹³ 量级，只是浮点数的舍入误差）。") print(" 而且后验只取决于成功和失败的**次数**，与顺序无关——所以老虎机算法只需要记录 (次数, 奖励总和)。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 3：先验的影响——一个新广告的点击率（真实值 3%）") print("=" * 78) priors = [("均匀：Beta(1,1)", 1, 1), ("了解行情：Beta(2,48)", 2, 48), ("错误但不固执：Beta(1,3)", 1, 3), ("错误且固执：Beta(50,50)", 50, 50)] clicks = rng.random(5000) < 0.03 print(" " + pad("先验（先验均值）", 34) + " n=0 n=10 n=100 n=1000 n=5000 ← 后验均值") for name, pa, pb in priors: row = [] for n in [0, 10, 100, 1000, 5000]: s = clicks[:n].sum() row.append((pa + s) / (pa + pb + n)) print(" " + pad(f"{name}（{pa / (pa + pb):.3f}）", 34) + "".join(f"{v:9.4f}" for v in row)) mle = [clicks[:n].mean() if n > 0 else float("nan") for n in [0, 10, 100, 1000, 5000]] print(" " + pad("不用先验：点击次数 / 展示次数", 34) + "".join(f"{v:9.4f}" for v in mle)) print(f" （前 10 次展示中有 {clicks[:10].sum()} 次点击，前 100 次中有 {clicks[:100].sum()} 次。）") print(" 观察：") print(" - 后验均值 = (α + 成功) / (α + β + n)：先验相当于 α + β 次\"虚拟观察\"。") print(" Beta(50,50) 相当于\"已经看过 100 次、一半被点击\"——需要几千次真实数据才能把它纠正过来。") print(" - \"了解行情\"的先验在数据很少时最准：n = 10 时恰好有 1 次点击，直接用比例估计是 10%（真实值的 3 倍多），") print(" 均匀先验给出 16.7%，而 Beta(2,48) 给出 5%。这就是先验知识的价值：**数据少时，先验很重要**。") print(" - 数据足够多时（n = 5000），除了最固执的先验，所有估计都收敛到真实值附近——**数据终将压倒先验**。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 4：可信区间可信吗？（真实点击率从 Beta(2,48) 中抽取，每次看 50 次展示；重复 2000 次）") print("=" * 78) n_trials, n_obs = 2000, 50 mu = rng.beta(2, 48, size=n_trials) # 每次试验一个新广告 s = rng.binomial(n_obs, mu) print(" " + pad("用来计算后验的先验", 26) + " 90% 区间包含真值的比例 区间平均宽度") for name, pa, pb in [("正确的先验 Beta(2,48)", 2, 48), ("均匀先验 Beta(1,1)", 1, 1), ("错误且固执 Beta(50,50)", 50, 50)]: samples = rng.beta(pa + s[:, None], pb + n_obs - s[:, None], size=(n_trials, 4000)) lo, hi = np.quantile(samples, [0.05, 0.95], axis=1) cover = np.mean((lo <= mu) & (mu <= hi)) print(" " + pad(name, 26) + f" {cover * 100:5.1f}% {np.mean(hi - lo):.4f}") half = np.sqrt(np.log(2 / 0.1) / (2 * n_obs)) # 霍夫丁 90% 区间的半宽 mean = s / n_obs cover_h = np.mean((mean - half <= mu) & (mu <= mean + half)) print(" " + pad("霍夫丁区间（不用先验）", 26) + f" {cover_h * 100:5.1f}% {2 * half:.4f}") print(" 观察：") print(" - 先验正确时，90% 可信区间恰好在约 90% 的情况下包含真值——贝叶斯推断是**校准**的。") print(" 注意这个\"90%\"是对\"从先验中抽出的所有问题\"平均而言的，不是对某一个固定的 μ。") print(" - 均匀先验：仍然大致可用，但区间偏宽（它不知道点击率通常很低）。") print(" - 先验错误且固执：区间被先验牢牢拉在 0.27～0.41 一带（真实点击率只有百分之几），**从不**包含真值——") print(" **错误的先验会让你\"自信地犯错\"**：后验看起来和正常的后验一样，你从它本身看不出任何问题。") print(" - 霍夫丁区间对任何分布都成立，所以总是\"覆盖\"得很好，但宽得多（这里比真实点击率的范围还宽）。") print(" 这就是频率学派与贝叶斯方法的根本取舍：前者不需要先验、有最坏情况保证，但保守；") print(" 后者能利用先验知识、更高效，但依赖先验的正确性。") if __name__ == "__main__": main()