漫谈统计学习之经验贝叶斯（Empirical Bayes）

2022 年 3 月 23 日 PaperWeekly

©作者 | 黄春喜

学校 | 香港科技大学

研究方向 | 智能交通

1.1 前言

在传统的贝叶斯方法中，我们在获取到任何观测数据之前就已经有了确定的先验分布，然后利用观测数据对先验分布进行修正得到后验概率。举个很简单的例子：我在乘坐飞机之前就知道飞机失事的概率非常小，这个概率在我乘坐飞机之前就已经在我心里了；但如果观察到最近国际局势紧张（观测数据），此时我所乘坐的飞机失事的概率（经过修正后的后验概率）可能就会增加。

与此不同的是，经验贝叶斯方法中先验分布并不提前给出，而是从获得的观测数据中去估计，这恰好印证了其命名中的“经验”二字。

在统计推断中，经验贝叶斯方法是非常重要的工具之一，能够帮助我们在缺乏先验分布的情况下做出推断，有时有些场景下的推断甚至性能很不错。在过去的研究中，已经有众多令人称奇的故事和结果出现在我们的视野里，在这里我打算分享一些关于经验贝叶斯的理解。

本文第一部分主要分享经验贝叶斯出现、发展历程中的一个重要节点：Robbins' formula，它证明简单，但意义重大。后面两部分将分享关于 the missing species problems 以及 James-Stein estimator 的相关内容。

1.2 Robbins’formula

Robbins, H. (1956). An empirical Bayes approach to statistics. In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1954–1955, Vol. I, Berkeley and Los Angeles, pp. 157–163. University of California Press.

考虑以下推断场景：

对于一个未知的概率密度（此处的密度也包含离散变量的概率），其产生了一个随机样本，其中包含，且这些都是无法观察到的。对于每一个，它们都根据已知的概率密度产生了一个可以观察到的随机变量：

且之间是相互独立的。从可观察到的，我们希望估计参数。

对于这个问题，不难发现，如果利用极大似然估计我们很容易得到，很明显这样的话每一个观测到的随机变量只对估计有贡献。

但在上述这篇文章中，Robbins提出了一个具有里程碑意义的观点：除了以外的其他观测值能够为我们估计感兴趣的参数提供帮助！

乍一听觉得有些不可思议，以下是一个可以用来阐明 Robbins 的观点的例子：

一个汽车保险公司总结了去年它的 9461 位顾客中每位顾客的理赔数量，并列出了理赔数和不同理赔数所对应的顾客数。如下表第一行和第二行所示，很明显，在 9461 位司机中，有 7840 位去年驾驶非常安全、一次理赔也没有提过，而与之相反的是有 1 位司机，竟然在短短一年的时间中理赔了 7 次。

为了制定下一年的汽车保险价格计划，该汽车保险公司现在希望能够预测一下这些顾客下一年大概要理赔的次数，这样的话可以针对性的提前调整保险价格来保证保险公司的收益率。

所以问题来了，对于去年一次都没有理赔过的 7840 位司机而言，在新的一年里他们平均会理赔多少次呢？对于去年理赔过一次的 1317 位司机来说，新的一年他们平均又会理赔多少次？……在这里，我们的目标是对进行一个预测。

在此，我们先给出 Robbins' formula，看看他的预测方式，再详细解释和说明：

举个简单的例子：对于那 7840 位去年没有理赔过的司机，根据 Robbins'formula，他们下一年平均会理赔的次数是，以此类推可得到上表中第三行中的数据。

刚看到这里时，可能会觉得有些诧异，我们感兴趣的是那 7840 位去年没有理赔过的司机，怎么跟理赔次数为 1 的 1317 位司机又有关系了呢？Robbins' formula 中的这种关系究竟是如何建立起来的呢？

1.3 Proof of Robbins' formula

以下简单推导 Robbins' formula：

首先假设之前提到的是泊松的，可以写成：

。

对于，根据贝叶斯定理可得：，其中为的边缘密度，是所有可能的值的集合。

代入泊松密度我们可得：

于是，我们可以通过积分求得：

在经验贝叶斯中本身和是未知的，但它们很好估计：在这里的期望值就是。所以，当我们将以及代入上式时，自然就得到了我们上述的 Robbins' formula：

证毕。

