认识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−1dxBmn01xm11xn1dx

之后，他发现，如果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,...,XnX1X2Xn;
3. 问第7大的数的概率分布?

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

我们先将之一般化，对于一般的情况XkXk的概率密度是什么呢？下面，我们尝试计算一下XkXk落在一个区间[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−kPxXkxΔxn1Δxn1k1xk11xΔxnk

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

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−kfxx0PxXkxΔxΔxn1n1k1xk11xnknk1nkxk11xnkΓn1ΓkΓnk1xk11xnk

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

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

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

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,...,XnX1X2Xn,我们要猜测第k大的数p=XkpXk;
3. 我们再生成m个数，Y1,Y2,...,Ym∼Uniform(0,1)Y1Y2YmUniform01, 其中有mimi个数比p小，m2m2比p大；
4. 求P(p‖Y1,Y2,...,Ym)PpY1Y2Ym的分布是什么。

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

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

对应到Beta分布，后验数据其实相当于是做了m次Bernoulli实验，其中m1m1次比p小，m2m2次比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=1Beta11Γ11Γ1Γ1x111x111

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

贝叶斯学派和频率学派的不同
>假设有一个不均匀的硬币，抛出正面的概率为p，抛掷mm次后，出现正面和翻面的次数分别为m1m1和m2m2，那么按照传统频率学派的观点，p的估计值应该为p̂ =m1mpm1m,而从贝叶斯学派的观点来看，开始对硬币不均匀性一无所知，所以应该假设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,...,XnX1X2Xn;
3. 求(Xk1,Xk2)Xk1Xk2的联合分布。

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

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

P(Xk1∈(x1,x1+Δx),Xk1+k2∈(x2,x2+Δx))=(n1)(n−11)(n−2k1−1,k2−1)xk1−11xk2−12xn−k1−k23(Δx)2PXk1x1x1Δ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−k23fx1x2x3nk11k21nk1k2x1k11x2k21x3nk1k2Gamman1Γ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−13fx1x2x3Γα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−1kDirpαΓ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)β−1dtEp01tBetatαβdt01tΓαβΓαΓβtα11tβ1dtΓαβΓαΓβ01tα1tβ1dt

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

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

而且我们有

∫10Beta(t‖α+1,β)dt=101Betatα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)