贝叶斯推断--Gibbs Sampling

简介

    当我们用贝叶斯模型去建模复杂的问题时，就越来越需要有效的推断方法来求解这些模型了。一言以蔽之，贝叶斯推断的目的就是在一组随即变量上维护一组完全后验概率分布。维护和应用这个概率分布常包含积分步骤，在面对较为复杂的模型时，是无法有效计算的。基于MCMC的抽样技术是一种可能的方法来求解和推断这些模型。

    MCMC抽样的基本思想是，通过均值遍历来估计任意期望分布。也就是说，只要我们能够得到一个后验的足够样本，就可以计算这个的任意统计量值。如下：

E|f(s)|p≃1N∑i=1Nf(s(i))

    其中 P 就是我们讨论的后验概率分布，f(s) 是要求的期望值，也就是目标统计量。f(s(i)) 是从后验分布中抽得的 第 i 个(模拟)样本。

    那么如何从后验分布得到所需的样本呢？Gibbs sampling就是适合这种需求的一种MCMC技术。Gibbs sampling的基本思想是，在每一个变量之间进行不断切换来得到他们的条件分布，在得到其中某一个变量的条件分布的时候，需要保持其他变量的值固定在当前抽样值。举例来说，假定有三个随即变量 X1，X2 和 X3 ，首先对它们进行初始化赋值为 x(0)1 ， x(0)2 和 x(0)3 (即由先验分布抽得的值，在LDA中就是由具有参数 λ 的dirichilet分布抽样值)。那么条件分布的抽取过程就是，在第 i 次迭代中，依次抽取

x(i)1∼p(X1=x1|X2=x(i−1)2,X3=x(i−1)3) ，
x(i)2∼p(X2=x2|X1=x(i−1)1,X3=x(i−1)3) ，
x(i)3∼p(X3=x3|X1=x(i−1)1,X2=x(i−1)2)

    重复上述过程直至算法收敛。这样得到的(模拟)样本值就好象真的是从后验分布抽得一样。算法一给出了一个上述问题的抽样过程：

算法 1 Gibbs sampler

  Initialize x(0)∼q(x)
  for iteration i=1,2,… do
    x(i)1∼p(X1=x1|X2=x(i−1)2,...,XD=x(i−1)D) ，
    x(i)2∼p(X2=x2|X1=x(i−1)1,X3=x(i−1)3,...,XD=x(i−1)D) ，
    x(i)D∼p(XD=xD|X1=x(i−1)1,X2=x(i−1)2,..,,...,XD−1=x(i−1)D−1)

    算法1没有直接从后验分布去抽样，而是模拟了从后验条件分布中进行抽样，一次抽取一个变量。由于初始化的时候随即变量的初值是随机赋予的，因而在早期迭代的抽样值是不准确的，无法对目标后验分布进行充分表达。好在MCMC的理论保证了在算法1中抽样结果收敛之后得到的在抽样样本上的稳定分布就是目标后验联合分布。因而Gibbs sampling一般会跑很多循环以保证抽样达到稳定。

接下来给出一个MCMC技术求解实例。该实例是用Gibbs sampler求解“变换点”模型。

应用实例

    假设我们观测到一个自然数序列 x1,x2,...xN ，其中前1,...,n 个具有一个相同均值,后 n+1,...N 个数具有另外一个相同均值. 自然地，我们用泊松变量表示序列中的自然数在时间 i 的值，并具有下述密度函数：

Poisson(x;λ)=e−λλxx!=exp(xlogλ−λ−log(x!))

其中λ 是分布的均值。进而我们令λ 是一个服从Gamma分布的随机变量，并具有密度函数：

Gamma(λ;a,b)=1Γbaλa−1exp(−bλ)=exp((a−1)logλ−bλ−logΓ(a)+alogb)

初始值λ1 (也就是前面说的第一个均值)经过一个随机n 步后，变为λ2 。因而可以得到下述生成过程。

n∼Uniform(1,2,...,N)
λi∼Gamma(λi;a,b)
x∼{Poisson(xi;λ1),1≤i≤nPoisson(xi;λ2)n<i≤N

    注意到，上述模型具有随机变量n , λ1,2 和 x1,...,N ，其中x1,...,N 是显式变量，其余为隐含变量。由于显式变量是可观测到的，或直接计算得到，所以感兴趣的变量为隐含变量(latent variable)。由贝叶斯理论可以得到：

p(λ1,λ2,n|x1:N)∝p(x1:N|λ1)p(xn+1:N|λ2)p(λ1)p(λ2)p(n)

条件分布

    如算法1中，我们需要得到每个变量的后验条件概率来进行Gibbs抽样。下面介绍得到这些后验条件概率的过程，这个过程首先是得到所有隐含变量的联合概率分布(见上式)，进而推导出每个隐含变量( λ1,2 和 n)的后验条件概率分布。

    基于通过上式中得到的联合概率分布公式，可以得到：

p(x1:N|λ1)p(xn+1:N|λ2)p(λ1)p(λ2)p(n)=(∏i=1np(xi|λ1))(∏i=n+1Np(xi|λ2))p(λ1)p(λ2)p(n)

    对上式取log，就可以得到下式：

logp(x1:N|λ1)+logp(xn+1:N)+logp(λ1)+logp(λ2)+logp(n)=∑ni=1(xilogλ1−λ1−log(x!i))+∑Ni=n+1(xilogλ2−λ2−log(x!i))+(a−1)logλ1−bλ1−logΓ(a)+alogb+(a−1)logλ2−bλ2−logΓ(a)+alogb−logN

     接下来就分别得到各个隐含变量的调降概率公式。其实方法非常简单，举例来说，如果要得到λ1 的条件概率公式，则只需要抽取出全概率公式中与变量λ1 相关的项，加和在一起即可。因而有：

logp(λ1|n,λ2,x1:N)=+∑ni=1(xilogλ1−λ1)+(a−1)logλ1−bλ1=(a+∑ni=1xi−1)logλ1−(a+b)λ1=+logΓ(a+∑ni=1xi,n+b)

    这里值得注意的的是，上式得到的后验概率与先验概率的公式形式的是相同的，即λi∼Γ(λi;a,b)，这是由于Gamma-Poisson是一对共轭分布造成的。这种trick在统计模型中经常用到，用以简化模型推导。