""" ============================================================================= 第 11 讲 · 共享模块：多臂老虎机、并行模拟器与各种选臂策略 ============================================================================= 对应讲义：lectures/lecture11_bandits_ucb.md ----------------------------------------------------------------------------- 【这个模块提供什么？】 ----------------------------------------------------------------------------- 1. 伯努利多臂老虎机：K 个"臂"（老虎机的拉杆），拉第 a 个臂，以概率 μ_a 得到奖励 1，否则得到 0。 μ_a 对智能体是**未知**的——它只能通过拉臂、看奖励来慢慢了解。 2. 一个**并行模拟器** simulate()：同时独立地跑 n_runs 次实验（相当于 n_runs 个随机种子）， 记录每次实验的累积遗憾（regret）。 3. 各种选臂策略：贪心、ε-贪心、递减 ε-贪心、乐观初始化、UCB1、KL-UCB、先探索后承诺（ETC）。 ----------------------------------------------------------------------------- 【为什么要"并行"模拟？】 ----------------------------------------------------------------------------- 老虎机实验要跑很多步（1 万～10 万步），还要重复很多次（几百次）才能得到可靠的平均。 如果用 Python 的 for 循环一次一次地跑，会非常慢。 这里的技巧：把 n_runs 次实验"叠在一起"，用形状为 (n_runs, K) 的数组同时保存所有实验的状态， 每一步用 numpy 一次性为所有实验选臂、发奖励。for 循环只剩下"时间步"这一层。 ----------------------------------------------------------------------------- 【策略的统一接口】 ----------------------------------------------------------------------------- 每个策略都是一个函数 policy(counts, sums, t, rng) → actions： counts : 形状 (n_runs, K)，counts[i, a] = 第 i 次实验中臂 a 被拉过的次数 N_a sums : 形状 (n_runs, K)，sums[i, a] = 第 i 次实验中臂 a 获得的奖励总和 t : 当前是第几步（从 1 开始） rng : 随机数生成器 返回 : 形状 (n_runs,) 的整数数组，第 i 个元素是第 i 次实验这一步要拉的臂 所有策略只根据"拉了几次、得了多少奖励"做决定——这正是老虎机问题中智能体能看到的全部信息。 为了统一接口，每个策略都写成"工厂函数"：先调用 greedy() 得到一个策略函数，再交给模拟器。 """ import unicodedata import numpy as np # 本讲默认使用的老虎机：5 个臂，最好的是臂 4（μ = 0.7） DEFAULT_MEANS = np.array([0.2, 0.4, 0.5, 0.6, 0.7]) # ============================================================================= # 第 1 部分：辅助函数 # ============================================================================= def argmax_random_tiebreak(values, rng): """ 对每一行取最大值的下标；如果有好几个并列最大，**随机**选一个。 【参数】values : 形状 (n_runs, K) 的数组（可以包含 np.inf） 【返回】形状 (n_runs,) 的整数数组 【为什么要随机打破平局？】伯努利奖励只有 0 和 1，经验均值经常出现并列（比如两个臂都是 1/2）。 np.argmax 遇到并列总是返回**第一个**，这会让策略系统性地偏向编号小的臂—— 一个隐蔽的 bug，会让实验结果依赖于"臂的编号顺序"。 【做法】给每个并列最大的位置一个随机数，其余位置设为 -1，再取 argmax。 """ is_max = values >= values.max(axis=1, keepdims=True) noise = rng.random(values.shape) return np.argmax(np.where(is_max, noise, -1.0), axis=1) def pad(text, width): """ 在字符串右边补空格，使它在终端里的**显示宽度**为 width。用来对齐含中文的表格。 【为什么需要？】Python 的 f"{s:<10s}" 按"字符个数"补齐，但一个中文字符在终端里占 2 个英文字符的宽度， 所以含中文的行会错位。unicodedata.east_asian_width(ch) 返回 "W"（宽）或 "F"（全角）时，该字符占 2 格。 【边界】text 已经比 width 宽时，原样返回（不截断）。 """ w = sum(2 if unicodedata.east_asian_width(ch) in "WF" else 1 for ch in text) return text + " " * max(0, width - w) def empirical_means(counts, sums): """ 经验均值 μ̂_a = sums / counts。没拉过的臂（counts = 0）返回 0。 【Python 小知识】np.divide(a, b, out=..., where=...) 只在 where 为真的位置做除法， 其余位置保留 out 中的默认值——避免"除以 0"的警告。 """ return np.divide(sums, counts, out=np.zeros_like(sums), where=counts > 0) def bernoulli_kl(p, q): """ 两个伯努利分布之间的 KL 散度：kl(p, q) = p log(p/q) + (1-p) log((1-p)/(1-q))。 【直觉】"真实成功率是 p，却以为是 q"时，每观察一次平均能获得多少"区分它们"的信息。 p 和 q 越接近，kl 越小，要区分它们需要的样本就越多（约 1/kl 次）。 【边界】p 或 q 等于 0 或 1 时 log 会出问题，所以先把它们裁剪到 [1e-12, 1-1e-12]。 【RL 联系】Lai–Robbins 遗憾下界（03 文件）和 KL-UCB 算法都用到它。 """ p = np.clip(p, 1e-12, 1 - 1e-12) q = np.clip(q, 1e-12, 1 - 1e-12) return p * np.log(p / q) + (1 - p) * np.log((1 - p) / (1 - q)) # ============================================================================= # 第 2 部分：并行模拟器 # ============================================================================= def simulate(policy, means, T, n_runs, seed, checkpoints=None): """ 同时独立地跑 n_runs 次老虎机实验，每次 T 步。 【参数】 policy : 选臂策略（见文件开头的"统一接口"） means : 每个臂的真实成功率 μ_a（长度 K） T : 每次实验的总步数 n_runs : 独立重复的次数（相当于 n_runs 个随机种子） seed : 随机种子（同样的 seed 得到同样的结果，方便复现） checkpoints : 在哪些时刻记录累积遗憾，例如 [100, 1000, 10000]；默认只记录第 T 步 【返回】一个字典： "checkpoints" : 记录的时刻 "regret" : 形状 (n_runs, 记录时刻数)，每次实验在每个记录时刻的累积遗憾 "counts" : 形状 (n_runs, K)，实验结束时每个臂被拉的次数 "first_pulls" : 形状 (n_runs, 前 20 步)，每次实验前 20 步拉的臂（用来观察策略的行为） 【遗憾的定义】这里记录的是**伪遗憾（pseudo-regret）**： Regret(T) = Σ_{t=1}^{T} (μ* - μ_{a_t}) = Σ_a Δ_a · N_a(T) 其中 μ* 是最好的臂的成功率，Δ_a = μ* - μ_a 是臂 a 的"差距"。 它衡量"因为没有一直拉最好的臂，期望上少拿了多少奖励"。 用期望奖励 μ 而不是实际拿到的 0/1 奖励来算，可以去掉奖励本身的随机噪声，曲线更平滑。 【边界】checkpoints 中大于 T 的时刻会被忽略。 """ means = np.asarray(means, dtype=float) K = len(means) rng = np.random.default_rng(seed) counts = np.zeros((n_runs, K)) sums = np.zeros((n_runs, K)) gap = means.max() - means # 每个臂的差距 Δ_a cum_regret = np.zeros(n_runs) checkpoints = sorted(c for c in set(checkpoints or [T]) if c <= T) record = np.zeros((n_runs, len(checkpoints))) first_pulls = np.zeros((n_runs, min(20, T)), dtype=int) rows = np.arange(n_runs) # 0, 1, ..., n_runs-1：用来同时索引每一次实验 ci = 0 for t in range(1, T + 1): a = policy(counts, sums, t, rng) # 每次实验各选一个臂 r = (rng.random(n_runs) < means[a]).astype(float) # 以概率 μ_a 得到 1 counts[rows, a] += 1 sums[rows, a] += r # 【Python 小知识】counts[rows, a] 是"花式索引"：同时取出 (0, a[0]), (1, a[1]), ... 这些位置。 # 因为每一行只有一个位置被选中，不会有重复下标，所以直接 += 是安全的。 cum_regret += gap[a] if t <= first_pulls.shape[1]: first_pulls[:, t - 1] = a if ci < len(checkpoints) and t == checkpoints[ci]: record[:, ci] = cum_regret ci += 1 return {"checkpoints": np.array(checkpoints), "regret": record, "counts": counts, "first_pulls": first_pulls} def lai_robbins_constant(means): """ Lai–Robbins 下界的常数：C = Σ_{a: Δ_a > 0} Δ_a / kl(μ_a, μ*)。 【含义】对任何"在所有伯努利老虎机上都表现合理"的算法，当 T → ∞ 时 Regret(T) ≥ (C - o(1)) · ln T 也就是说：遗憾**至少**要按 ln T 增长，而且增长的系数不可能小于 C。详见 03 文件。 """ means = np.asarray(means, dtype=float) best = means.max() sub = means < best return float(np.sum((best - means[sub]) / bernoulli_kl(means[sub], best))) # ============================================================================= # 第 3 部分：选臂策略 # ============================================================================= def greedy(): """ 贪心：先把每个臂各拉一次，之后永远拉"经验均值最高"的臂。 【实现】没拉过的臂的值设为 +∞，保证它们先被拉到（这是"每个臂先拉一次"的常用写法）。 【问题】完全不探索。如果最好的臂一开始运气不好（比如第一次拉得到 0），它的经验均值是 0， 可能永远不会再被拉——贪心会**永远卡在**一个次优的臂上，遗憾随 T **线性**增长。 """ def policy(counts, sums, t, rng): values = np.where(counts > 0, empirical_means(counts, sums), np.inf) return argmax_random_tiebreak(values, rng) return policy def epsilon_greedy(eps): """ ε-贪心：以概率 ε 随机拉一个臂（探索），否则拉经验均值最高的臂（利用）。 【参数】eps : 探索概率，0 ≤ eps ≤ 1。eps = 0 就是贪心；eps = 1 就是完全随机。 【问题】探索的比例**永远**是 ε——即使已经非常确定哪个臂最好，每一步仍有 ε 的概率随机拉， 而随机拉到的大多是次优的臂。所以遗憾至少是 ε · (平均差距) · T，**线性**增长。 """ def policy(counts, sums, t, rng): n_runs, K = counts.shape values = np.where(counts > 0, empirical_means(counts, sums), np.inf) greedy_a = argmax_random_tiebreak(values, rng) random_a = rng.integers(K, size=n_runs) explore = rng.random(n_runs) < eps return np.where(explore, random_a, greedy_a) return policy def decaying_epsilon_greedy(c, d): """ 递减 ε-贪心（Auer, Cesa-Bianchi & Fischer, 2002）：ε_t = min(1, c·K / (d² · t))。 【参数】 c : 常数（越大探索越多） d : 对"最小差距 Δ_min"的猜测。理论保证要求 0 < d ≤ Δ_min。 【直觉】探索的概率按 1/t 递减，总探索次数约为 Σ_t c·K/(d²·t) ≈ (c·K/d²) · ln T——对数增长。 【问题】需要事先知道（或猜对）最小差距 Δ_min。d 猜得太大（以为臂之间差别很大）→ 探索不够， 可能卡在次优臂上；d 猜得太小 → 探索太多，浪费。实际问题中我们恰恰不知道 Δ_min。 """ def policy(counts, sums, t, rng): n_runs, K = counts.shape eps = min(1.0, c * K / (d * d * t)) values = np.where(counts > 0, empirical_means(counts, sums), np.inf) greedy_a = argmax_random_tiebreak(values, rng) random_a = rng.integers(K, size=n_runs) return np.where(rng.random(n_runs) < eps, random_a, greedy_a) return policy def optimistic_greedy(init_value=1.0, pseudo_count=1.0): """ 乐观初始化 + 贪心：假装每个臂在开始之前已经被拉了 pseudo_count 次、每次都得到 init_value。 【估计】Q_a = (sums_a + init_value · pseudo_count) / (counts_a + pseudo_count) 【直觉】所有臂一开始都"看起来很好"（奖励最多为 1，init_value = 1 是最乐观的）。 贪心地拉一个臂，它的真实表现不如预期，估计值下降，于是换别的臂试试—— 乐观的初始值自动驱动了探索。第 4 讲的练习中我们在 Q-learning 里见过同样的技巧。 【问题】乐观只在开始时存在，随着拉的次数增多被"冲淡"，之后退化为贪心。 如果最好的臂在乐观被冲淡之前运气不好，仍然可能被永远放弃。pseudo_count 越大，乐观持续越久。 """ def policy(counts, sums, t, rng): values = (sums + init_value * pseudo_count) / (counts + pseudo_count) return argmax_random_tiebreak(values, rng) return policy def ucb1(c=2.0): """ UCB1（Auer, Cesa-Bianchi & Fischer, 2002）：拉"置信上界"最高的臂。 UCB_a(t) = μ̂_a + sqrt( c · ln t / N_a ) 【直觉：面对不确定性时保持乐观】 第二项是"探索奖励"：拉得越少（N_a 小），这一项越大，这个臂"可能有多好"的上限就越高。 UCB 的做法是：**假设每个臂都和它"可能的最好情况"一样好**，然后选最好的。 - 如果一个臂的上界高是因为均值高 → 拉它是"利用"，很好； - 如果上界高是因为不确定 → 拉它是"探索"，拉完之后不确定性变小，上界降下来。 两种情况都不亏：要么拿到高奖励，要么学到东西。 【为什么是 sqrt(ln t / N_a)？】来自霍夫丁不等式（02 文件）： 以很大的概率（约 1 - 1/t^{c/2}），真实的 μ_a 落在 μ̂_a ± sqrt(c · ln t / (2 N_a)) 之内。 c = 2 是 UCB1 原论文中的选择（对应 sqrt(2 ln t / N_a)）。 ln t 缓慢增长保证了：即使一个臂看起来很差，它的上界也会慢慢升高，最终被再试一次—— 所以永远不会"永久放弃"任何一个臂。 【参数】c : 探索强度。c 越大探索越多。理论分析用 c = 2；实践中更小的 c 往往更好（02 演示 5）。 【边界】没拉过的臂 N_a = 0，上界为 +∞，保证每个臂先被拉一次。 """ def policy(counts, sums, t, rng): mean = empirical_means(counts, sums) bonus = np.sqrt(c * np.log(t) / np.maximum(counts, 1)) values = np.where(counts > 0, mean + bonus, np.inf) return argmax_random_tiebreak(values, rng) return policy def kl_ucb(n_bisect=16): """ KL-UCB（Garivier & Cappé, 2011）：用 KL 散度而不是霍夫丁不等式来构造置信上界。 UCB_a(t) = max { q ∈ [μ̂_a, 1] : N_a · kl(μ̂_a, q) ≤ ln t } 【直觉】霍夫丁不等式对所有 [0,1] 上的分布都成立，所以对具体的伯努利分布偏保守： 比如 μ̂_a = 0.05 时，真实成功率不太可能是 0.5，但霍夫丁给出的上界却可能超过 0.5。 KL 散度"知道"伯努利分布的形状，给出更紧的上界。 它的遗憾恰好达到 Lai–Robbins 下界的常数——渐近最优（03 文件）。 【实现】kl(μ̂, q) 随 q 在 [μ̂, 1] 上单调递增，所以可以用**二分法**找满足条件的最大 q。 n_bisect 次二分后精度约为 2^{-n_bisect}。所有实验、所有臂同时做二分（向量化）。 """ def policy(counts, sums, t, rng): mean = empirical_means(counts, sums) n = np.maximum(counts, 1) lo, hi = mean.copy(), np.ones_like(mean) target = np.log(t) / n for _ in range(n_bisect): mid = (lo + hi) / 2 ok = bernoulli_kl(mean, mid) <= target # mid 仍在置信范围内 → 上界至少是 mid lo = np.where(ok, mid, lo) hi = np.where(ok, hi, mid) values = np.where(counts > 0, lo, np.inf) return argmax_random_tiebreak(values, rng) return policy def explore_then_commit(m): """ 先探索后承诺（Explore-Then-Commit, ETC）：先把每个臂各拉 m 次，然后**只根据这段探索的结果** 选出经验均值最高的臂，此后永远拉它，不再改变主意。 【参数】m : 每个臂的探索次数。 【直觉】最朴素的"先做实验，再做决定"——就像 A/B 测试：先各放 m 个用户看效果，再全量上线好的版本。 【问题】最佳的 m 取决于差距 Δ 和总步数 T（大约 m ∝ ln(T Δ²) / Δ²），而 Δ 恰恰是未知的。 m 太小 → 可能选错，之后每一步都亏 Δ（线性遗憾）；m 太大 → 在差的臂上浪费太多。 【实现上的一个坑】如果承诺之后每一步都重新选"经验均值最高的臂"，那么承诺的臂的均值会继续更新， 选错了还能被纠正过来——那就变成了"强制探索 + 贪心"，不再是 ETC。 这里的做法：探索刚结束的那一步（所有臂都恰好拉了 m 次）按经验均值做决定； 之后被承诺的臂拉的次数一定最多，所以直接选"拉得最多的臂"，就能一直坚持最初的决定。 """ def policy(counts, sums, t, rng): still_exploring = counts.min(axis=1) < m explore_a = argmax_random_tiebreak(-counts, rng) # 拉得最少的臂（轮流拉） deciding_now = counts.max(axis=1) == m # 探索刚结束：所有臂都恰好 m 次 decide_a = argmax_random_tiebreak(empirical_means(counts, sums), rng) committed_a = argmax_random_tiebreak(counts, rng) # 之前承诺的臂 = 拉得最多的臂 commit_a = np.where(deciding_now, decide_a, committed_a) return np.where(still_exploring, explore_a, commit_a) return policy if __name__ == "__main__": # 直接运行本模块：做一个小的自检 res = simulate(ucb1(), DEFAULT_MEANS, T=1000, n_runs=100, seed=0, checkpoints=[100, 1000]) print("默认老虎机的真实成功率：", DEFAULT_MEANS) print("UCB1 在 t = 100、1000 时的平均累积遗憾：", res["regret"].mean(axis=0).round(2)) print("Lai–Robbins 常数 C =", round(lai_robbins_constant(DEFAULT_MEANS), 2))