认识beta和dirichlet分布

认识Beta/Dirichlet分布

本文主要是对rickjin的《LDA数学八卦》以及PRML一书中关于Beta分布和Dirichlet分布的知识整理。

我们知道Gamma函数 

Γ(x)=∫∞0tx−1e−tdt Γx0∞tx1etdt

是阶乘运算在实数集上延伸。它具有如下性质

Γ(x+1)=xΓ(x) Γx1xΓx

所以，我们其实有 Γ(n)=(n−1)! Γnn1

是不是感觉有些奇怪？为什么不是  Γ(n)=n! Γnn 而是  Γ(n)=(n−1)! Γnn1？  
今天读了rickjin的《LDA数学八卦》才知道，原来欧拉是研究了Beta函数 

B(m,n)=∫10xm−1(1−x)n−1dx Bmn01xm11xn1dx

之后，他发现，如果Gamma函数的定义选取满足 Γ(n)=(n−1)! Γnn1,那么Beta函数会有一个很漂亮的对称形式 

B(m,n)=Γ(m)Γ(n)Γ(m+n) BmnΓmΓnΓmn

而如果选取 Γ(n)=n! Γnn的定义，则有 

B′(m,n)=Γ(m)Γ(n)Γ(m+n+1) BmnΓmΓnΓmn1

这个形式显然不如B(m,n)优美，而数学家总是很在乎美感的。

讲完了Gamma函数，我们再来扒一扒Beta函数，他又有什么物理意义呢？

Beta分布

我之前找工作面试的时候，曾经被面过这么一个题目

1. X∼Uniform(0,1) XUniform01;
2. 随机生成10个数，把这10个数排序后得到的顺序统计量是 X1,X2,...,Xn X1X2Xn;
3. 问第7大的数的概率分布?

那时候我是不知道Beta分布，否则肯定不会被虐的这么惨::>_<:: 

我们先将之一般化，对于一般的情况 Xk Xk的概率密度是什么呢？下面，我们尝试计算一下 Xk Xk落在一个区间 [x,x+Δx] xxΔx的概率值

P(x≤Xk≤x+Δx)=? PxXkxΔx

如上图所示，我们把[0,1]区间分成三段 [0,x),[x,x+Δx],(x+Δx,1] 0xxxΔxxΔx1三段。我们假定， Δx Δx足够小，只能够容纳一个点,则由排列组合理论可得

P(x≤Xk≤x+Δx)=(n1)Δx(n−1k−1)xk−1(1−x−Δx)n−k PxXkxΔxn1Δxn1k1xk11xΔxnk

所以我们可以得到 Xk Xk的概率密度函数为

f(x)=limx→0P(x≤Xk≤x+Δx)Δx=(n1)(n−1k−1)xk−1(1−x)n−k=n!(k−1)!(n−k)!xk−1(1−x)n−k=Γ(n+1)Γ(k)Γ(n−k+1)xk−1(1−x)n−k fxx0PxXkxΔxΔxn1n1k1xk11xnknk1nkxk11xnkΓn1ΓkΓnk1xk11xnk

我们取 α=k,β=n−k+1 αkβnk1,于是

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

这就是Beta分布！
回到上面那个面试题，把 n=10,k=7 n10k7带入其中，得到密度函数

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

Beta-Binomial共轭

上边的面试题还有几个衍化版本，我们先看第一个衍化版本：

1. X∼Uniform(0,1) XUniform01;
2. 随机生成n个数,由小到大排序后为 X1,X2,...,Xn X1X2Xn,我们要猜测第k大的数 p=Xk pXk;
3. 我们再生成m个数， Y1,Y2,...,Ym∼Uniform(0,1) Y1Y2YmUniform01, 其中有 mi mi个数比p小， m2 m2比p大；
4. 求 P(p‖Y1,Y2,...,Ym) PpY1Y2Ym的分布是什么。

容易看出，我们一共生成了 m+n mn个数，而 p=Xk pXk在最终生成的m+n个数中，是第 k+m1 km1大的。按照我们之前讲过的Beta分布的逻辑，p其实应该服从 α=k+m1,β=n−k+1+m2 αkm1βnk1m2的Beta分布。我们知道贝叶斯学派进行参数估计的基本过程是

先验分布 + 后验数据 = 后验分布

对应到Beta分布，后验数据其实相当于是做了m次Bernoulli实验，其中 m1 m1次比p小， m2 m2次比p大，相当于

Beta(p‖α,β)+BinomCount(m1,m2)=Beta(p‖α+m1,β+m2) BetapαβBinomCountm1m2Betapαm1βm2

上面这个式子描述的就是 Beta−Bonomial共轭 BetaBonomial.

共轭的意思是，参数的先验分布和后验分布都能保持Beta分布的形式，这样的好处是，我们能够在先验分布中赋予参数明确的物理意义，并且这个物理意义可以通过后验数据，延续到后验分布中进行解释。

由上边的解释可知，Beta分布重的参数 α,β αβ其实都可以理解为物理技术，这两个参数也经常被称为伪计数（pseudo-count)。所以， Beta(α,β) Betaαβ也可以理解为

