LDA-math-认识Beta/Dirichlet分布

2. 认识Beta/Dirichlet分布
2.1 魔鬼的游戏—认识Beta 分布

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

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

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

潘多拉魔盒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 

E={ X 1 ∈[x,x+Δx],X i ∈[0,x)(i=2,⋯,k),X j ∈(x+Δx,1](j=k+1,⋯,n)}  

事件 E 

则有

P(E) =∏ i=1 n P(X i )=x k−1 (1−x−Δx) n−k Δx=x k−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]  ，

E ′ ={ X 1 ,X 2 ∈[x,x+Δx],X i ∈[0,x)(i=3,⋯,k),X j ∈(x+Δx,1](j=k+1,⋯,n)}  

事件E’

则有

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

从以上分析我们很容易看出，只要落在 [x,x+Δx]  内的数字超过一个，则对应的事件的概率就是 o(Δx)  。于是

 P(x≤X (k) ≤x+Δx)=n(n−1k−1)P(E)+o(Δx)=n(n−1k−1)x k−1 (1−x) n−k Δx+o(Δx)  

所以，可以得到 X (k)   的概率密度函数为

f(x) =lim Δx→0 P(x≤X (k) ≤x+Δx)Δx =n(n−1k−1)x k−1 (1−x) n−k =n!(k−1)!(n−k)! x k−1 (1−x) n−k x∈[0,1]  

利用Gamma 函数，我们可以把 f(x)  表达为

f(x)=Γ(n+1)Γ(k)Γ(n−k+1) x k−1 (1−x) n−k  

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

f(x)=Γ(α+β)Γ(α)Γ(β) x α−1 (1−x) β−1  (1) 

这个就是一般意义上的 Beta 分布！可以证明，在 α,β  取非负实数的时候，这个概率密度函数也都是良定义的。

好，我们回到魔鬼的游戏，这 n=10,k=7  这个具体的实例中，我们按照如下密度分布的峰值去猜测才是最有把握的。

f(x)=10!(6)!(3)! x 6 (1−x) 3 x∈[0,1] 

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

2.2 Beta-Binomial 共轭

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

1. X 1 ,X 2 ,⋯,X n ∼ iid Uniform(0,1)  ，对应的顺序统计量是 X (1) ,X (2) ，⋯,X (n)   , 我们要猜测 p=X (k)   ；
2. Y 1 ,Y 2 ,⋯,Y m ∼ iid Uniform(0,1)  , Y i   中有 m 1   个比 p  小， m 2   个比 p  大；
3. 问 P(p|Y 1 ,Y 2 ,⋯,Y m )  的分布是什么。

由于 p=X (k)   在 X 1 ,X 2 ,⋯,X n   中是第 k  大的，利用 Y i   的信息，我们容易推理得到 p=X (k)   在 X 1 ,X 2 ,⋯,X n ,Y 1 ,Y 2 ,⋯,Y m ∼ iid Uniform(0,1)  这 (m+n)  个独立随机变量中是第 k+m 1   大的，于是按照上一个小节的推理，此时 p=X (k)   的概率密度函数是 Beta(p|k+m 1 ,n−k+1+m 2 )  。按照贝叶斯推理的逻辑，我们把以上过程整理如下：

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

贝努利实验

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

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

以上贝叶斯分析过程的简单直观的表述就是

Beta(p|k,n−k+1)+Count(m 1 ,m 2 )=Beta(p|k+m 1 ,n−k+1+m 2 ) 

其中 (m 1 ,m 2 )  对应的是二项分布 B(m 1 +m 2 ,p)  的计数。更一般的，对于非负实数 α,β  ，我们有如下关系

Beta(p|α,β)+Count(m 1 ,m 2 )=Beta(p|α+m 1 ,β+m 2 ) (2) 

这个式子实际上描述的就是   Beta-Binomial 共轭，此处共轭的意思就是，数据符合二项分布的时候，参数的先验分布和后验分布都能保持Beta 分布的形式，这种形式不变的好处是，我们能够在先验分布中赋予参数很明确的物理意义，这个物理意义可以延续到后验分布中进行解释，同时从先验变换到后验过程中从数据中补充的知识也容易有物理解释。

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

Beta(p|1,1)+Count(α−1,β−1)=Beta(p|α,β)(∗∗∗) 

其中 Beta(p|1,1)  恰好就是均匀分布Uniform(0,1)。

