错过的蝴蝶与经验贝叶斯

我最早在 Efron 的教科书 Computer-Age Statistical Inference 看到这个故事. 作为经济学从业者, 我读过很多故事, "错过的蝴蝶"仍是我读过的最好的故事之一.

科贝特捉蝴蝶

二战期间, 年轻的博物学家科贝特仍在英属马来亚 (今天的马来西亚) 的小岛捉蝴蝶. 这个小岛有很多科贝特从未见过的珍稀蝴蝶物种. 科贝特在马来亚捉了两年蝴蝶, 一共捕捉到了 501 个珍稀蝴蝶物种.¹

我们用 k=1,2,...,501 来标记蝴蝶种类. 对物种 k, 科贝特捕获到 $x_k \ge 1$ 个标本. 因此, 科贝特这两年的收获可以表示为 $\{ x_k \}_{k=1}^{501}$.

对任意珍稀蝴蝶物种 k, 科贝特捕获的标本数都不超过 24 ($\forall k$: $1 \le x_k \le 24$). 符号 $y_x$ 表示恰好被捕获 $x$ 次的物种个数: $$y_x = \# \{k : x_k = x\}, \quad x = 1, 2, ..., 24$$

下表记录了科贝特这两年捕蝶生涯的收获. 可以验证, $\sum_{x=1} ^{24} y_x = 501$.

x  1   2  3  4  5  6  7  8  9  10 11 12
y  118 74 44 24 29 22 20 19 20 15 12 14

x  13 14 15 16 17 18 19 20 21 22 23 24
y  6  12 6  9  9  6  10 10 11 5  3  3

两年的捕蝶生涯过后, 科贝特不得不离开马来亚返回英国. 回到英国后, 科贝特仍对马来亚的蝴蝶恋恋不舍, 他经常问自己:

如果当初再多待一年, 可以发现多少新的珍稀蝴蝶物种呢?

时光不能倒流, 科贝特永远无法知道这个问题的真正答案. 为此, 科贝特请教了很多人, 包括当时在剑桥工作的统计学家费雪. 费雪为科贝特的这个问题给出了满意的答案.

科贝特捕蝶: 初步分析

记科贝特所在的岛上一共有 S 个珍稀蝴蝶物种. 参数 S (>501) 是未知的, 因为它包括科贝特观察到的和没观察到的蝴蝶物种.

我们用泊松过程来描述"捕蝶"这个动态过程. 由于科贝特在岛上待了两年, 用两年作为一个"时间单位", 物种 $k \in \{ 1,...,S \}$ 在一个时间单位内被捕获的次数 $x_k$ 服从参数为 $\theta_k$ 的泊松分布. 其中, 参数 $θ_k$ 表示物种 k 蝴蝶的到达速率 (或捕捉速率), 并且 $𝔼 [x_k] = θ_k$.

考虑一个额外的捕捉时间段 $t$. (科贝特之问中, $t = 1 / 2$ 时间单位 (即一年), 不过我们下面的分析仍然把 $t$ 看成变量, 允许其取值为任意正数) 令 $x_k(t)$ 表示新时间段内物种 $k$ 被捕获的次数. $x_k(t)$ 服从到达速率为 $\theta_k t$ 的泊松分布, 且 $x_k(t)$ 独立于此前的捕捉数 $x_k$.

小结: 我们使用泊松分布来描述物种 $k$ 的捕捉次数, 模型的参数为 $\{ S, θ_1 , ..., θ_S \}$.

• 第一个时间段物种 k 捕捉数 $x_k \sim \text{Poi} (θ_k)$,
• 第二个时间段物种 k 捕捉数 $x_k(t) \sim \text{Poi} (θ_k t)$.

新物种数: 事前与事中

接下来, 我们将科贝特之问中的参数表示为 $θ_1 , ..., θ_S$ 的函数.

• 记事件 $\text{Event}_k =$ {物种 k 在初始捕捉时间段内未被捕获, 但在新时间段内被捕获}.
• 由独立性, $p_k := \Pr [\text{Event}_k] = e^{-\theta_k}(1 - e^{-\theta_k t})$
• 令 $ξ_k$ 为 $\text{Event}_k$ 的示性函数: 若事件 $\text{Event}_k$ 发生, $ξ_k=1$; 否则, $ξ_k=0$.
• 第二个时间段内捕捉的新物种数为 $ξ = \sum_{k=1}^S ξ_k$. 其期望为: $$\sum_{k=1}^{S} p_k = \sum_{k=1}^{S} e^{-\theta_k}(1 - e^{-\theta_k t}) \equiv E(t)\qquad (1)$$

记上式的值为函数 $E(t)$, $t > 0$. 如果我们称第一个时间段开始前为事前 (ex-ante), 第一个时间段和第二个时间段之间为事中 (interim), 第二个时间段之后为事后 (ex-post), 这里的 $E(t)$ 是我们对科贝特之问的事前回答.