1.4 关于Robbins' formula的几点说明

1. 在以上这个汽车保险公司的例子中：在没有任何先验分布/信息的情况下，简单的一行数据（）中就提供了所有用于估计一个贝叶斯统计值（）的信息。某种意义上来说，这种不依赖先验分布/信息的估计是完完全全的频率学派的估计方法，这也恰恰印证了这种方法的大名：经验贝叶斯。不管怎样，任它频率学派还是贝叶斯学派，事实上 Robbins' formula 就是在没有贝叶斯信息（先验分布）的情况下得到了一个贝叶斯的结果。这也给了后续的研究一个革命性的启发：一个足够大的数据集中隐含了其本身的先验分布。

2. Robbins'formula：，很明显，在上述例子中对于一位司机而言，这个估计不仅受到他/她自己过去的经验的影响，还受到别的司机过去经验的影响。从公式中可以看到，一个很安全的司机（理赔数为 0）的估计只受比他/她不那么安全一点（理赔数为 1）的司机的经验影响。这里其实有一个隐含的假设：别的信息要有一定的相关程度才会被用于帮助预测，以后有空可以再写篇文章探索一下这里。

3. 从本质上来说，Robbins's formula 是一个非参数式的方法：它并不依赖任何参数式形式的先验分布。事实上在上文推导过程中，我们也从未去估计这个先验分布。

4. 有一点需要注意的是，在上文汽车保险公司例子的表中，第三行所示为利用 Robbins' formula 得出的估计结果，我们发现当的值变小的时候（表的右侧），这种估计好像出问题了。个人认为这与数据分布本身有关，毕竟从概念上来说，经验贝叶斯的先验分布来源于观测数据本身。而与之相反的则是表中第四行假设先验分布为 gamma 分布得出的估计值，相较之下更为稳定。

2.1 前言

上一章节我介绍了关于 Robbins‘ formula 的相关内容。很明显，所谓的经验贝叶斯实际上还是在贝叶斯框架下做事，唯一的区别就在于先验分布的给出上：传统的贝叶斯方法直接给出了明确的先验，而经验贝叶斯则依靠观测数据来估计先验分布、进而实现贝叶斯推理。

虽然 Robbins 正式命名了经验贝叶斯方法、并将其概念阐述得很清楚，但实际上，在 Robbins 之前的一篇更早的论文可以视为经验贝叶斯的里程碑。

Fisher, R., A. Corbet, and C. Williams (1943). The relation between the number of species and the number of individuals in a random sample of an animal population. J. Anim. Ecol. 12, 42–58.

这就是本文即将讲述的有关经验贝叶斯的另一个成功故事：the missing species problem。此外，我将以过去的研究为例，简单介绍经验贝叶斯框架下的两种方法：f-modeling 和 g-modeling。

在接下来的第三篇文章，我将继续讲述关于经验贝叶斯的另一个故事：James-Stein estimator。

2.2 The missing species problem

回到问题本身，the missing species problem 主要考虑的场景如下：

Corbet 是一个自然学家，他曾经花了两年时间在马来西亚捕获蝴蝶物种标本。以下表格是他所收集到的部分实验数据：其中代表的是捕获到某种类蝴蝶的次数，所代表的是捕获次数为的蝴蝶的种类数。

从上表中不难看出：在两年的时间里，有 118 种蝴蝶非常非常稀有，它们只被捕获过 1 次；有 74 种蝴蝶总共在两年内被捕获了两次；与之相反的是，有 3 种蝴蝶非常常见，在两年时间里被捕获了 24 次。与我在第一篇文章中用来推导 Robbins‘ formula 的汽车保险公司的例子相比，这里最大的区别就是：缺少了 x=0 时的数据。x=0 意味着什么呢？意味着那些我们没有捕获到的、没有见过的蝴蝶品种，也就是问题中所说的“missing species”。

花了两年时间弄到这些数据，Corbet 心里充满了好奇，自然界如此之大，到底还有多少蝴蝶是没见过的呢？于是他问了 20 世纪最著名的统计学家之一 Fisher 一个问题：如果我继续再在马来西亚多待一年的时间，还能捕获到多少之前没见过的新品种的蝴蝶呢？这个问题乍一看好像很无厘头，似乎没法回答。

