错过的蝴蝶与经验贝叶斯

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

科贝特捉蝴蝶

二战期间, 年轻的博物学家科贝特仍在英属马来亚 (今天的马来西亚) 的小岛捉蝴蝶. 这个小岛有很多科贝特从未见过的珍稀蝴蝶物种. (注: 我查阅了维基百科, 蝴蝶在生物分类学中是指鳞翅目中名为凤蝶总科的一个总科级演化支, 又名真蝶总科. 因此, 蝴蝶系"科"级分类, 其下有不同的(物)种.)

科贝特在马来亚捉了两年蝴蝶, 一共捕捉到了 501 个珍稀蝴蝶物种. 我们用 k=1,2,...,501 来标记蝴蝶种类. 对物种 k, 科贝特捕获到 \(x_k \ge 1\) 个标本. 因此, 科贝特这两年的收获可以表示为 \(\{ x_k \}_{k=1}^{501}\).

对任意珍稀蝴蝶物种 k, 科贝特捕获的标本数都不超过 24 (即 \(1 \le x_k \le 24\), \(\forall k\)). 用 \(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}) =: 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.}] =: E\text{-interim}(t) \]

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

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

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

后文中, 我们默认科贝特之问讨论的是事前期望 E(t), 而不是事中期望 E-interim(t). 我不知道为什么费雪选择的是事前期望 E(t), 我猜是因为它比 E-interim(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 (*) \]

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 (**) \end{aligned} \]

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

4. 比较方程 \((*)\) 和 \((**)\), 发现 \(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 估计量.

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