""" ============================================================================= 第 2 讲 · 01：Large Language Monkeys 的 power law 从哪里来？ ============================================================================= 对应讲义：lectures/lecture02_test_time_compute.md 第 3 节 ----------------------------------------------------------------------------- 【一个谜题】 ----------------------------------------------------------------------------- Large Language Monkeys 发现：对一组题目重复采样 k 次， −log(coverage@k) ≈ A · k^(−b) （power law，在 log-log 图上是直线） 但对**单独一道**单次成功率为 p 的题，k 次都失败的概率是 (1 − p)^k —— 这是**指数**衰减， 比任何 power law 都快得多。指数衰减的题目加在一起，怎么会变成 power law？ Schaeffer et al. 2025（How Do Large Language Monkeys Get Their Power (Laws)?）的回答： **看单次成功率 p 在题目之间的分布。** 如果有一小部分题的 p 非常接近 0（分布的"左尾"很重）， 它们需要极多的采样才能被解决，正是这些题把整体的曲线"拖"成了 power law。 ----------------------------------------------------------------------------- 【两个演示】 ----------------------------------------------------------------------------- 演示 1 纯数学：假设 p 服从 Beta(a, β) 分布（密度在 0 附近像 p^(a−1)）， 可以精确算出 1 − coverage@k = B(a, β + k) / B(a, β) ≈ 常数 × k^(−a)。 指数 b 就是 a：**左尾越重（a 越小），coverage 增长得越慢**。 如果没有任何题的 p 接近 0（例如所有 p ≥ 0.05），整体最终还是指数衰减，power law 消失。 演示 2 toy world：1B 模型，100 道题，每题采样 10000 次得到"真实"的 coverage 曲线。 只用每题前 100 个样本，比较两种预测 coverage@10000 的方法： · 直接拟合：在 k ≤ 100 上拟合 power law，然后外推； · 分布拟合：用 100 个样本估计 p 的分布（beta-binomial 最大似然），再用演示 1 的公式算。 运行时间：约 10 秒（演示 2 要生成一百万个样本；只用标准库）。 """ import math import random from tts_tools import ToyLLM, make_problem_set, mean, pad # ============================================================================= # 演示 1：Beta 分布下的精确公式 # ============================================================================= def log_beta(a, b): """ log B(a, b)，其中 B 是 Beta 函数：B(a, b) = Γ(a)Γ(b) / Γ(a + b)。 【思路】直接算 Γ 函数很容易溢出（Γ(171) 就超出 float 的范围了），所以在对数空间里算： log B(a, b) = log Γ(a) + log Γ(b) − log Γ(a + b) 【参数】a、b：正数 【返回】float 【Python 小知识】math.lgamma(x) 返回 log|Γ(x)|，对很大的 x 也不会溢出。 """ return math.lgamma(a) + math.lgamma(b) - math.lgamma(a + b) def beta_failure(a, b, k): """ p ~ Beta(a, b) 时，"k 次全失败"的平均概率 E[(1 − p)^k] = B(a, b + k) / B(a, b)。 【思路】(1 − p)^k 乘上 Beta 的密度 p^(a−1)(1−p)^(b−1)，正好变成 p^(a−1)(1−p)^(b+k−1)， 积分就是 B(a, b + k)。所以 1 − coverage@k 有精确的公式。 当 k 很大时，B(a, b + k) ≈ Γ(a) · k^(−a)，于是 1 − coverage@k ∝ k^(−a)：power law，指数是 a。 【参数】a、b：Beta 分布的参数（> 0）；k：采样次数（≥ 0） 【返回】0~1 的 float """ return math.exp(log_beta(a, b + k) - log_beta(a, b)) def uniform_failure(lo, k): """ p ~ Uniform[lo, 1] 时的 E[(1 − p)^k] = (1 − lo)^(k+1) / ((k + 1)(1 − lo))。 【思路】这是一个"没有任何题接近 0"的分布：最难的题也有 lo 的成功率。 积分 ∫_lo^1 (1−p)^k dp / (1 − lo) = (1 − lo)^(k+1) / ((k+1)(1 − lo))， 含有 (1 − lo)^k 这一项 —— k 大了以后是指数衰减，不是 power law。 【参数】lo：最小的成功率（0 < lo < 1）；k：采样次数 """ return (1 - lo) ** (k + 1) / ((k + 1) * (1 - lo)) def neg_log_coverage(failure): """ −log(coverage) = −log(1 − failure)。failure 很小时约等于 failure 本身。 【Python 小知识】math.log1p(x) = log(1 + x)。x 很小（例如 1e-20）时，直接算 1 + x 会被舍入成 1.0， log 的结果变成 0；log1p 专门为这种情况设计，能保留精度。 """ return -math.log1p(-failure) if failure < 1 else math.inf def demo_beta_math(): print("=" * 78) print("演示 1：每道题是指数衰减，整体却是 power law —— 关键在 p 的分布") print("=" * 78) ks = (1, 10, 100, 1000, 10000) print("(a) 单独一道题：−log(pass@k)，k 次至少一次成功的概率取负对数") print(pad("k", 8) + "".join(pad(f"p = {p}", 14, "right") for p in (0.3, 0.03, 0.003))) for k in ks: print(pad(str(k), 8) + "".join(pad(f"{neg_log_coverage((1 - p) ** k):.2e}", 14, "right") for p in (0.3, 0.03, 0.003))) print("\n(b) 一组题，p 服从不同的分布：−log(coverage@k)，括号里是相邻两行之间 log-log 的斜率") dists = (("Beta(0.3, 3)", lambda k: beta_failure(0.3, 3.0, k)), ("Beta(1, 3)", lambda k: beta_failure(1.0, 3.0, k)), ("Uniform[0.05,1]", lambda k: uniform_failure(0.05, k))) print(pad("k", 8) + "".join(pad(name, 24, "right") for name, _ in dists)) prev = None for k in ks: row = pad(str(k), 8) cur = [] for j, (_, f) in enumerate(dists): v = neg_log_coverage(f(k)) cur.append(v) if prev is None or v <= 0 or prev[j] <= 0: row += pad(f"{v:.2e}", 24, "right") else: slope = (math.log(v) - math.log(prev[j])) / (math.log(k) - math.log(ks[ks.index(k) - 1])) row += pad(f"{v:.2e} ({slope:+.2f})", 24, "right") prev = cur print(row) print("\n观察：") print(" · (a) 单独一道题：k 超过 1/p 以后，k 每扩大 10 倍，−log(pass@k) 就缩小很多个数量级（指数衰减）。") print(" p = 0.003 的题在 k = 100 时只有约 26% 的机会被做对（e^−1.35 ≈ 0.26），") print(" 到 k = 10000 时，失败的概率已经小到 1e−13。（0.00e+00 是数值太小，超出了 float 的范围。）") print(" · (b) Beta 分布：斜率趋近于 −a（Beta(0.3, 3) 约 −0.3，Beta(1, 3) 约 −1）—— 一条 power law。") print(" a 越小，0 附近的题越多，曲线越平：总有一批\"还没被碰到\"的难题在等着。") print(" · Uniform[0.05, 1]：起初也像 power law（斜率约 −1），但 k 一旦超过 1/0.05 = 20 的若干倍，") print(" 斜率就越来越陡，最后是指数衰减：没有\"左尾\"，就没有 power law。") print(" 这解释了论文里观察到的\"偏离 power law\"的现象：它取决于这组题里最难的题有多难。\n") # ============================================================================= # 演示 2：在 toy world 里预测 coverage # ============================================================================= def fail_at_k_stable(n, c, k): """ 1 − pass@k 的无偏估计（"随机抽 k 个样本全部失败"的概率）：Π_{i = n−c+1}^{n} (1 − k / i)。 【思路】为什么不直接算 1 − pass_at_k_stable(...)？当失败率极小（例如 1e−21）时，pass@k 在 float 里 会被舍入成 1.0，1 − 1.0 = 0，信息全丢了（叫"灾难性抵消"）。直接算连乘就没有这个问题。 【参数】同 pass_at_k_stable """ if n - c < k: return 0.0 prod = 1.0 for i in range(n - c + 1, n + 1): prod *= 1.0 - k / i return prod def pass_at_k_stable(n, c, k): """ pass@k 的无偏估计，用连乘写成数值稳定的形式（Chen et al. 2021 附录里的写法）： 1 − C(n−c, k)/C(n, k) = 1 − Π_{i = n−c+1}^{n} (1 − k / i) 【思路】和 toy_llm.pass_at_k 是同一个公式，但 n = 10000 时组合数非常大（几千位的整数）， 直接算组合数很慢。连乘只需要 c 次乘法，而且每一项都在 0~1 之间，不会溢出。 【参数】n ≥ 1；0 ≤ c ≤ n；1 ≤ k ≤ n 【返回】float """ return 1.0 - fail_at_k_stable(n, c, k) def fit_power_law_on_coverage(ks, coverages): """ 拟合 −log(coverage@k) = A · k^(−b)：两边取对数，log(−log cov) = log A − b · log k，是一条直线。 【参数】ks、coverages：等长的列表，coverage 必须在 (0, 1) 之间 【返回】predict(k) 函数和指数 b 【边界】coverage = 1 时 −log(cov) = 0，取不了对数 —— 调用者要保证传进来的点都 < 1。 """ xs = [math.log(k) for k in ks] ys = [math.log(-math.log(c)) for c in coverages] mx, my = mean(xs), mean(ys) slope = math.fsum((x - mx) * (y - my) for x, y in zip(xs, ys)) / math.fsum((x - mx) ** 2 for x in xs) intercept = my - slope * mx return (lambda k: math.exp(-math.exp(intercept + slope * math.log(k)))), -slope def fit_beta_binomial(counts, n): """ 用 beta-binomial 最大似然，从"每道题 n 次里成功了 c 次"估计 p 的分布 Beta(a, β)。 【思路】假设每道题的 p 来自 Beta(a, β)，那么"n 次里成功 c 次"的概率是 C(n, c) · B(c + a, n − c + β) / B(a, β) 把所有题的对数似然加起来，在 (a, β) 的网格上找最大值。 注意 c = 0 的题（n 次全错）也提供了信息：它们告诉我们 0 附近有多少质量。 【参数】counts：每道题的成功次数（list）；n：每题的样本数 【返回】(a, β) 【边界】网格是对数均匀的：a ∈ [0.02, 5]，β ∈ [0.1, 200]，各 60 个点。真实的最优值如果在网格外， 会被截在边界上（这时应该扩大网格）。 """ grid_a = [0.02 * (5 / 0.02) ** (i / 59) for i in range(60)] grid_b = [0.1 * (200 / 0.1) ** (i / 59) for i in range(60)] best = None for a in grid_a: for b in grid_b: ll = math.fsum(log_beta(c + a, n - c + b) for c in counts) - len(counts) * log_beta(a, b) if best is None or ll > best[0]: best = (ll, a, b) return best[1], best[2] def demo_forecast(rng): print("=" * 78) print("演示 2：toy world —— 只用每题 100 个样本，能不能预测 coverage@10000？") print("=" * 78) problems = make_problem_set(random.Random(30), n_problems=100, min_steps=3, max_steps=8) model = ToyLLM(1e9) n_big, n_small = 10000, 100 big_counts, small_counts = [], [] for p in problems: flags = [model.generate(p, rng).correct for _ in range(n_big)] big_counts.append(sum(flags)) small_counts.append(sum(flags[:n_small])) # 前 100 个样本：假装我们只采了这么多 # (a) 几道单独的题：失败率是否 ≈ (1 − p)^k？ print("(a) 三道单独的题（p 用 10000 个样本估计）：1 − pass@k 的实测值 vs (1 − p)^k") picks = sorted(range(len(problems)), key=lambda i: big_counts[i]) chosen = [picks[len(picks) // 2], picks[len(picks) // 8], picks[3]] print(pad("k", 8) + "".join(pad(f"题{problems[i].pid}(p={big_counts[i] / n_big:.4f})", 26, "right") for i in chosen)) for k in (1, 10, 100, 1000): row = pad(str(k), 8) for i in chosen: measured = fail_at_k_stable(n_big, big_counts[i], k) theory = (1 - big_counts[i] / n_big) ** k row += pad(f"{measured:.2e} / {theory:.2e}", 26, "right") print(row) # (b) 预测整体的 coverage fit_ks = [5, 10, 20, 50, 100] fit_cov = [mean(pass_at_k_stable(n_small, c, k) for c in small_counts) for k in fit_ks] direct, b_direct = fit_power_law_on_coverage(fit_ks, fit_cov) a_hat, beta_hat = fit_beta_binomial(small_counts, n_small) print(f"\n(b) 只用前 {n_small} 个样本：直接拟合的指数 b = {b_direct:.3f}；" f"beta-binomial 拟合 a = {a_hat:.3f}，β = {beta_hat:.2f}") print(f" 前 {n_small} 个样本里全错的题：{sum(1 for c in small_counts if c == 0)} 道；" f"10000 个样本里全错的题：{sum(1 for c in big_counts if c == 0)} 道") print(pad("k", 8) + pad("真实 coverage", 16, "right") + pad("直接拟合外推", 16, "right") + pad("分布拟合", 14, "right") + pad("真实 1−cov", 14, "right")) for k in (10, 100, 1000, 3000, 10000): truth = mean(pass_at_k_stable(n_big, c, k) for c in big_counts) pred_direct = direct(k) pred_dist = 1 - beta_failure(a_hat, beta_hat, k) print(pad(str(k), 8) + pad(f"{truth:.4f}", 16, "right") + pad(f"{pred_direct:.4f}", 16, "right") + pad(f"{pred_dist:.4f}", 14, "right") + pad(f"{1 - truth:.4f}", 14, "right")) print("\n观察：") print(" · (a) 每道题单独看，失败率和 (1 − p)^k 基本吻合：每道题都是指数衰减。") print(" k = 1000 那一行实测值偏小一些：k 已经是样本数的 1/10，\"不放回地抽 k 个\"的无偏估计") print(" 和\"有放回\"的 (1 − p)^k 不再相同（第 1 讲练习 1）—— 这是估计方法的性质，不是模型的性质。") print(" · (b) 在 k ≤ 100（样本覆盖的范围）内，两种方法都和真实值接近。") print(" 再往外推，两种方法都偏悲观，而且彼此差不多：k = 1000 时真实的 1 − coverage 是 0.027，") print(" 两种预测都在 0.07 左右；k = 10000 时所有题其实都被解决了，预测却还剩 2~3% 没解决。") print(" · 原因：前 100 个样本里有 21 道题一次都没做对。它们的 p 到底是 0.005 还是 0.000001，") print(" 100 个样本给不出这个信息。拟合出的 Beta(a ≈ 0.44) 假设 p 的分布一直延伸到 0，") print(" 但这组题里最难的一道 p 也有约 0.001 —— 左尾在这里\"截断\"了，就像演示 1 的 Uniform 例子。") print(" （真实的 benchmark 也可能反过来：有些题模型根本做不对，p = 0，coverage 永远到不了 1。）") print(" · 结论：任何外推都依赖对\"没见过的那部分难题\"的假设。分布的视角告诉我们该关心什么") print(" （左尾的形状），但没有方法能凭空得到从未解出过的题的信息。\n") def main(): demo_beta_math() demo_forecast(random.Random(31)) if __name__ == "__main__": main()