而 Fisher 给出了他的答案。从我们今天的角度来看，他的答案无疑就是参数化的经验贝叶斯方法。以下：

2.3 Fisher‘s answer

Fisher 的答案是这样给出的：

首先，我们用来对蝴蝶物种进行编号的话，其中是包括观察到的和没有观察到的蝴蝶物种总数。假设蝴蝶物种被捕获的过程服从一个参数为的泊松过程。

令代表一个单位时间（在本例中为两年）内蝴蝶物种所被捕获的次数，令代表未来 t 个单位时间内蝴蝶物种所被捕获的次数。根据以上的泊松过程假设，我们有：

其中后面一个式子代表的是物种在原先的捕获阶段没有发现、但在未来的捕获阶段中被发现了的概率，此时的物种所代表的不就是 Corbet 所想要的新的没有见过的“missing species”吗？

先验分布可以被看成一个经验分布：个参数中的每一个都有相同的概率。这样的话我们可以把在额外的 t 个单位时间内的能捕获到的新物种的期望值写为：

此处忽略本身在这个例子中的离散性。

为了解决这个问题，Fisher 假设泊松过程参数的先验分布是一个 scaled gamma form。可以写成：

其中。

实际上就是，是一个形状参数为的标准 gamma 变量。在这种情况下，把参数从利用来替换可以将参数变为，那么边际概率密度可以写成：

将上表中的数据作为一个多项式样本来考虑，我们可以利用极大似然估计来最大化 log likelihood：（本例中从 1 取到 24）

利用上式即可从观测到的数据中估计出，的值，利用的估计值可以算出的估计值，进而得出一整个先验分布。在本例中，，。

再将这个由观测数据所估计得来的先验分布代入

即可计算出。在本例中，Fisher 利用 gamma 模型作为先验得出的为：

通过公式，Fisher 就成功的回答了 Corbet 的问题，如果再多待一年的话，t=0.5（一个单位时间为两年），则有。换句话说，根据 Fisher 的方法，如果 Corbet 再多待一年的时间，他将有可能再多捕获 40+ 种新的蝴蝶物种。

Fisher 的方法中，他预先假设了先验分布，然后利用观测数据去估计了先验分布中的参数，得到估计的先验分布后利用贝叶斯定理解决了问题。这无疑是典型的经验贝叶斯方法。

然而，针对这一问题后来又有研究者给出了另一种经验贝叶斯的解决方法。

Good, I. and G. Toulmin (1956). The number of new species, and the increase in population coverage, when a sample is increased. Biometrika 43, 45–63.

2.4 Good and Toulmin’s answer

与 Fisher 的解决方法完全不同的是，Good 和 Toulmin的解决方法是完全非参数化的方法。换句话说，他们的方法并不对先验分布做预先的假设。据 Good 本人的说法，这项工作也得到了 Alan Turing 的帮助和启发。

首先，他们同样假设蝴蝶物种被捕获的过程服从一个参数为的泊松过程。所以有以下关系成立：

同样地，他们将先验分布看成一个经验分布：个参数中的每一个都有相同的概率。这样的话在额外的t个单位时间内的能捕获到的新物种的期望值可写为：

此处对进行泰勒展开，可以得到：

而与此同时，对于总共个 species，将对于每一个 species 可以累加得到的期望值，可以写成：

比较以上两式可以将写为：

在此式中，可以由来替代（是的一个无偏估计），所以我们有：

利用上式，我们同样可以回答 Corbet 的问题，再多待一年的话，

如果再多待两年的话，依然可以利用上式轻松算出。

很明显可以看出，在 Good 和 Toulmin 的方法中，他们没有去对先验分布进行假设、利用数据进行估计，而是直接通过非参数式的方法解决了这个问题。

2.5 关于两种方法的说明

下图为两种方法的 solution path 的比较：

红色虚线代表的是 Fisher 的方法，黑色实线代表的则是 Good 等人的方法。两种不同方法给出的答案竟然在的范围内有着惊人的相似！

有一点需要说明的是，黑色实线所代表的非参数式的经验贝叶斯方法在时会存在类似收敛的行为，这是有缺陷的，但红色虚线所代表的参数式的经验贝叶斯方法并不会存在此类问题。事实上，在上一篇文章中针对 Robbins' formula 的例子中，我们也发现，利用非参数式的方法到了表格右侧（数据集末尾）时同样出现了问题，无法给出合理的结果。