对于(***) 式，我们其实也可以纯粹从贝叶斯的角度来进行推导和理解。 假设有一个不均匀的硬币抛出正面的概率为 p  ,抛 m  次后出现正面和反面的次数分别是 m 1 ,m 2   ，那么按传统的频率学派观点， p  的估计值应该为 p ^ =m 1 m   。而从贝叶斯学派的观点来看，开始对硬币不均匀性一无所知，所以应该假设 p∼Uniform(0,1)  , 于是有了二项分布的计数 (m 1 ,m 2 )  之后，按照贝叶斯公式如下计算 p  的后验分布

P(p|m 1 ,m 2 ) =P(p)⋅P(m 1 ,m 2 |p)P(m 1 ,m 2 ) =1⋅P(m 1 ,m 2 |p)∫ 1 0 P(m 1 ,m 2 |t)dt =(mm 1 )p m 1  (1−p) m 2  ∫ 1 0 (mm 1 )t m 1  (1−t) m 2  dt =p m 1  (1−p) m 2  ∫ 1 0 t m 1  (1−t) m 2  dt   

计算得到的后验分布正好是 Beta(p|m 1 +1,m 2 +1)  。

百变星君Beta分布

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

在上一个小节中，我们从二项分布推导Gamma 分布的时候，使用了如下的等式

P(C≤k)=n!k!(n−k−1)! ∫ 1 p t k (1−t) n−k−1 dt,C∼B(n,p) (3) 

现在大家可以看到，左边是二项分布的概率累积，右边实际上是 Beta(t|k+1,n−k)  分布的概率积分。这个式子在上一小节中并没有给出证明，下面我们利用和魔鬼的游戏类似的概率物理过程进行证明。

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

贝努利实验最多成功 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)=∫ 1 p Beta(t|k+1,n−k)dt=n!k!(n−k−1)! ∫ 1 p t k (1−t) n−k−1 dt  

最后我们再回到魔鬼的游戏，如果你按出的5个随机数字中，魔鬼告诉你有2个小于它手中第7大的数，那么你应该
按照如下概率分布的峰值做猜测是最好的

Beta(x|9,7)=15!(8)!(6)! x 8 (1−x) 6 x∈[0,1] 

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

2.3 Dirichlet-Multinomial 共轭

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

1. X 1 ,X 2 ,⋯,X n ∼ iid Uniform(0,1)  ，
2. 排序后对应的顺序统计量 X (1) ,X (2) ，⋯,X (n)   ,
3. 问 (X (k 1 ) ,X (k 1 +k 2 ) )  的联合分布是什么；

游戏3

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

(X (k 1 ) ,X (k 1 +k 2 ) )  的联合分布推导

 P(X (k 1 ) ∈(x 1 ,x 1 +Δx),X (k 1 +k 2 ) ∈(x 2 ,x 2 +Δx))=n(n−1)(n−2k 1 −1,k 2 −1)x k 1 −1 1 x k 2 −1 2 x n−k 1 −k 2  3 (Δx) 2 =n!(k 1 −1)!(k 2 −1)!(n−k 1 −k 2 )! x k 1 −1 1 x k 2 −1 2 x n−k 1 −k 2  3 (Δx) 2   

于是我们得到 (X (k 1 ) ,X (k 1 +k 2 ) )  的联合分布是

f(x 1 ,x 2 ,x 3 ) =n!(k 1 −1)!(k 2 −1)!(n−k 1 −k 2 )! x k 1 −1 1 x k 2 −1 2 x n−k 1 −k 2  3 =Γ(n+1)Γ(k 1 )Γ(k 2 )Γ(n−k 1 −k 2 +1) x k 1 −1 1 x k 2 −1 2 x n−k 1 −k 2  3   

熟悉 Dirichlet的同学一眼就可以看出，上面这个分布其实就是3维形式的 Dirichlet 分布 Dir(x 1 ,x 2 ,x 3 |k 1 ,k 2 ,n−k 1 −k 2 +1)  。令 α 1 =k 1 ,α 2 =k 2 ,α 3 =n−k 1 −k 2 +1  ,于是分布密度可以写为

f(x 1 ,x 2 ,x 3 )=Γ(α 1 +α 2 +α 3 )Γ(α 1 )Γ(α 2 )Γ(α 3 ) x α 1 −1 1 x α 2 −1 2 x α 3 −1 3  (4) 

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

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

不同 α  下的Dirichlet 分布

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

1.   X 1 ,X 2 ,⋯,X n ∼ iid Uniform(0,1)  ，排序后对应的顺序统计量 X (1) ,X (2) ，⋯,X (n)  
2. 令 p 1 =X (k 1 ) ,p 2 =X (k 1 +k 2 ) ,p 3 =1−p 1 −p 2   (加上 p 3   是为了数学表达简洁对称),我们要猜测 p → =(p 1 ,p 2 ,p 3 )  ；
3. Y 1 ,Y 2 ,⋯,Y m ∼ iid Uniform(0,1)  , Y i   中落到 [0,p 1 ),[p 1 ,p 2 ),[p 2 ,1]  三个区间的个数分别为 m 1 ,m 2 ,m 3   ， m=m 1 +m 2 +m3  ；
4.  问后验分布 P(p → |Y 1 ,Y 2 ,⋯,Y m )  的分布是什么。

