(转)认识Beta/Dirichlet分布 #

Date: 2015/6/17

本文来自 我爱自然语言处理的博客我爱自然语言处理的博客more魔鬼的游戏—认识Beta 分布

统计学就是猜测上帝的游戏,当然我们不总是有机会猜测上帝，运气不好的时候就得揣度魔鬼的心思。有一天你被魔鬼撒旦抓走了，撒旦说：”你们人类很聪明，而我是很仁慈的，和你玩一个游戏，赢了就可以走，否则把灵魂出卖给我。游戏的规则很简单，我有一个魔盒，上面有一个按钮，你每按一下按钮，就均匀的输出一个[0,1]之间的随机数，我现在按10下，我手上有10个数，你猜第7大的数是什么，偏离不超过0.01就算对。“ 你应该怎么猜呢？

从数学的角度抽象一下，上面这个游戏其实是在说随机变量X1,X2,⋯,Xniid∼Uniform(0,1)，把这n 个随机变量排序后得到顺序统计量 X(1),X(2)，⋯,X(n), 然后问 X(k) 的分布是什么。

对于不喜欢数学的同学而言，估计每个概率分布都是一个恶魔，那在概率统计学中，均匀分布应该算得上是潘多拉魔盒，几乎所有重要的概率分布都可以从均匀分布Uniform(0,1)中生成出来;尤其是在统计模拟中，所有统计分布的随机样本都是通过均匀分布产生的。

pandora潘多拉魔盒Uniform(0,1)

对于上面的游戏而言 n=10,k=7, 如果我们能求出 X(7) 的分布的概率密度，那么用概率密度的极值点去做猜测就是最好的策略。对于一般的情形，X(k) 的分布是什么呢？那我们尝试计算一下X(k) 落在一个区间 [x,x+Δx] 的概率，也就是求如下概率值
P(x≤X(k)≤x+Δx)=?

把 [0,1] 区间分成三段 [0,x),[x,x+Δx],(x+Δx,1],我们先考虑简单的情形，假设n 个数中只有一个落在了区间 [x,x+Δx]内，则因为这个区间内的数X(k)是第k大的，则[0,x)中应该有 k−1 个数，(x,1] 这个区间中应该有n−k 个数。不失一般性，我们先考虑如下一个符合上述要求的事件E

\begin{align*}  E = { & X_1 \in [x, x+\Delta x], \  & X_i \in [0,x)\quad (i=2,\cdots,k), \  & X_j \in (x+\Delta x,1] \quad (j=k+1,\cdots,n)}  \end{align*}

beta-game-1事件 E

则有
P(E)=n∏i=1P(Xi) =xk−1(1−x−Δx)n−kΔx =xk−1(1−x)n−kΔx+o(Δx)

o(Δx)表示Δx的高阶无穷小。显然，由于不同的排列组合，即n个数中有一个落在 [x,x+Δx]区间的有n种取法，余下n−1个数中有k−1个落在[0,x)的有(n−1k−1)种组合，所以和事件E等价的事件一共有 n(n−1k−1)个。继续考虑稍微复杂一点情形，假设n 个数中有两个数落在了区间 [x,x+Δx]，

\begin{align*}  E’ = { & X_1,X_2\in [x, x+\Delta x], \  & X_i \in [0,x) \quad (i=3,\cdots,k), \  & X_j \in (x+\Delta x,1] \quad (j=k+1,\cdots,n)}  \end{align*}

beta-game-2事件E’

则有
P(E′)=xk−2(1−x−Δx)n−k(Δx)2=o(Δx)

还记得神奇的 Gamma 函数可以把很多数学概念从整数集合延拓到实数集合吧。我们在上式中取α=k,β=n−k+1, 于是我们得到
f(x)=Γ(α+β)Γ(α)Γ(β)xα−1(1−x)β−1

好，我们回到魔鬼的游戏，这n=10,k=7这个具体的实例中，我们按照如下密度分布的峰值去猜测才是最有把握的。
f(x)=10!(6)!(3)!x6(1−x)3x∈[0,1]

然而即便如此，我们能做到一次猜中的概率也不高，很不幸，你第一次没有猜中，魔鬼微笑着说：“我再仁慈一点，再给你一个机会，你按5下这个机器，你就得到了5个[0,1]之间的随机数，然后我可以告诉你这5个数中的每一个，和我的第7大的数相比，谁大谁小，然后你继续猜我手头的第7大的数是多少。”这时候我们应该怎么猜测呢？Beta-Binomial 共轭

魔鬼的第二个题目，数学上形式化一下，就是