对于这种现象，我本人的理解是这样的：首先对于参数式的经验贝叶斯方法，先验分布被假设后，其中的参数被观测数据很好的估测出来了，此时的先验分布有着非常明确的模型表达，在用于计算 posterior、乃至其他值时有非常好的性质，可以给出明确的模型表达，所以当我们对变量进行变化时，由参数式经验贝叶斯得出的结果都是看起来比较 reasonable 的。

与之相反的是，非参数式的经验贝叶斯或利用观测数据绕过先验分布直接去找边际分布、或利用一些 trick 来得出后验的一些表达，但这些都是限制在我们已观测到的数据上的，一旦超出我们已观测数据的范围，这种非参数式的方法就很容易出问题。换句话说，当我们没有足够的数据支持时，所谓的经验贝叶斯自然就不那么 reliable 了。

上图中时就是一个很好的例子，当时，我们已有观测到的一个单位时间内（两年）的数据，所以非参数式的方法有足够的数据支持，表现得很好，但当变量的范围超出了我们手中的观测数据时，经验贝叶斯也就无处可以提取经验了。这一问题，后续我还会继续深入了解，也欢迎各位伙伴一起讨论。

事实上，Fisher 和 Good 的方法实际上对应了经验贝叶斯方法中的两个类别，Fisher 这种利用数据估计先验概率的参数式方法被称为 g-modeling，而 Good 这种更关注边际分布的非参数式的方法被称为 f-modeling。个人感觉 f-modeling 很需要灵感和问题的特别性，而 g-modeling 更像一种易用的方法，有着清晰的框架和解决问题的思路。

3.1 前言

前面两个章节中，我分别介绍了经验贝叶斯方法的两个重要历史节点：Robbins' formula 以及 The missing species problem。

以上两个极具代表性的早期经验贝叶斯方法主要是基于假设观测值是泊松分布的，而本文即将涉及的 James-Stein estimator 则是基于观测值是正态分布的假设。

James, W. and C. Stein (1961). Estimation with quadratic loss. In Proc. 4th Berkeley Sympos. Math. Statist. Prob., Vol. I, pp. 361–379. Berkeley: University of California Press.

作为曾经一度颠覆极大似然估计、刷新统计学界认知的方法，James-Stein estimator 距今已有几十年的历史，但其思想和拓展依旧在当今的统计学习方法中熠熠生辉。基于此，本文主要介绍 James-Stein estimator 的基本概念、推导过程以及 James-Stein Theorem 的证明。

3.2 What's it? - The baseball player example

James-Stein estimator 到底做了一件什么事？不同在哪里？可以考虑以下统计推断场景：，其中，是观测值，是未知参数。

在这种情况下，如何根据去估计呢？

传统最直接的方法就是利用 maximum likelihood estimator，此时易得：。

然而 James-Stein estimator 给出的估计却是这样的：

。

在上述推断场景中，简单比较以上两种估计方法不难发现以下两点区别：

极大似然估计给出的估计结果中每一个只用到了其本身对应的一个观测值，而 James-Stein estimator 给出的估计结果中每一个用到了所有的观测值；
相较于极大似然估计来说，James-Stein estimator 给出的估计结果小于本身，换句话说，括号中减去的那一项起到了一个缩小估计量的作用。

以下利用一个实例来简单展示 James-Stein estimator 的魅力所在：（推导和证明见下文）

下表是 1970 年的赛季中 18 位棒球击球手的平均击球数据：

其中第一列是 18 位击球手在该赛季前 45 次击球中的平均击球数（45 次击球大约占击球手一个赛季击球次数的 1/10 左右），而平均击球数指的是平均每击一次球成功的次数。以第一位击球手 Clemente 为例：在总共 45 次击球中，他成功了 18 次，所以他的平均击球数就是 18/45=0.4。在这种情况下，所观察到的 45 次击球中的平均击球数可以作为该击球手在该赛季剩余击球次数中平均击球数的极大似然估计，也就是。

而下表第二列则是利用 James-Stein estimator 所得到的估计结果：。

最后，表格第三列中，是该赛季结束后所统计到的最终每一位击球手在该赛季中剩余的击球次数中各自的平均击球数。