Beta(α,β)=Beta(1,1)+BinomCount(α−1,β−1) BetaαβBeta11BinomCountα1β1

其中

Beta(1,1)=Γ(1+1)Γ(1)Γ(1)x1−1(1−x)1−1=1 Beta11Γ11Γ1Γ1x111x111

这恰好就是均匀分布 Uniform(0,1) Uniform01。

贝叶斯学派和频率学派的不同
>假设有一个不均匀的硬币，抛出正面的概率为p，抛掷 m m次后，出现正面和翻面的次数分别为 m1 m1和 m2 m2，那么按照传统频率学派的观点，p的估计值应该为 p̂ =m1m pm1m,而从贝叶斯学派的观点来看，开始对硬币不均匀性一无所知，所以应该假设p服从均匀分布 Uniform(0,1) Uniform01,也就是 Beta(1,1) Beta11,于是在有了后验数据之后，我们得出p其实应该服从 Beta(p|m1+1,m2+1) Betapm11m21.

百变星君Beta分布

Beta分布的概率密度如上图， α,β αβ的不同，他可以是凹的、凸的、单调上升的、单调下降的，可以是曲线也可以是直线； 而且，如前所述，均匀分布也是Beta分布的一种特殊形式。正是由于Beta分布能够你和如何之多的形状，因此他经常被贝叶斯学派用作先验分布。

Dirichlet-Multinomial共轭：Beta分布的高维推广

更一步的问题

1. X∼Uniform(0,1) XUniform01;
2. 随机生成n个数，排序后为 X1,X2,...,Xn X1X2Xn;
3. 求 (Xk1,Xk2) Xk1Xk2的联合分布。

同推导Beta分布类似,我们取 Δx Δx足够小，只能容纳一个点.

由于 Δx Δx足够小，我们有 x1+x2+x3=1 x1x2x31. 

P(Xk1∈(x1,x1+Δx),Xk1+k2∈(x2,x2+Δx))=(n1)(n−11)(n−2k1−1,k2−1)xk1−11xk2−12xn−k1−k23(Δx)2 PXk1x1x1ΔxXk1k2x2x2Δxn1n11n2k11k21x1k11x2k21x3nk1k2Δx2

于是我们得到 (Xk1,Xk2) Xk1Xk2的联合分布为

f(x1,x2,x3)=n!(k1−1)!(k2−1)!(n−k1−k2)!xk1−11x2k2−1xn−k1−k23=Gamma(n+1)Γ(k1)Γ(k2)Γ(n−k1−k2+1)xk1−11xk2−12xn−k1−k23 fx1x2x3nk11k21nk1k2x1k11x2k21x3nk1k2Gamman1Γk1Γk2Γnk1k21x1k11x2k21x3nk1k2

令 α1=k1,α2=k2,α3=n−k1−k2+1 α1k1α2k2α3nk1k21,我们得到

f(x1,x2,x3)=Γ(α1+α2+α3)Γ(α1)Γ(α2)Γ(α3)xα1−11xα2−12xα3−13 fx1x2x3Γα1α2α3Γα1Γα2Γα3x1α11x2α21x3α31

上边这个分布其实就是一个三维形式的Dirichlet分布 Dir(α1,α2,α3) Dirα1α2α3.同Beta分布类似，Dirichlet分布也是一个百变星君，下图为不同 α α值时Dirichlet分布的图像。

一般形式的Dirichlet分布定义如下

Dir(p⃗ ‖α⃗ )=Γ(∑Kk=1αk)∏k=1KΓ(αk)∏k=1Kpαk−1k DirpαΓk1Kαkk1KΓαkk1Kpkαk1

Dirichlet分布也是Binomial共轭的

Dir(p⃗ ‖α⃗ )+MultCount(m⃗ )=Dir(p⃗ ‖α⃗ +m⃗ ) DirpαMultCountmDirpαm

我们同样也有

Dir(p⃗ ‖α⃗ )=Dir(p⃗ ‖1⃗ )+MultCount(m⃗ −1⃗ ) DirpαDirp1MultCountm1

Beta分布和Dirichlet分布的性质

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

E(p)=∫10t∗Beta(t‖α,β)dt=∫10t∗Γ(α+β)Γ(α)Γ(β)tα−1(1−t)β−1dt=Γ(α+β)Γ(α)Γ(β)∫10tα(1−t)β−1dt Ep01tBetatαβdt01tΓαβΓαΓβtα11tβ1dtΓαβΓαΓβ01tα1tβ1dt

上式右边的积分对应到概率分布 Beta(t‖α+1,β) Betatα1β

Beta(t‖α+1,β)=∫10t∗Γ(α+β+1)Γ(α+1)Γ(β)tα(1−t)β−1dt Betatα1β01tΓαβ1Γα1Γβtα1tβ1dt

而且我们有

∫10Beta(t‖α+1,β)dt=1 01Betatα1βdt1

所以我们有

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

把上式带入E(p)中得

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

同样的，对于Dirichlet分布我们可以得到

E(p⃗ )=(α1∑i=1Kαi,α2∑i=1Kαi,...,αK∑i=1Kαi)