X1,X2,⋯,Xniid∼Uniform(0,1)，对应的顺序统计量是 X(1),X(2)，⋯,X(n), 我们要猜测 p=X(k)；
Y1,Y2,⋯,Ymiid∼Uniform(0,1), Yi中有m1个比p小，m2个比p大；
问 P(p|Y1,Y2,⋯,Ym) 的分布是什么。
由于p=X(k)在 X1,X2,⋯,Xn中是第k大的，利用Yi的信息，我们容易推理得到 p=X(k) 在X1,X2,⋯,Xn,Y1,Y2,⋯,Ymiid∼Uniform(0,1) 这(m+n)个独立随机变量中是第 k+m1大的，于是按照上一个小节的推理，此时p=X(k) 的概率密度函数是 Beta(p|k+m1,n−k+1+m2)。按照贝叶斯推理的逻辑，我们把以上过程整理如下：

p=X(k)是我们要猜测的参数，我们推导出 p 的分布为 f(p)=Beta(p|k,n−k+1),称为 p 的先验分布；
数据Yi中有m1个比p小，m2个比p大，Yi相当于是做了m次贝努利实验，所以m1 服从二项分布 B(m,p)；
在给定了来自数据提供的(m1,m2)的知识后，p 的后验分布变为 f(p|m1,m2)=Beta(p|k+m1,n−k+1+m2)

coin-toss贝努利实验

我们知道贝叶斯参数估计的基本过程是

先验分布 + 数据的知识 = 后验分布

以上贝叶斯分析过程的简单直观的表述就是
Beta(p|k,n−k+1)+Count(m1,m2)=Beta(p|k+m1,n−k+1+m2)

而我们从以上过程可以看到，Beta 分布中的参数α,β都可以理解为物理计数，这两个参数经常被称为伪计数(pseudo-count)。基于以上逻辑，我们也可以把Beta(p|α,β)写成下式来理解
Beta(p|1,1)+Count(α−1,β−1)=Beta(p|α,β)(∗∗∗)

对于(*) 式，我们其实也可以纯粹从贝叶斯的角度来进行推导和理解。 假设有一个不均匀的硬币抛出正面的概率为p,抛m次后出现正面和反面的次数分别是m1,m2，那么按传统的频率学派观点，p的估计值应该为 ˆp=m1m。而从贝叶斯学派的观点来看，开始对硬币不均匀性一无所知，所以应该假设p∼Uniform(0,1), 于是有了二项分布的计数(m1,m2)
之后，按照贝叶斯公式如下计算p 的后验分布
$$\begin{align}
P(p|m_1,m_2) & = \frac{P(p)\cdot P(m_1,m_2|p)}{P(m_1,m_2)} \
& = \frac{1\cdot P(m_1,m_2|p)}{\int_0^1 P(m_1,m_2|t)dt} \
& = \frac{\binom{m}{m_1}p^{m_1}(1-p)^{m_2}}{\int_0^1 \binom{m}{m_1}t^{m_1}(1-t)^{m_2}dt} \
& = \frac{p^{m_1}(1-p)^{m_2}}{\int_0^1 t^{m_1}(1-t)^{m_2}dt}
\end{align}$$
计算得到的后验分布正好是 Beta(p|m1+1,m2+1)。

beta-distribution 百变星君Beta分布

Beta 分布的概率密度我们把它画成图，会发现它是个百变星君，它可以是凹的、凸的、单调上升的、单调下降的；可以是曲线也可以是直线，而均匀分布也是特殊的Beta分布。由于Beta 分布能够拟合如此之多的形状，因此它在统计数据拟合中被广泛使用。

在上一个小节中，我们从二项分布推导Gamma 分布的时候，使用了如下的等式
P(C≤k)=n!k!(n−k−1)!∫1ptk(1−t)n−k−1dt,C∼B(n,p)

我们可以如下构造二项分布，取随机变量 X1,X2,⋯,Xniid∼Uniform(0,1),一个成功的贝努利实验就是 Xi<p,否则表示失败,于是成功的概率为p。C用于计数成功的次数，于是C∼B(n,p)。

beta-binomial 贝努利实验最多成功k次

显然我们有如下式子成立

P(C≤k)=P(X(k+1)>p)

此处X(k+1)是顺序统计量，为第k+1大的数。等式左边表示贝努利实验成功次数最多k次，右边表示第 k+1 大的数必然对应于失败的贝努利实验，从而失败次数最少是n−k次，所以左右两边是等价的。由于X(k+1)∼Beta(t|k+1,n−k), 于是
P(C≤k)=P(X(k+1)>p) =∫1pBeta(t|k+1,n−k)dt =n!k!(n−k−1)!∫1ptk(1−t)n−k−1dt