而表中最后一列则是每位选手在该赛季中剩余的击球次数。

基于以上数据，我们可以简单计算出 maximum likelihood estimator 和 James-Stein estimator 各自的 prediction error，如下：

这差距！很明显，在这个例子中，James-Stein estimator 把 prediction error 降低了两倍多！或许对于这个例子中预测棒球选手的击球数来说，这显得无足轻重，可如果这是在预测某种疾病治愈的概率呢？这时候这种差距就是生与死的区别了。毫无疑问，James-Stein estimator 的预测效果是我们需要重点研究的对象。

James-Stein theorem 对于预测效果给出了答案：

当时，

没有其他任何假设，没有其他任何条件，在 square loss 的情况下，James-Stein estimator 是更好的选择。

那 James-Stein estimator 又是如何做到这一点的呢？为什么 James-Stein estimator 被视作经验贝叶斯方法中的一种呢？以下 James-Stein estimator 的推导以及 James-Stein theorem 的证明将回答以上问题。

3.3 Derivation of James-Stein estimator

以下对 James-Stein estimator 进行简要推导：

首先，我们有

是已知的。，其中。

假设先验分布为，则的边际密度分布可写为：。

根据贝叶斯定理和高斯分布的相关性质（marginal, conditional, and partitioned Gaussians），可以简单的总结为：万物皆高斯，换句话说，如果边际分布是高斯的，其条件分布、边际分布、条件分布也全都是高斯的。由此可以得出的后验分布为：

故

至此，我们写出了 James-Stein estimator 的估计式，然而其中有一个未知的参数，而这个参数正是先验分布中的参数，可以从我们的观测数据中估计而来，这也就提供了 James-Stein estimator 作为经验贝叶斯方法的重要证据：先假设先验分布，并利用观测数据估计先验分布中的参数，过程中使用贝叶斯定理进行统计推断。

上文中，我们有，对其标准化我们有：

其中是自由度为的卡方分布。

，那么就服从自由度为的逆卡方分布，根据逆卡方分布的性质，我们有

最终，我们可以用作为的估计，代入上式中，即可得到：

推导完成。

此外，在估计先验分布参数时，也可以直接使用极大似然估计，此时得出的结果是：

可以发现两者唯一的区别在于和。在实际应用中，当大到一定程度（>10）时，两者差距并不大，都可以被视作是合理的估计，所获得的最终估计结果相差不大。

3.4 Proof of James-Stein theorem

以下简述 James-Stein theorem 的证明思路：

首先，我们有：

对两边取期望有：（注意此时是的函数）

其中

利用 Stein identity：

（可由分部积分推得），可得：

把和分别代入可得：

由以上两式，可以看出，当时：

证毕。

另一方面，我们也可以通过 simulation 来直观的展示出这个 pattern，codes 很简单就不放在这了，本文的标题图即为 JSE 和 MLE 在 N 值不同情况下的对比。

3.5 总结

本文介绍了 James-Stein estimator 的来源，并分享了简单的推导过程和证明。从上面的推导中我们发现，James-Stein estimator 实际上是一种 shrinkage estimator，如下图所示，红色点是上文棒球例子中的观测到的45次击球次数中的平均击球数（此处可代表）,而绿色点则代表，蓝色点是该赛季剩余击球次数中的平均击球数（可代表 ground truth 的值）。

提到 shrinkage，不免会想到 Lasso 和 Ridge，实际上，James-Stein estimator 在某种程度上在实行统计推断时具有某种隐含的正则化的作用。

在以上三个章节中，我分别分享了 Robbins's formula、the missing species problem 以及 James-Stein estimator，其中经验贝叶斯的理念一直贯穿其中：假设一个先验，并用观测数据去估计这个先验（中的参数）。然而，这三个故事中，使用的先验分布涉及到了泊松分布、正态分布等不同的分布，不由得会使人好奇：如果假设别的先验，经验贝叶斯还会 work 吗？答案当然是肯定的，事实上我在以上三个章节中分享的只是皮毛，Tweedie's formula 会给我们一个更 general 的答案。此外，如果在相同推断场景对 James-Stein estimator 和 Tweedie's formula 进行对比的话，它们的 shrinkage 行为与 Lasso 和 Ridge 有着惊人的相似。如果有伙伴感兴趣的话下次抽空再跟大家一起探讨。