游戏4

为了方便，我们记

m − →  =(m 1 ,m 2 ,m 3 ),k → =(k 1 ,k 2 ,n−k 1 −k 2 +1) 

由游戏中的信息，我们可以推理得到 p 1 ,p 2   在 X 1 ,X 2 ,⋯,X n ,  Y 1 ,Y 2 ,⋯,Y m   ∼ iid Uniform(0,1)  这 m+n  个数中分别成为了第 k 1 +m 1 ,k 2 +m 2   大的数，于是后验分布 P(p → |Y 1 ,Y 2 ,⋯,Y m )  应该是 Dir(p → |k 1 +m 1 ,k 1 +m 2 ,n−k 1 −k 2 +1+m 3 )  ,即 Dir(p → |k → +m − →  )  。按照贝叶斯推理的逻辑，我们同样可以把以上过程整理如下：

1.  我们要猜测参数 p → =(p 1 ,p 2 ,p 3 )  ，其先验分布为 Dir(p → |k → )  ；
2.  数据 Y i   落到 [0,p 1 ),[p 1 ,p 2 ),[p 2 ,1]  三个区间的个数分别为 m 1 ,m 2 ,m 3   ，所以 m − →  =(m 1 ,m 2 ,m 3 )  服从多项分布 Mult(m − →  |p → ) 
3.  在给定了来自数据提供的知识 m − →    后， p →   的后验分布变为 Dir(p → |k → +m − →  ) 

贝叶斯推理过程

以上贝叶斯分析过程的简单直观的表述就是

Dir(p → |k → )+MultCount(m − →  )=Dir(p → |k → +m − →  ) 

令 α → =k →   ,把 α →   从整数集合延拓到实数集合，更一般的可以证明有如下关系

Dir(p → |α → )+MultCount(m − →  )=Dir(p|α → +m − →  ) (5) 

以上式子实际上描述的就是  Dirichlet-Multinomial 共轭，而我们从以上过程可以看到，Dirichlet 分布中的参数 α →   都可以理解为物理计数。类似于 Beta 分布，我们也可以把 Dir(p → |α → )  作如下分解

Dir(p → |1 → )+MultCount(m − →  −1 → )=Dir(p → |α → ) 

此处 1 → =(1,1,⋯,1)  。自然，上式我们也可以类似地用纯粹贝叶斯的观点进行推导和解释。

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

Dir(p → |α → )=Γ(∑ K k=1 α k )∏ K k=1 Γ(α k ) ∏ k=1 K p α k −1 k  (6) 

对于给定的 p →   和 N  ,多项分布定义为

Mult(n → |p → ,N)=(Nn → )∏ k=1 K p n k  k  (7) 

而 Mult(n → |p → ,N)  和 Dir(p → |α → )  这两个分布是共轭关系。

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

2.4 Beta/Dirichlet 分布的一个性质

如果 p∼Beta(t|α,β)  , 则

E(p) =∫ 1 0 t∗Beta(t|α,β)dt=∫ 1 0 t∗Γ(α+β)Γ(α)Γ(β) t α−1 (1−t) β−1 dt=Γ(α+β)Γ(α)Γ(β) ∫ 1 0 t α (1−t) β−1 dt  

上式右边的积分对应到概率分布 Beta(t|α+1,β)  ，对于这个分布，我们有

∫ 1 0 Γ(α+β+1)Γ(α+1)Γ(β) t α (1−t) β−1 dt=1 

把上式带入 E(p)  的计算式，得到

E(p) =Γ(α+β)Γ(α)Γ(β) ⋅Γ(α+1)Γ(β)Γ(α+β+1) =Γ(α+β)Γ(α+β+1) Γ(α+1)Γ(α) =αα+β  (8) 

这说明，对于Beta 分布的随机变量，其均值可以用 αα+β   来估计。Dirichlet 分布也有类似的结论，如果 p → ∼Dir(t → |α → )  ，同样可以证明

E(p → )=(α 1 ∑ K i=1 α i  ,α 2 ∑ K i=1 α i  ,⋯,α K ∑ K i=1 α i  ) (9) 

以上两个结论很重要，因为我们在后面的 LDA 数学推导中需要使用这个结论。