事中回答. 我们用编号 $k=1,...,501$ 表示科贝特第一个时间段捕到的蝴蝶种类, $k=502,...,S$ 表示科贝特没有捕到的蝴蝶种类. 科贝特之问的事中回答是: $$\sum_{k=502}^S 𝔼 [(1 - e^{-\theta_k t}) \mid k \text{ unobserved.}] \equiv \bar{E} (t)$$

无论是事前回答 $E(t)$ 还是事中回答 $\bar E(t)$, 根据已有的数据都不好估计. 这是因为, 科贝特之问关心的是新的珍稀蝴蝶物种, 这涉及到上表中缺失的条目 $x = 0$:

x 0     1   2  ...
y S-501 118 74 ...

我们称这个问题为缺失物种 (missing species) 问题: 表中没有提供 x=0 对应的数据 ($y_0$).

后文中, 我们默认科贝特之问讨论的是事前期望 E(t), 而不是事中期望 $\bar E(t)$. 我不知道为什么费雪选择的是事前期望 E(t), 我猜是因为它比 $\bar E(t)$ 好估一些.

估计 E(t): 非参数方法

我们首先介绍一种估计 E(t) 的非参数方法: 可以绕过 $\{ θ_1,...,θ_S \}$, 直接用已有数据 $\{ y_1,...,y_{24} \}$ 来估计 E(t).

1. 首先, 将 E(t) 的表达式从求和改成积分形式: $$E(t) = S \sum_{k=1}^{S} [e^{-\theta_k}(1 - e^{-\theta_k t}) \frac 1S]$$ $$= S \int_{0}^{\infty} e^{-\theta}(1 - e^{-\theta t}) d G(\theta) \qquad (2)$$

   • 其中 $G(\theta)$ 表示参数 $θ$ 的经验分布, 这个分布给每个情况 $[θ=\theta_k]$ 赋予概率 $1 / S$. 由于 G 是离散的, 我们这里的积分是勒贝格-斯蒂尔杰斯积分.
   • 离散分布 G 理论上不存在密度函数. 但出于符号的方便, 我们仍用 $g(\cdot)$ 表示其对应的概率密度, 并称 $g(θ)$ 为参数 $θ$ 的经验密度. 读者将后面的 "$g(\theta) d \theta$" 替换为 "$d G(\theta)$" 即可.

2. 对 $1 - e^{-\theta t}$ 进行级数展开: $$E(t) = S \int_{0}^{\infty} e^{-\theta}\left[\theta t - (\theta t)^2 / 2! + (\theta t)^3 / 3! - \cdots\right] g(\theta) d \theta \qquad \text{(△)}$$

3. $y_x$ 的期望为: $$\begin{aligned} e_x &= E [y_x ] = \sum_{k=1}^{S} e^{-\theta_k} \theta_k^x / x! \\ &= S \int_{0}^{\infty}\left[e^{-\theta} \theta^x / x!\right] g(\theta) d \theta \qquad \text{(□)} \end{aligned}$$

   • 注: 这里我们再次用到了经验密度 g 这个代数小技巧.

4. 比较方程 ($\text{△}$) 和 ($\text{□}$), 发现 $E(t)$ 可以表示为 $$E(t) = e_1 t - e_2 t^2 + e_3 t^3 - \cdots$$

上面这个式子意义重大:

• 为了回答科贝特之问, 我们真正感兴趣的参数是 E(t), 但 E(t) 的值依赖经验分布 $g(θ)$, 直接估计 $g(θ)$ 难度较大.
• 现在, 新的表达式中 $E(t)$ 只依赖 $e_1$, $e_2$, ... 并且, 我们可以使用替换原则, 用科贝特表中的 $y_x$ 替换 $e_x$: $$\hat{E}(t) = y_1 t - y_2 t^2 + y_3 t^3 - \cdots$$

科贝特之问中, $t = 1 / 2$. 因此, $\hat{E}(1 / 2) = 118 \times (1 / 2) - 74 \times (1 / 2)^2 + 44 \times (1 / 2)^3 - \dots = 45.2$.

作为学术统计学家, 费雪自然不会只满足于给出参数的点估计, 至少也要给出替换估计量 $\hat{E}$ 的标准误吧. 如果 $x_k$ 相互独立, 可以用下面的值作为 $\hat{E}(t)$ 的近似标准误: $$\widehat{\mathrm{sd}}(t) = \left(\sum_{x=1}^{24} y_x t^{2x}\right)^{1 / 2}$$

对于 t=0.5, 计算得到的标准误为 11.2. 如果我们想得到一个区间估计, 可以用 $\hat{E}(0.5) = 45.2 \pm 11.2$.

• 也就是说, 如果当初科贝特再多待一年, 他会多捕获 34--56 个新的蝴蝶物种!

小结:

• 我们首先用了一个代数小技巧, 将问题的参数由 $S, θ_1, θ_2, ..., θ_S$ 转化为 $θ$ 的经验分布 $g(θ)$.
• 接下来, 我们绕过 $g(θ)$, 将感兴趣的参数 $E(t)$ 直接表示为 $e_x$ 的函数. 这样做的好处是, 我们可以直接用现有数据 $y_x$ 替换掉 $e_x$, 得到替换估计量 $\hat E(t)$.

费雪的方法: 参数化 g(θ)