4.1 前言

上面分别介绍了几个经验贝叶斯发展历程中的重要节点，包括 Robbins' formula、The missing species problem、James-Stein estimator 等。在经验贝叶斯方法发展的初期，Robbins、Fisher 等大佬在先验上用泊松分布用的比较多，后面 James-Stein estimator 又往经验贝叶斯的先验池子里加入了正态分布，然而这些还远远不够。

对于贝叶斯相关的方法，有一个经久不衰的疑问：为什么用这个先验？用其他先验行不行？

针对这个问题，在经验贝叶斯的框架中 Tweedie's formula 就在一定的推断场景下给出了更普遍的结果：不对先验的形式做限制，只要有就行。本文将主要介绍 Tweedie's formula 的推断场景、其利用 exponential family 及 cumulant generating function 的推导过程以及其与我们熟知的 Ridge regression、Lasso 之间的联系。

4.2 Tweedie's formula

考虑 James-Stein estimator 中的推断场景：

，，我们希望从观测值中估计。

在推导 James-Stein estimato 的过程中，我们假设了先验是正态分布：，在利用观测数据估计了其中的参数后最终写出了的后验分布并用其期望作为估计值。

然而在 Tweedie's formula 中，不对先验的形式做特殊要求，只假设其存在：，目标同样是求后验分布的期望。Tweedie's formula 给出的后验分布的期望为：

其中是的边缘分布。

很明显，Tweedie's formula 给出的估计中第一项就是 MLE 的估计，而第二项可以被视为经验贝叶斯在 MLE 估计的基础上做的贝叶斯修正。

我第一次看见这个结论的时候，除了 amazing 没有别的词可以形容了。这个公式意味着：在当前考虑的推断场景下，不对先验做形式上的要求，只要其存在，我们就可以求出后验分布的期望，是不是很神奇？

Tweedie's formula 到底是怎么得到的呢？

4.3 Tweedie's formula的推导

在推导 Tweedie's formula 之前，需要强调的是 Tweedie's formula 给出的答案是后验分布的期望，这与获得整个后验分布是完全不同的两码事，两者的区别有点像我们常说的判别模型和生成模型之间的区别。

以下进入 Tweedie's formula 的推导：

首先我们有，。

由此可以写出的边际分布为：

，其中为标准正态分布的概率密度函数。

根据贝叶斯公式我们可以得到的后验分布可写为：

在此处，先验分布形式不清楚的情况下，我们无法给出后验分布的解析形式。但此处我们的目标是后验分布的期望，可以通过 exponential family 及其中的 cumulant generating function 巧妙地推导而来。

对于 exponential family 中的分布而言，其概率密度函数可以写成以下形式：

其中被称为 natural parameter，被称为 cumulant generating function。

我们常见的正态分布、泊松分布就是属于 exponential family 中的分布，例如：

对于正态分布而言：

很明显，此时，。

对于泊松分布而言：

显然，此时，。

对于 exponential family 分布中的概率密度函数，其满足下式：

对此式两边取导数可写为：

将

代入上式可简化得：

上式等号左侧的第二项可以简化写为：

由此，我们得到了关于 cumulant generating function 的一条重要性质：。

回顾我们的目标：后验分布的期望。从上式显然期望已经出现在我们的视野中，剩下的就是将后验分布写成 exponential family 的形式，找到 cumulant generating function 求导即可。

以下为根据后验分布找出 cumulant generating function 的过程：

在上式中，可以被视为 natural parameter，cumulant generating function 则为：

对其求导，利用我们有：

Tweedie's formula 推导完毕。

4.4 JSE、Tweedie's formula与Ridge、Lasso的联系

Ridge 和 Lasso 是我们很熟悉的正则化方法，分别利用 L2 和 L1 范数实现了对参数的 shrinkage 以实现限制模型复杂度避免过拟合的效果。其中 Ridge 对所有参数做普遍的 shrinkage，得到的是 dense 的、数值上都较小的参数；而 Lasso 则对参数做 soft shrinkage，一部分参数直接 shrinkage 到 0，得到了 sprase 的结果。

