""" ============================================================================= 第 12 讲 · 代码 2：汤普森采样——从后验中抽一个"可能的世界"，然后在其中做最优选择 ============================================================================= 对应讲义：lectures/lecture12_thompson_sampling.md 第 3–5 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 有了每个臂的后验，怎么选臂？ - 选后验均值最高的（贝叶斯贪心）：把不确定性丢掉了，又回到贪心，会卡住。 - 选后验"上界"最高的：就是 UCB 的贝叶斯版本（Bayes-UCB）。 - **汤普森采样**：从每个臂的后验中各抽一个样本，选样本最大的臂。 汤普森采样的效果是**概率匹配**：每个臂被选中的概率 = 根据目前数据它是最好的臂的后验概率。 - 几乎肯定是最好的臂 → 几乎总是选它（利用）； - 有一定可能是最好的臂 → 按这个可能性去试一试（探索）； - 几乎不可能是最好的臂 → 几乎不再浪费机会。 本文件： 1. 验证概率匹配：汤普森采样的选择频率 = 后验"是最好的臂"的概率（后者用数值积分独立计算） 2. 一次实验中，各臂的后验如何演化 3. 频率学派的遗憾：汤普森采样 vs KL-UCB vs UCB1 vs 贝叶斯贪心 4. 贝叶斯遗憾与先验：正确的先验有多大帮助？错误的先验（过于乐观 / 过于悲观）有什么后果？ 5. 批量反馈：奖励每隔一段时间才汇总一次时，随机化的汤普森采样比确定性的 UCB 稳健得多 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- cd code/lecture12 python 02_thompson_sampling.py （约 80 秒；KL-UCB 比较慢） """ import numpy as np from bayes_bandits import (CTR_MEANS, DEFAULT_MEANS, bayes_greedy, beta_pdf, kl_ucb, lai_robbins_constant, pad, simulate, sparkline, thompson, ucb1) np.set_printoptions(precision=3, suppress=True) def prob_best_numeric(alphas, betas, n_grid=20001): """ 用数值积分计算"臂 a 是最好的臂"的后验概率： P(a 最好) = ∫ p_a(x) · Π_{b ≠ a} F_b(x) dx 其中 p_a 是臂 a 的后验密度，F_b 是臂 b 的后验累积分布函数（"μ_b ≤ x 的概率"）。 【直觉】臂 a 最好 ⟺ μ_a = x 且其他所有臂都小于 x；对所有可能的 x 求和（积分）。 【实现】在 [0,1] 的细网格上算密度，累积求和得到 CDF，再做数值积分。 这个计算完全不依赖随机抽样——用来**独立地**验证汤普森采样的选择频率。 """ x = np.linspace(0, 1, n_grid) dx = x[1] - x[0] pdfs = np.array([beta_pdf(x, a, b) for a, b in zip(alphas, betas)]) cdfs = np.cumsum(pdfs, axis=1) * dx cdfs = np.minimum(cdfs / cdfs[:, -1:], 1.0) # 归一化，消除离散误差 out = [] for a in range(len(alphas)): others = np.prod(np.delete(cdfs, a, axis=0), axis=0) out.append(np.sum(pdfs[a] * others) * dx) out = np.array(out) return out / out.sum() def run_single_ts(means, T, snapshots, seed): """ 单独跑一次汤普森采样（普通 for 循环，方便观察），在 snapshots 时刻记录每个臂的后验参数 (α, β)。 【返回】列表：(t, α 数组, β 数组) """ rng = np.random.default_rng(seed) K = len(means) alpha, beta = np.ones(K), np.ones(K) out = [] for t in range(1, T + 1): a = int(np.argmax(rng.beta(alpha, beta))) r = rng.random() < means[a] alpha[a] += r beta[a] += 1 - r if t in snapshots: out.append((t, alpha.copy(), beta.copy())) return out def main(): # --------------------------------------------------------------------- print("=" * 78) print("演示 1：概率匹配——汤普森采样选每个臂的频率 = 它是最好的臂的后验概率") print("=" * 78) t_snap, alpha, beta = run_single_ts(DEFAULT_MEANS, 60, {60}, seed=1)[0] counts, succ = alpha + beta - 2, alpha - 1 print(f" 取一次汤普森采样实验在第 {t_snap} 步的状态（真实 μ = {DEFAULT_MEANS}）：") p_best = prob_best_numeric(alpha, beta) rng = np.random.default_rng(2) picks = np.argmax(rng.beta(alpha, beta, size=(200000, len(alpha))), axis=1) freq = np.bincount(picks, minlength=len(alpha)) / len(picks) ucb = succ / counts + np.sqrt(2 * np.log(t_snap) / counts) print(" 臂 拉的次数 成功 后验 后验均值 P(它是最好的臂) 汤普森采样选它的频率") print(" （数值积分） （抽 20 万次）") for a in range(len(alpha)): print(f" {a} {int(counts[a]):4d} {int(succ[a]):4d} Beta({int(alpha[a]):2d},{int(beta[a]):2d})" f" {alpha[a] / (alpha[a] + beta[a]):.3f} {p_best[a]:.4f} {freq[a]:.4f}") post_mean = alpha / (alpha + beta) print(f" 贝叶斯贪心会选：臂 {int(np.argmax(post_mean))}（后验均值最高）——每次都选它，不会探索其他臂。") print(f" UCB1（c = 2）会选：臂 {int(np.argmax(ucb))}（上界最高）——确定性的选择。") print(" 观察：") print(" - 两列概率几乎完全相同（差异只来自抽样误差）：汤普森采样**恰好**按\"是最好的臂的概率\"来分配机会。") print(" - 这就是\"概率匹配\"。它不需要显式地计算这些概率——只要抽一次样本、取最大，就自动实现了。") print(" - 注意：汤普森采样是**随机**的策略，同样的数据下每次可能选不同的臂；UCB 和贝叶斯贪心是确定性的。") print(" - 三种方法在这个状态下的选择各不相同：贝叶斯贪心只看\"最可能的值\"，选臂 4；") print(" UCB1 选了只拉过 2 次、2 次都失败的臂 1——它的上界最高，因为它最不确定；") print(" 汤普森采样则大部分时候选臂 4，偶尔（按概率）试试其他臂，选臂 1 的概率只有约 2%。") print(" UCB 按\"可能有多好\"探索，汤普森采样按\"有多大可能是最好的\"探索。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 2：一次实验中后验的演化（每行是一个臂的后验密度形状，横轴 0 → 1）") print("=" * 78) for t, alpha, beta in run_single_ts(DEFAULT_MEANS, 10000, {10, 100, 1000, 10000}, seed=3): print(f"\n t = {t}：") for a in range(len(alpha)): n = int(alpha[a] + beta[a] - 2) print(f" 臂 {a}（μ = {DEFAULT_MEANS[a]:.1f}，拉了 {n:5d} 次）|{sparkline(alpha[a], beta[a])}|") print("\n 观察：") print(" - 最好的臂（臂 4）拉得最多，后验越来越窄，集中在 0.7 附近。") print(" - 差的臂拉得少，后验一直比较宽——汤普森采样不在乎把差的臂估计得准不准，") print(" 只要它们的后验\"几乎不可能超过\"最好的臂就够了。这一点与 UCB 一样：探索是为了决策，不是为了精确估计。") print(" - 与最好的臂最接近的臂 3 被拉得比其他次优臂多：它是最难排除的\"竞争者\"。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 3：频率学派的遗憾（固定的老虎机，T = 10⁵，重复 100 次）") print("=" * 78) cps = [1000, 10000, 100000] agents = [("汤普森采样（均匀先验）", thompson()), ("KL-UCB", kl_ucb()), ("UCB1（c = 0.5）", ucb1(0.5)), ("UCB1（c = 2）", ucb1(2.0)), ("贝叶斯贪心（后验均值）", bayes_greedy())] for label, m in [("默认老虎机", DEFAULT_MEANS), ("点击率老虎机", CTR_MEANS)]: C = lai_robbins_constant(m) print(f"\n {label}：μ = {m}，Lai–Robbins 常数 C = {C:.1f}") print(" " + pad("策略", 24) + " R(10³) R(10⁴) R(10⁵) R(10⁵)/ln(10⁵) 最差 1% 实验") for name, pol in agents: res = simulate(pol, m, T=100000, n_runs=100, seed=7, checkpoints=cps) r = res["regret"].mean(axis=0) print(" " + pad(name, 24) + "".join(f"{x:9.1f}" for x in r) + f" {r[2] / np.log(1e5):6.1f} {np.percentile(res['regret'][:, -1], 99):9.1f}") print("\n 观察：") print(" - 汤普森采样在两个老虎机上的遗憾都是最小的，比渐近最优的 KL-UCB 还小约三分之一。") print(" 理论上两者渐近地同样最优（Kaufmann et al., 2012；Agrawal & Goyal, 2012）；") print(" 在有限的 T 下，汤普森采样往往表现更好——这是它在工业界流行的主要原因之一（另一个原因是简单）。") print(" - 它的 R(T)/ln T 甚至低于 C：和第 11 讲一样，下界是 T → ∞ 的渐近结论，有限的 T 下可以低于它。") print(" - 贝叶斯贪心用的是**同一个后验**，只是取均值而不是抽样本——遗憾却线性增长，最差的实验卡死在次优臂上。") print(" **抽样本这一步，就是全部的探索**。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 4：贝叶斯遗憾与先验（500 个随机的点击率老虎机，每个 10 个臂，T = 10⁴）") print("=" * 78) rng = np.random.default_rng(0) M = rng.beta(2, 48, size=(500, 10)) # 每个老虎机的真实点击率从 Beta(2,48) 抽取（平均 4%） print(f" 真实点击率从 Beta(2,48) 中抽取：平均 {M.mean():.3f}，每个老虎机中最好的臂平均 {M.max(axis=1).mean():.3f}") agents = [("TS，正确的先验 Beta(2,48)", thompson(2, 48)), ("TS，均匀先验 Beta(1,1)", thompson(1, 1)), ("TS，过于乐观 Beta(20,20)", thompson(20, 20)), ("TS，过于悲观 Beta(1,500)", thompson(1, 500)), ("KL-UCB（不用先验）", kl_ucb()), ("UCB1（c = 0.5）", ucb1(0.5))] print(" " + pad("策略", 28) + " R(10²) R(10³) R(10⁴) 最差 1% 实验") for name, pol in agents: res = simulate(pol, M, T=10000, n_runs=500, seed=6, checkpoints=[100, 1000, 10000]) r = res["regret"].mean(axis=0) print(" " + pad(name, 28) + "".join(f"{x:9.1f}" for x in r) + f" {np.percentile(res['regret'][:, -1], 99):7.1f}") print(" 观察：") print(" - \"贝叶斯遗憾\" = 对\"从先验中抽出的问题\"平均的遗憾。在这个设定下，用正确先验的汤普森采样是最好的。") print(" - 均匀先验的汤普森采样稍差，但仍然优于 KL-UCB、远优于 UCB1——即使先验不准，汤普森采样依然很强。") print(" - **过于乐观**的先验（以为点击率在 50% 左右）：每个臂都要先\"失败\"几百次，后验才会降到真实水平——") print(" 探索过多，遗憾很大，但很稳定（和第 11 讲 UCB 的 c 太大时一样）。") print(" - **过于悲观**的先验（以为点击率只有 0.2%）：前期看起来很好（很少探索），但探索不足——") print(" 最差的实验遗憾很大，和贪心\"卡住\"一样。错误的先验，就是错误的探索强度。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 5：批量反馈——每 B 步才汇总一次奖励（默认老虎机，T = 20000，重复 100 次）") print("=" * 78) agents = [("汤普森采样", thompson()), ("KL-UCB", kl_ucb()), ("UCB1（c = 0.5）", ucb1(0.5))] print(" " + pad("策略", 16) + "".join(f" B = {B:<5d}" for B in [1, 10, 100, 1000])) for name, pol in agents: vals = [simulate(pol, DEFAULT_MEANS, T=20000, n_runs=100, seed=5, batch_size=B)["regret"].mean() for B in [1, 10, 100, 1000]] print(" " + pad(name, 16) + "".join(f"{v:11.1f}" for v in vals)) print(" 观察：") print(" - 真实系统中，奖励往往不是立刻反馈的：日志每小时汇总一次，模型每天更新一次。") print(" - B = 100 时，UCB 类方法的遗憾明显上升，汤普森采样几乎不受影响。原因：") print(" UCB 是**确定性**的——同一批里看到的数据都一样，它会把整批 100 次全部分给同一个臂；") print(" 汤普森采样是**随机**的——同一批里每次抽样都不同，自然地把机会分散给几个\"可能最好\"的臂。") print(" - B = 1000 时（相当于一次只更新 20 次），所有方法都明显变差：信息更新得太少，谁都没办法。") print(" - 这是 Chapelle & Li (2011) 在雅虎新闻推荐中发现的：汤普森采样在真实系统中常常比 UCB 更好。") if __name__ == "__main__": main()