""" ============================================================================= 第 4 讲 · 代码 1：Rasch 模型 —— 用户的"胃口"、选项的"吸引力"，以及它们怎样在比较中分离 ============================================================================= 对应讲义：lectures/lecture04_factor_models.md 第 2 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- Rasch（1960）为考试设计了一个模型：学生 i 答对题目 j 的概率 = σ(能力_i − 难度_j)。 换到偏好数据：用户 i 接受选项 j 的概率 = σ(a_i + z_j)，a_i 是用户的"胃口"（有人什么都点赞， 有人很挑剔），z_j 是选项的"吸引力"。这是最简单的"用户 × 选项"模型。 演示 1：从接受 / 拒绝数据里同时估计 300 个用户的胃口和 20 个回答的吸引力。 演示 2：Rasch 的"特定客观性"—— 只看"两个选项恰好接受了一个"的用户，选的是哪一个的概率 等于 σ(z_j − z_k)，和用户的胃口**无关**。这就是第 1 讲"打分的宽严在比较中抵消"的严格形式， 也是"成对比较只能学到选项之间的差"的原因。 演示 3：一个反直觉的统计陷阱（"附带参数问题"）—— 每个用户只答 2 道题时，同时估计用户和题目的 "联合最大似然"会把题目之间的差**估成真实值的 2 倍**；只用演示 2 的"条件"比较就没有这个问题。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 01_rasch.py # 约 5 秒 只用到 Python 标准库，以及第 1 讲的 toy_world.py、本目录的 factor_models.py。 """ import random from factor_models import fit_mf, predict_prob from toy_world import fit_bradley_terry, mean, pad, print_header, sigmoid def simulate_rasch(users, appetite, appeal, rng, items=None): """ 按 Rasch 模型模拟接受 / 拒绝：P(Y_ij = 1) = σ(a_i + z_j)。 【参数】users：用户列表；appetite：dict {用户: a_i}；appeal：dict {选项: z_j}；rng； items：每个用户回答哪些选项（默认全部） 【返回】list[(用户, 选项, y)] """ items = items or sorted(appeal) obs = [] for i in users: for j in items: y = 1 if rng.random() < sigmoid(appetite[i] + appeal[j]) else 0 obs.append((i, j, y)) return obs def correlation(xs, ys): """皮尔逊相关系数（−1 ~ 1）。""" mx, my = mean(xs), mean(ys) sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) sxx = sum((x - mx) ** 2 for x in xs) syy = sum((y - my) ** 2 for y in ys) return sxy / (sxx * syy) ** 0.5 def slope(est, true): """ 把"估计值"对"真实值"做最小二乘回归（两者都先减去平均值），返回斜率。 【思路】斜率 = 1：尺度估对了；> 1：估计值的差距被放大；< 1：被压缩。 """ me, mt = mean(est), mean(true) return sum((e - me) * (t - mt) for e, t in zip(est, true)) / sum((t - mt) ** 2 for t in true) # ============================================================================= # 演示 1：同时估计胃口和吸引力 # ============================================================================= def demo_fit(rng): print_header("演示 1：从接受 / 拒绝数据里，同时估计用户的胃口和回答的吸引力") users = [f"u{i:03d}" for i in range(300)] items = [f"r{j:02d}" for j in range(20)] appetite = {i: rng.gauss(0, 1.0) for i in users} appeal = {j: rng.gauss(0, 1.0) for j in items} obs = simulate_rasch(users, appetite, appeal, rng) model = fit_mf(obs, users, items, k=0, l2=0.1, rng=rng) # 只有 a + z 可识别：所有 a 加 c、所有 z 减 c，预测不变。比较之前把两边都平移到平均值为 0。 print(f"300 个用户，每人对 20 个回答各给一次接受 / 拒绝（共 {len(obs)} 条）。Rasch 模型 P(接受) = σ(a_i + z_j)。") print(f" 回答的吸引力：估计值和真实值的相关系数 {correlation([model.z[j] for j in items], [appeal[j] for j in items]):.3f}，" f"斜率 {slope([model.z[j] for j in items], [appeal[j] for j in items]):.3f}") print(f" 用户的胃口： 估计值和真实值的相关系数 {correlation([model.a[i] for i in users], [appetite[i] for i in users]):.3f}，" f"斜率 {slope([model.a[i] for i in users], [appetite[i] for i in users]):.3f}") mean_a = mean([model.a[i] for i in users]) print(f" 估计的胃口平均值 {mean_a:+.3f}、吸引力平均值 {mean([model.z[j] for j in items]):+.3f}" "（只有 a + z 可识别，平均值本身没有意义）") print(""" 观察： · 回答的吸引力估得很准（每个回答有 300 个用户的数据）；用户的胃口估得粗一些（每个用户只有 20 条数据）。 这是 Rasch 模型（以及所有因子模型）的普遍现象：数据多的一方估得准，数据少的一方估得粗。 · 胃口的斜率略小于 1：L2 正则化（l2 = 0.1）把每个用户的胃口往 0 拉，数据越少拉得越明显（“收缩”）。 不加正则化时方向正好相反：每人只有很少的数据，最大似然会把胃口和吸引力都估得**更极端** —— 演示 3。 · Rasch 模型在教育测量里叫“项目反应理论”（IRT）的单参数模型，是计算机自适应考试的基础（第 7 讲）。""") # ============================================================================= # 演示 2：特定客观性 —— 用户胃口在比较中抵消 # ============================================================================= def demo_objectivity(rng): print_header("演示 2：“恰好接受了一个”的时候，选哪个和胃口无关") users = [f"u{i:04d}" for i in range(6000)] appetite = {i: rng.gauss(0, 1.5) for i in users} appeal = {"甲": 0.8, "乙": -0.4} obs = simulate_rasch(users, appetite, appeal, rng) got = {} for i, j, y in obs: got.setdefault(i, {})[j] = y groups = sorted(users, key=lambda i: appetite[i]) third = len(groups) // 3 print("6000 个用户（胃口 ~ N(0, 1.5²)）各看两个回答：甲（吸引力 0.8）、乙（吸引力 −0.4）。") print("按胃口把用户分成三组，看“甲、乙恰好接受了一个”的用户里，接受的是甲的比例：\n") print(f" {pad('用户组', 14)} {pad('接受率（甲）', 12, 'right')} {pad('恰好接受一个的人数', 18, 'right')} {pad('其中接受甲的比例', 16, 'right')}") for name, grp in [("胃口最低的 1/3", groups[:third]), ("中间的 1/3", groups[third:2 * third]), ("胃口最高的 1/3", groups[2 * third:])]: rate = mean([got[i]["甲"] for i in grp]) one = [i for i in grp if got[i]["甲"] + got[i]["乙"] == 1] frac = mean([got[i]["甲"] for i in one]) print(f" {pad(name, 14)} {rate:>12.3f} {len(one):>18d} {frac:>16.3f}") print(f" 理论值 σ(0.8 − (−0.4)) = {sigmoid(1.2):.3f}") print(""" 观察： · 三组人的接受率差得很远（挑剔的人很少接受），但在“恰好接受了一个”的人里，接受甲的比例都约等于 σ(1.2)。 推导一行：P(接受甲、拒绝乙) / P(恰好接受一个) = σ(a+z甲)(1−σ(a+z乙)) / [同上 + (1−σ(a+z甲))σ(a+z乙)] = e^{z甲} / (e^{z甲} + e^{z乙}) = σ(z甲 − z乙)，a 完全消掉了。 · Rasch 把这叫“特定客观性”（specific objectivity）：比较两道题的难度，不需要知道是谁在答题； 比较两个人的能力，也不需要知道是哪道题。 · 对偏好学习的含义（教材第 1.10 节）： 逐项数据（接受 / 拒绝）能估计用户的胃口 a_i 和选项的吸引力 z_j； 两两比较只能估计选项之间的差 z_j − z_k —— 用户的胃口在比较里消掉了。 第 1 讲代码 01 的“打分的宽严在比较中抵消”就是这件事：宽严 = 胃口。""") # ============================================================================= # 演示 3：附带参数问题 # ============================================================================= def demo_incidental(rng, n_users=3000): print_header("演示 3：每人只答几道题时，“联合最大似然”会把差距估大") print(f"""{n_users} 个用户（胃口 ~ N(0, 1)），每人回答同样的 n 个选项，吸引力在 [−1.5, 1.5] 里等距分布。 两种估计选项吸引力的方法： · 联合最大似然（JMLE）：同时估计所有用户的胃口和所有选项的吸引力（去掉“全接受 / 全拒绝”的用户， 他们的胃口是 ±∞，不提供选项的信息）。 · 条件最大似然（CMLE 的简化版）：只用演示 2 的比较 —— 对每个用户的每一对选项，如果恰好接受了一个， 记一次“接受的 ≻ 拒绝的”，再用 BT 拟合。用户的胃口根本不出现。 报告“估计的吸引力”对“真实吸引力”的斜率（1 表示尺度正确）： """) print(f" {pad('每人回答几个', 12)} {pad('JMLE 的斜率', 12, 'right')} {pad('条件比较的斜率', 14, 'right')} " f"{pad('理论值 n/(n−1)', 14, 'right')}") for n in [2, 5, 20]: items = [f"r{j:02d}" for j in range(n)] appeal = {j: -1.5 + 3.0 * t / (n - 1) for t, j in enumerate(items)} users = [f"u{i:04d}" for i in range(n_users)] appetite = {i: rng.gauss(0, 1.0) for i in users} obs = simulate_rasch(users, appetite, appeal, rng) score = {} for i, j, y in obs: score[i] = score.get(i, 0) + y keep = [i for i in users if 0 < score[i] < n] kept_obs = [(i, j, y) for i, j, y in obs if i in set(keep)] jmle = fit_mf(kept_obs, keep, items, k=0, l2=1e-6, rng=rng, sweeps=30) got = {} for i, j, y in obs: got.setdefault(i, {})[j] = y comps = [] for i in users: for x in range(n): for w in range(x + 1, n): a, b = items[x], items[w] if got[i][a] + got[i][b] == 1: comps.append((a, b) if got[i][a] == 1 else (b, a)) cml = fit_bradley_terry(comps, items, l2=0.0, tol=1e-10) true = [appeal[j] for j in items] s_j = slope([jmle.z[j] for j in items], true) s_c = slope([cml[j] for j in items], true) print(f" {n:>12d} {s_j:>12.3f} {s_c:>14.3f} {n / (n - 1):>14.3f}") print(""" 观察： · 每人只答 2 道题时，JMLE 把两道题的吸引力之差估成了真实值的**两倍**（Andersen 1973 证明了这个精确的结果）； 答 5 道题时大约高估 1/4；答 20 道题时只高估约 5%。经验修正因子 (n−1)/n 正是从这里来的。 · 原因：每个用户的胃口只由他自己的几条数据估计，永远估不准（数据量不随用户数增加）， 这个误差又“传染”给了选项的估计 —— 这叫“附带参数问题”（Neyman & Scott 1948）： 参数个数随数据量一起增长时，最大似然可能**不一致**（数据再多也不收敛到真值）。 · 条件比较完全不受影响：它根本不估计胃口，只用“恰好接受了一个”的信息。 · 对偏好学习的启示：给每个标注员（或每个用户）一个参数，而每人的数据又很少时，要小心这类偏差。 正则化（第 5 讲）和把用户参数当作随机效应积分掉（第 6 讲的贝叶斯视角）是常见的补救。""") def main(): rng = random.Random(1960) # 1960：Rasch 出版他的测验模型专著的年份 demo_fit(rng) demo_objectivity(rng) demo_incidental(rng) if __name__ == "__main__": main()