而经验贝叶斯中的两个重要节点：之前所提到 James-Stein estimator 的和本文分享的 Tweedie's formula 则惊奇地有与 Ridge 和 Tweedie's formula 相似的效果。如下图所示，黑色为 JSE 和 Tweedie's formua 的估计值，红色线为。

这联系多美！

4.5 总结

相较之前的几种方法对先验形式的要求，Tweedie's formula 算是在普适性上一定程度的进步。但这种普适性仍然是有条件的，在推导过程中，不难发现，exponential family 的 cumulant generating function 扮演了很重要的角色。换句话说，想要使用 Tweedie's formula，需要先看看推断目标的分布（本章节中为的后验分布）能否被写成 exponential family 的形式。尽管有这样的限制，这丝毫不影响 Tweedie's formula 成为经验贝叶斯历史上浓墨重彩的一笔。

本质上来说，JSE 和 Tweedie's formula 都属于 shrinkage estimator，能与 Lasso 和 Ridge 联系起来自然也就不奇怪了。很多方法在运用的过程中，明面上模型中没有正则化的项，但在求解过程中的一些步骤会起到隐含的正则化的效果。有时间可以再写写最近学习的跟正则化相关的内容~

5.1 前言

在学习的过程中不禁感叹经验贝叶斯的强大之处，能够为那么多令人震惊的结论提供解释和方法论的支持。

在 James-Stein estimator 中，我们惊叹于这样一个看上去有着令人着迷性质的频率学派的工具、到后来竟然可以用经验贝叶斯的理论来解释和理解。与之类似地，在本文中我将分享经验贝叶斯发展历程中的另一个令人称奇的故事，它同样之前也是在频率学派中首先被使用、但后来被经验贝叶斯完美的解释，它的名字就是：False Discovery Rate（FDR）。

1995 年 Benjamini和 Hochberg 首次提出了对 Family Wise Error Rate（FWER）进行改进、用于有效控制 FDR 的经典算法，其流程极为简单、但功效强大，感兴趣的伙伴可以参阅原文：

Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal statistical society: series B (Methodological),57(1), 289-300.

要想讲清楚 FDR 与经验贝叶斯的联系，故事得先从 large-scale hypothesis testing 开始说起。

5.2 Large-scale hypothesis testing & FDR

以传统的假设检验为例，设传统假设检验中有 M 个 null hypothesis，一般来说 M 的值通常较小，如：M<20。然而在现代的各种丰富的高维数据场景下，大规模假设检验变成了必须要解决的一个问题。与传统假设检验不同的是，大规模假设检验通常 M 的值较大，如在许多生物信息学的研究中涉及基因组相关的假设检验 M 的值动辄可达几十万乃至上百万。

与小规模的假设检验相比，大规模假设检验的主要问题就在于：如何增加 statistical power （如减少 false negative）？如何控制第一类误差（false postive）？false positive rates 是统计和统计学习中很基础的概念，可以简单写为，其也被称为 type I err r。与之相对应的概念则是 false negative rates，也被称为 type II error。

假设我们有 M 个相互独立的测试，每一个测试都提供了一个反对它的 null hypothesis 的 p 值，这样一来 M 个 p 值就组成了一个向量。令作为一个以为输入、以接受或拒绝 null hypothesis 的决策为输出的规则。假设我们观察到该规则拒绝了 M 个测试中的 R 个，那么 false discovery proportion（Fdp）可以定义为：

，其中是个拒绝的测试中 null hypothesis 实际是正确的数量，换句话说，也就是 false discoveries 的数量。

基于此，规则的 false discovery rate 就被定义为 Fdp 的期望值，写成：，当然，时 Fdp 的值则为 0。

上面提到，我们希望在大规模假设检验中控制 false positive，可怎样的方法能够控制呢？乍一看会觉得没头绪，Benjamini 和 Hochberg 给出他们的经典答案：BH procedures。

5.3 Benjamini-Hochberg (BH) Procedures & Example

BH procedures 能达到什么样的效果呢？预设一个较小的 q 值（如 0.1），BH procedures 可以保证把 Fdr 控制在的水平之内。乍一看觉得很无厘头，到底是怎么做到的呢？事实上很简单，仅以下两步：

把 M 个测试的 p 值从小到大排序：。
令为满足的最大的：对于所有的规则拒绝其 null hypothesis ；对于所有的规则接受其 null hypothesis ；

