""" ============================================================================= 第 12 讲 · 代码 2：有特征的老虎机 —— 个性化剂量（模拟的华法林问题）与逻辑斯蒂 Thompson 采样 ============================================================================= 对应讲义：lectures/lecture12_thompson_acquisition.md 第 4 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 华法林是一种常用的抗凝药，合适的剂量因人而异：太少没效果，太多会出血。合适的剂量和年龄、体重、 两个基因（VKORC1、CYP2C9）、是否同时服用胺碘酮等有关（国际华法林药物基因组学联盟 IWPC 2009 给出了剂量公式）。 Bastani & Bayati（2020）把它做成一个老虎机问题：病人一个一个来，医生（算法）从"低 / 中 / 高"三档剂量里选一档， 选对了奖励 1，选错了 0。 **本文件用的是模拟数据**（仿照上面的设定：每周剂量 < 21 mg 为低、21~49 为中、> 49 为高），不是真实病人数据。 算法（每个剂量一个逻辑斯蒂回归 P(选对 | 病人特征)，用拉普拉斯近似得到后验，每 10 个病人重新拟合一次）： · 固定中剂量（不看特征；临床上常见的起始剂量） · 贪心：选"预测选对概率"最大的剂量 · ε-贪心：10% 的病人随机选 · UCB：选 logit 的"均值 + 1 倍标准差"最大的剂量 · Thompson 采样：从每个剂量的后验各抽一组权重，选抽到的 logit 最大的剂量 · 最优的个性化策略（知道真实的剂量模型）：选对的概率的上限 演示 1：整体的选对率随病人数的变化；相对最优策略的累计遗憾。 演示 2：需要高剂量的少数病人：各算法给他们高剂量的比例。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 02_logistic_thompson.py # 约 15 秒 只需要 numpy。依赖本目录的 ts_tools.py。 """ import numpy as np from ts_tools import LogisticArm, norm_cdf, pad, print_header DOSES = ["低", "中", "高"] NOISE = 7.0 # 同样特征的病人，真实需要的剂量的个体差异（标准差，mg/周） ALGOS = {"fixed": "固定中剂量", "greedy": "贪心", "eps": "ε-贪心（ε=0.1）", "ucb": "UCB", "ts": "Thompson 采样", "oracle": "最优的个性化策略"} def make_patients(n, rng): """ 模拟 n 个病人。 【返回】(X, need, best)： X ：(n, 6) 特征（常数项、年龄、体重、VKORC1 变异个数、CYP2C9 变异、是否服用胺碘酮，已中心化） need：每个病人真正需要的剂量档（0 低 / 1 中 / 2 高） best：知道真实模型时最好的猜测（"最优的个性化策略"会选的档） 【模型】需要的剂量（mg/周）= 55 − 2.2·年龄（十年为单位）+ 0.4·(体重 − 75) − 10·VKORC1 − 9·CYP2C9 − 9·胺碘酮 + 噪声 噪声 ~ N(0, 7²)：同样特征的人也会不同，所以即使知道模型，也不可能总是选对。 """ age = rng.integers(2, 10, n) weight = rng.normal(75, 15, n) vk = rng.choice(3, n, p=[0.35, 0.45, 0.20]) cyp = (rng.random(n) < 0.2).astype(float) amio = (rng.random(n) < 0.08).astype(float) mu = 55 - 2.2 * age + 0.4 * (weight - 75) - 10.0 * vk - 9.0 * cyp - 9.0 * amio dose = mu + rng.normal(0, NOISE, n) need = np.where(dose < 21, 0, np.where(dose <= 49, 1, 2)) p_low = norm_cdf((21 - mu) / NOISE) p_high = 1 - norm_cdf((49 - mu) / NOISE) best = np.argmax(np.stack([p_low, 1 - p_low - p_high, p_high], axis=1), axis=1) X = np.stack([np.ones(n), (age - 5.5) / 2.3, (weight - 75) / 15, vk - 0.85, cyp - 0.2, amio - 0.08], axis=1) return X, need, best def run(algo, rep, N=5000, refit=10): """ 给 N 个病人依次选剂量。 【返回】(correct, chose, best)：每个病人是否选对、选了哪一档、最优策略会选哪一档 """ rng = np.random.default_rng(rep) X, need, best = make_patients(N, rng) arms = [LogisticArm(X.shape[1]) for _ in DOSES] correct = np.zeros(N) chose = np.zeros(N, dtype=int) for t in range(N): x = X[t:t + 1] if algo == "fixed": a = 1 elif algo == "oracle": a = best[t] else: if algo == "ts": score = [arm.sample_logit(x, rng)[0] for arm in arms] elif algo == "ucb": score = [arm.mean_logit(x)[0] + arm.std_logit(x)[0] for arm in arms] elif algo == "eps" and rng.random() < 0.1: score = list(rng.random(3)) else: score = [arm.mean_logit(x)[0] for arm in arms] a = int(np.argmax(np.array(score) + 1e-9 * rng.random(3))) correct[t] = float(a == need[t]) chose[t] = a if algo not in ("fixed", "oracle"): arms[a].add(x[0], correct[t]) if (t + 1) % refit == 0: for arm in arms: arm.refit() return correct, chose, best def main(reps=10): rng = np.random.default_rng(0) _, need, best = make_patients(50000, rng) print_header("演示 1：整体的选对率") print("模拟病人（不是真实数据）。需要低 / 中 / 高剂量的比例：" + " / ".join(f"{p:.2f}" for p in np.bincount(need) / len(need)) + f"；“最优的个性化策略”最多能选对 {np.mean(best == need):.3f}（个体差异让它不可能总对）。") print(f"每种算法 {reps} 次独立实验，每次 5000 个病人：\n") windows = [(0, 200), (200, 1000), (1000, 3000), (3000, 5000)] print(f" {pad('算法', 18)}" + "".join(pad(f"第 {a + 1}~{b} 个", 16, "right") for a, b in windows) ) results = {} for algo, name in ALGOS.items(): res = [run(algo, rep) for rep in range(reps)] results[algo] = res C = np.array([r[0] for r in res]) # (实验数, 病人数)：每个病人是否选对 results[algo + "_C"] = C print(f" {pad(name, 18)}" + "".join(f"{C[:, a:b].mean():>16.3f}" for a, b in windows)) oracle_total = results["oracle_C"].sum(axis=1).mean() print("\n 累计遗憾（相对最优的个性化策略少选对了多少个病人，5000 个病人里）：") for algo, name in ALGOS.items(): if algo == "oracle": continue print(f" {pad(name, 18)} {oracle_total - results[algo + '_C'].sum(axis=1).mean():>8.1f}") print(""" 观察： · 固定中剂量对大部分病人（需要中剂量的约 69%）是对的，但对需要低或高剂量的病人总是错的。 · 几种学习算法都很快超过了固定剂量。出人意料的是，**贪心的整体选对率和 UCB、Thompson 采样相差不大**： 病人的特征千差万别，每来一个新病人，贪心的“最好猜测”自然就落在了不同的地方 —— 特征的多样性 本身提供了探索（Bastani、Bayati、Khosravi 2021：“几乎不需要探索”的情境老虎机）。 · 累计遗憾 UCB 最小；Thompson 采样一开始更“爱冒险”（前 200 个病人选对率最低），之后追上来， 累计遗憾和贪心差不多。在这个问题上，“探索”的好处没有第 11 讲那么明显 —— 但看演示 2。 · 所有学习算法都停在最优策略以下：“选中剂量对不对”和特征的关系不是单调的（中剂量只在“不高不低”时对）， 逻辑斯蒂回归的线性 logit 表示不了 —— 模型设定错误给学习设了上限（第 8 讲的教训在这里再次出现）。""") return results def demo_minority(results): print_header("演示 2：需要高剂量的少数病人") print("只看第 3001~5000 个病人里“最优策略会给高剂量”的那些（约 7%），各算法给了他们哪一档：\n") print(f" {pad('算法', 18)} {pad('给了低', 8, 'right')} {pad('给了中', 8, 'right')} {pad('给了高', 8, 'right')}") for algo, name in ALGOS.items(): frac = np.zeros(3) for correct, chose, best in results[algo]: m = best[3000:] == 2 frac += np.bincount(chose[3000:][m], minlength=3) / m.sum() frac /= len(results[algo]) print(f" {pad(name, 18)} {frac[0]:>8.2f} {frac[1]:>8.2f} {frac[2]:>8.2f}") print(""" 观察： · 贪心几乎从不给高剂量：一开始“高剂量”的模型数据很少、预测很悲观，贪心就不选它；不选它，就永远学不到 “哪些病人需要高剂量” —— 第 11 讲“贪心卡住”的情境版本。整体选对率看不出这个问题，因为这些病人只占 7%。 · UCB 和 Thompson 采样对“不确定”的剂量给了机会，给高剂量的比例明显更高，但 5000 个病人之后仍远没学会。 · **平均指标会掩盖对少数群体的系统性失败**。在医疗、推荐、内容审核里，这种失败恰恰最需要关注（第 17 讲的公平）。 评估一个在线学习系统，要按子群体看，而不只看平均。""") if __name__ == "__main__": demo_minority(main())