最后我们再回到魔鬼的游戏，如果你按出的5个随机数字中，魔鬼告诉你有2个小于它手中第7大的数，那么你应该
按照如下概率分布的峰值做猜测是最好的
Beta(x|9,7)=15!(8)!(6)!x8(1−x)6x∈[0,1]

很幸运的，你这次猜中了，魔鬼开始甩赖了：这个游戏对你来说太简单了，我要加大点难度，我们重新来一次，我按魔盒20下生成20个随机数，你同时给我猜第7大和第13大的数是什么，这时候应该如何猜测呢？Dirichlet-Multinomial 共轭

对于魔鬼变本加厉的新的游戏规则，数学形式化如下：

X1,X2,⋯,Xniid∼Uniform(0,1)，
排序后对应的顺序统计量 X(1),X(2)，⋯,X(n),
问 (X(k1),X(k1+k2))的联合分布是什么；
游戏3

完全类似于第一个游戏的推导过程，我们可以进行如下的概率计算(为了数学公式的简洁对称，我们取x3满足x1+x2+x3=1,但只有x1,x2是变量)

dirichlet-game
(X(k1),X(k1+k2))的联合分布推导

P(X(k1)∈(x1,x1+Δx),X(k1+k2)∈(x2,x2+Δx)) =n(n−1)(n−2k1−1,k2−1)xk1−11xk2−12xn−k1−k23(Δx)2 =n!(k1−1)!(k2−1)!(n−k1−k2)!xk1−11xk2−12xn−k1−k23(Δx)2

这个就是一般形式的3维 Dirichlet 分布，即便 →α=(α1,α2,α3) 延拓到非负实数集合，以上概率分布也是良定义的。

从形式上我们也能看出，Dirichlet 分布是Beta 分布在高维度上的推广，他和Beta 分布一样也是一个百变星君，密度函数可以展现出多种形态。

dirichlet-distribution不同 α 下的Dirichlet 分布

类似于魔鬼的游戏2，我们也可以调整一下游戏3，从魔盒中生成m个随机数Y1,Y2,⋯,Ymiid∼Uniform(0,1) 并让魔鬼告诉我们Yi和(X(k1),X(k1+k2))相比谁大谁小。于是有如下游戏4

X1,X2,⋯,Xniid∼Uniform(0,1)，排序后对应的顺序统计量 X(1),X(2)，⋯,X(n)
令p1=X(k1),p2=X(k1+k2),p3=1−p1−p2(加上p3是为了数学表达简洁对称),我们要猜测 →p=(p1,p2,p3)；
Y1,Y2,⋯,Ymiid∼Uniform(0,1), Yi中落到[0,p1),[p1,p2),[p2,1] 三个区间的个数分别为 m1,m2,m3，m=m1+m2+m3；
问后验分布 P(→p|Y1,Y2,⋯,Ym) 的分布是什么。
游戏4

为了方便，我们记
→m=(m1,m2,m3),→k=(k1,k2,n−k1−k2+1)

我们要猜测参数 →p=(p1,p2,p3)，其先验分布为Dir(→p|→k)；
数据Yi落到[0,p1),[p1,p2),[p2,1]三个区间的个数分别为 m1,m2,m3，所以→m=(m1,m2,m3) 服从多项分布Mult(→m|→p)
在给定了来自数据提供的知识→m后，→p 的后验分布变为 Dir(→p|→k+→m)
贝叶斯推理过程

以上贝叶斯分析过程的简单直观的表述就是
Dir(→p|→k)+MultCount(→m)=Dir(→p|→k+→m)

以上的游戏我们还可以往更高的维度上继续推，譬如猜测 X(1),X(2)，⋯,X(n) 中的4、5、…等更多个数，于是就得到更高纬度的 Dirichlet 分布和 Dirichlet-Multinomial 共轭。一般形式的 Dirichlet 分布定义如下
Dir(→p|→α)=Γ(∑Kk=1αk)∏Kk=1Γ(αk)K∏k=1pαk−1k

Beta-Binomail 共轭和 Dirichlet-Multinomail 共轭都可以用纯粹数学的方式进行证明，我们在这两个小节中通过一个游戏来解释这两个共轭关系，主要是想说明这个共轭关系是可以对应到很具体的概率物理过程的。

2.4 Beta/Dirichlet 分布的一个性质

如果 p∼Beta(t|α,β), 则
$$\begin{align*}
E(p) & = \int_0^1 tBeta(t|\alpha,\beta)dt \
& = \int_0^1 t \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} t^{\alpha-1}(1-t)^{\beta-1}dt \
& = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \int_0^1 t^{\alpha}(1-t)^{\beta-1}dt
\end{align*}上式右边的积分对应到概率分布$Beta(t|α+1,β)$，对于这个分布，我们有