关于 Fdr 有一个定理：规则的 Fdr 为，其中，是正确的 null hypothesis 的数量。本身是不可观测到的，但在大规模假设检验中，，所以在任何情况下上述定理通常可以描述为规则将 Fdr 水平控制在 q 水平处。

举个很简单的例子：

在关于某种疾病的生物信息学研究中，我们希望比较 50 名患者和 52 名非患者之间 6033 个基因的活动区别，针对每一个基因的活动区别进行一次假设检验，这样就一共产生了 6033 个 z 值，null hypothesis 是 z 值服从标准正态分布。下图展示了 6033 个 z 值的分布以及标准正态分布的对比：（其中黑色实线是由 6033 个 z 值的分布拟合而来、红色虚线代表标准正态分布）

理论上来说，如果所有的 null hypothesis 都是正确的话，上图中 6033 个 z 值的 histogram 应该服从标准正态分布，但很明显上图中红色虚线在中间更高一些、在尾部更低一些。

对以上数据右侧的单侧 p 值利用 BH procedures，设置 q=0.1 我们可以作出下图：

其中黑色点代表从小到大排列的 p 值，红色实线代表的是时的的值（图中 z 值总数用 N 表示）。很明显可以看到，这些 p 值在第 28 个后即穿过了那条的控制线，所以根据 BH procedures，对于拒绝 null hypothesis；对于，我们接受 null hypothesis。利用 BH procedures 所获得的拒绝阈值是。

值得一提的是，如果针对这个问题我们采用的是平时比较熟悉的 Bonferroni bound（），我们获得的拒绝阈值就是，这样的话就只拒绝前4个基因了，换句话说能得到的发现就非常有限了。所以相对来说，Bonferroni bound 太过保守，而 BH procedures 正好能帮助我们在大规模假设检验中有效控制 FDR 的同时有更多的发现。

说到这里，利用 BH procedures 来控制 FDR 怎么看都是很频率学派风格的工具，它究竟是怎么跟经验贝叶斯所联系在一起的呢？

5.4 FDR & Empirical Bayes

事实上，通过 BH procedures 控制 FDR 可以用一个很简单的贝叶斯模型来解释。

在大规模假设检验中，我们假设每一个 z 值对应的假设都是被随机选择成为“null”和“non-null”的。其中被选择成为“null”的先验概率为，而被选择成为“non-null”的概率为。“null”的 z 值们服从正确的 null distribution：标准正态分布，记其 pdf 为，与之相对地“non-null”的 pdf 为。

令和为对应的右侧 cdf，则可写为：以及，且。

则根据贝叶斯定理，我们可以把定义为观测到的 z 值超过 z 时其为 null 的后验概率：

尽管已知，也约等于 1，但这个边际 cdf 仍然是未知的。然而，有一个很明显的非参数式的估计：，代入上式，即可得到：

在实际应用中，。

如果，我们可以决定拒绝 null hypothesis，而根据上式，

，换句话说，对于 null hypothesis ，当时我们拒绝 null hypothesis。

很明显，就是对应的 p 值，，代入上式易得：！这不就是 BH procedures 中看似与贝叶斯毫无关系的操作吗？

不难看出，通过 BH procedures 对 FDR 进行控制可以算是经验贝叶斯的另外一个经典故事：在推导过程中巧妙地利用了观测数据来得到边际 cdf 的估计值。归根结底来说，就是通过从观测数据中估计 null distribution。

5.5 总结

本章节由大规模假设检验入手，介绍了经验贝叶斯中的另外一个成功的故事：BH procedures 控制下的 FDR，这样一个看起来完全频率学派的工具竟然可以用经验贝叶斯来解释和推导，不禁感叹其美妙。

当然，比起频率学派中那种准确的 FDR 控制，从经验贝叶斯的角度出发控制 FDR 采用了一些估计方法并不那么准确，但也打开了新的想象空间和更宽广的可能性。有一个重要的区别在于频率学派的 FDR 控制要求观测值之间是独立的，而经验贝叶斯的结果就不需要。事实上，经验贝叶斯在 FDR 中还有另一个迷人的地方本文尚未提及：local FDR，感兴趣的伙伴可以自行查阅了解。

更多阅读