在费雪对"缺失物种问题"的回答中, 他实际上使用了一个参数化模型, 其中密度函数 $g(\theta)$ 为伽马函数形式: $$g (θ) = \frac {θ^{\nu - 1} e ^ {-θ/σ}}{σ^\nu Γ(\nu)}$$

费雪进一步证明了, E(t) 可以表示为 $e_1$ 和伽马函数参数的表达式: $$E(t) = e_1 (1 - (1 + γt )^{-\nu}) / (γ \nu)$$ 其中 $γ = σ/(1+σ)$. 费雪用 $y_1$ 来替换 $e_1$, 用对应的极大似然估计来替换 $\nu$ 和 $σ$, 进而得到了 E(t) 的估计.

对比前面提到的非参方法和费雪的参数化方法:

• 非参方法没有规定密度函数 $g(\theta)$ 的形式, 因此其假设更少, 模型自由度更大. 并且, 非参方法和参数方法在 $t \in [0,1]$ 时的预测结果很接近.
• 对 $t > 1$ 的情形, 参数化方法的预测结果相对更"合理". 这是因为非参方法的预测结果 $\hat{E}(t)$ = $y_1 t - y_2 t^2 + y_3 t^3 - \dots$ 在 $t>1$ 时取值很不稳定.
• 越灵活的模型, 估计 (or 预测) 结果的标准误一般也越大.

经验贝叶斯

请读者思考, 本文上面介绍的非参方法, 属于贝叶斯方法吗?

• 注: 贝叶斯学派和频率学派的区别在于: 前者假设参数是服从某个先验分布的随机变量; 后者眼中参数就是参数, 不是随机的.
• 有博弈论基础的读者, 可以看这篇博文

抛开观念上的争论, 只看数学本身的话, 这里的非参方法不是贝叶斯方法. 我们的 g(θ) 是参数 θ 的经验分布, 而不是先验分布: E(t) 原本的求和表达式 (1) 中根本没用到积分和 g(θ), 将求和表达式 (1) 转化为积分表达式 (2) 这一步是纯代数的, 不必引入 "参数 θ 服从某个先验分布" 这么庞大的假设.

但是, 这里把 g(θ) 解读为参数 θ 的先验密度似乎也没问题. 并且, 这个解释在数学上更自然, 能更好地看到 E(t) 积分表达式的动机.

费雪的参数化方法中, 假设 g(θ) 函数为伽马函数形式, 这就更有"贝叶斯味"了. 这类方法一般被称作经验贝叶斯 (empirical bayes) 方法.

(传统)贝叶斯方法一般包括如下步骤:

1. 选择某个统计模型来描述数据 x 的生成机制, 这个模型通常可描述为概率分布 p(x|θ).
2. 假设参数 θ 服从某个先验分布 g (或使用共轭先验)
3. 使用贝叶斯公式计算参数 θ 的后验, 即: 后验 p(θ|x) $∝$ 先验 g(θ) $×$ 似然 p(x|θ). 利用后验分布 p(θ|x) 进行预测或推断.
   • 贝叶斯方法通常计算量较大, 且对于较复杂的模型难以得到解析解. 比较流行的做法是用 MCMC 进行数值求解.

不同于贝叶斯方法, 经验贝叶斯中先验分布的具体函数形式和数据有关. 当然, 如何利用数据来估计先验分布就见仁见智了. 费雪假设先验分布为伽马分布, 并用极大似然法来估计这个伽马分布的参数.

• 我个人的理解是, 贝叶斯方法是一套完整的工作流: 从先验和似然性出发, 到计算后验分布. 而 "经验贝叶斯" 更像是一种"后门" (backdoor) 技巧, 它允许先验和(经验)数据有关.

最后, "经验贝叶斯"这个叫法和经济学中的"垄断竞争"很像, 二者都属于矛盾修辞 (oxymoron). 贝叶斯方法的核心在于先验分布, 而"先验" = "先于经验". 如果康德看到 empirical prior 这个表达, 肯定会斥之邪道: 先验的事物是不能受到经验数据"污染"的. (当然, 今天的统计学似乎是"实用主义"至上, 好用的模型才是好模型)

进一步阅读

本文主要内容都搬运自文章开头提到的 efron 教科书 Ch6. 除了这个捉蝴蝶的故事外, efron 还提供了 (1) 保险公司预测来年理赔数量和 (2) 鉴定莎士比亚真迹的例子. 其中, 保险公司的预测问题中所有投保人数 (即捉蝴蝶问题中的所有珍稀蝴蝶物种数) S 是已知的. 感兴趣的读者应该读一读这个例子, 并思考保险公司预测问题中给出的估计是"事前估计"还是"事中估计".

杨灿教授有一篇关于经验贝叶斯的中文介绍, 用的例子是 James-Stein 估计量.

如果你有幸能访问互联网, 维基关于经验贝叶斯的介绍很好. (我个人的观察是, 维基上统计和机器学习的词条质量普遍都不错, 比{数学,经济学,博弈论}的质量高很多. 当然, 作为非专业统计学研究人员, 我的判断肯定是有偏的)