MCMC采样技术及其在贝叶斯推断上的应用

房爱东，谢士春

(宿州学院信息工程学院，安徽 宿州 234000)

从理论上讲，贝叶斯统计是非常简单的——后验分布正比例于样本似然分布与先验分布.比例因子往往很难得到，所以不能直接用于精确推断.大多数情况下，比例因子要求求解一个数值积分，然而当有多个参数时可能是困难的.尽管统计学家们从他们的决策理论研究中发现贝叶斯统计在理论上提供了真正优势，但在实践中这些优势并不是真的可用，计算边缘后验分布密度函数的困难成为阻碍贝叶斯方法广泛应用的最大障碍.随着计算机算力的增强，计算贝叶斯统计技术的出现改变了这一切——统计推断可以基于后验分布抽取的随机样本，而且从采样样本中计算出的估计量可以通过设置足够大的样本大小来实现任何所需的精度[1].马氏链式蒙特卡罗采样法的基本思想是：构造一个马氏链，并将后验分布作为其稳态分布[2-4].只要我们让马氏链运行足够长，从链中提取的样本可以被认为是从目标(后验)分布中的样本抽取，然后基于这些样本就可以做出各种贝叶斯统计推断.

1 蒙特卡罗类方法背后的思想

引入采样技术的目的是获取指定概率分布(密度)的采样点.大数定理指出，如果能从未知的概率分布p(x)中获取独立同分布的采样点xi，那么随着采集样本数的增加，样本分布将收敛到真正分布[5].图1表明，随着采样点的增多，每个样本点出现的频率渐近于该点的概率.

图1 采样点出现的频率渐近于该点的实际概率(右侧图)

也就是说，样本点x出现的频率PN(x)随着采样次数的增多近似于其概率PN(x)[6].

(1)

在贝叶斯推断中，常常需要计算各种统计量，例如求解函数f(x)关于随机变量x的期望，且p(x):

(2)

(3)

2 独立采样法

利用采样技术的目的是采集指定概率分布p(x)(目标分布)的样本点，以用于近似计算.针对不同形式的概率分布，有不同的解决方案：

(1)如果能够直接从目标分布抽样，例如离散型随机模型，或者连续型密度函数p(x)的累积分布函数存在解析形式的逆函数，借助于类似轮盘赌法或逆采样法即可获得独立同分布样本点.

(2)如果不能，只能求数值解，可以考虑用拒绝/接受或者重要采样法.事先选定一个已知的且易抽样(概率统计教科书中介绍的常见分布)的提议分布q(x)，从中抽取样本点x(i)作为候选样本点，并应用某种判优规则： 若判定结果略优于均匀采样，那么此样本就可以接受为p(x)的样本点，否则拒绝.这种抽样的准确率很大程度上取决于p(x)和q(x)的相似度，若相差很大，拒绝率会很高，特别是高维分布.

图2 拒绝采样法

算法思想：寻找易采样的提议分布q(x)，确定一个常数M，使得Mq(x)完全“裹”住p(x).重复如下步骤，直至获取N个样本点：

(1)采样x(i)～q(x),u～Unif(0,1)；(2)如果u接受x(i)；否则拒绝，转(1).

图2阐释了拒绝采样的基本思想[7]：对于q(x)产生的采样点x(i)，由于p(x(i))≥uMq(x(i))，所以x(i)作为目标分布p(x)的高密度点(High Density Area，高密区)被选中，作为低密点的x(j)却被拒绝.事实上，低密点代表稀有事件，高密点才是较具代表性的样本点.但由于u是均匀分布，一定概率的低密度点同样也可能被选中，以满足充分采样的需要，否则我们采样得到的样本平均值将高于从后验分布中抽取独立样本的均值.这种采样法的缺点一是M较难确定，二是要求提议分布q(x)的目标分布p(x)形状相似(两者的高密区低密区相似)，否则拒绝率较高.

3 马氏链式蒙特卡罗(MCMC)类采样法

MCMC类采样法通过设计转移矩阵，生成一条马尔可夫链使其极限分布等于目标分布.主要包括MH算法和Gibbs算法.MH算法适用于一维随机变量的采样，对于高维分布需要使用Gibbs算法，该方法是前者的特例.

3.1 马尔可夫链

马尔可夫链(简称马氏链)是一种描述随机过程的数学模型.从t0时刻起，以一定的概率转移到状态空间S={s1,s2,…,sn}的某一个，以后各个时刻均如此.马氏链具有“无记忆”特点，即给定过程的过去和现在状态，将来状态只取决于现在的状态，又称之为马尔可夫属性.由于马氏链的转移概率仅取决于当前的状态，可为之建立一步状态转移概率矩阵T(si→sj).我们把t0时刻初始状态和以后各时刻的状态拉成一条链，称为马氏链.

我们关心一类特殊的马氏链.如果马氏链可以从每个其他状态到达每个状态，称为不可约马尔可夫链.一个具有所有正返回状态的不可约马氏链称为可遍历的马氏链[8].我们将只使用可遍历马氏链，因为它们具有唯一的稳态分布.

假设p(x)是我们的目标分布函数，我们拟获取一系列服从该分布的采样点.MCMC类算法将采样点的产生过程构成一个可遍历马尔可夫链，它呈平稳分布且收敛于稳态p(x)，所以采样点序列服从p(x).p(x)是已知的，问题的关键是如何构造可遍历马尔可夫链.

平稳状态及细节平衡条件.假如我们模拟大量物理粒子在马氏链中的随机行为[9]，粒子不断地从一个状态跃迁到另一状态，最后这些移动将达到一个动态平衡，称为平稳状态.在达到平稳状态后，从某个状态x流出去的粒子数，将会等于流回该状态的粒子数，这时系统满足 “一般平衡条件”(公式(4)的第一个等式)：

(4)

公式(4)推导过程中的最后一个等式表明此时系统处在平稳状态，即分布p(x)对马氏链是不变的，这意味着转移概率不会改变p(·)分布.为方便实现，我们试着把条件加强，由此导出公式(5)的细节平衡条件(detailed balance)：

P(x)T(x→y)=P(y)T(y→x)

(5)

可以证明，满足细节平衡条件必然满足一般平衡条件.

MCMC类算法有两种主要方法：Metropolis-Hastings算法和吉布斯抽样算法.

3.2 Metropolis算法

算法思想：引进一个提议分布q(x)，抽取候选状态(样本点)x*～q(x*丨x(t)).评估目标分布是否在候选状态x*附近具有足够大的密度，若满足，则接受x*，并将其设置为下一个状态.如果候选状态密度低于当前状态x(t)，那么很可能(但不保证)它将被拒绝.接受或拒绝候选状态的标准由以下启发式定义：

(1)如果p(x*)≥p(x(t))，接受x*作为目标分布的样本点，并作为提议分布的新状态，意味着上述操作将鼓励跳转到目标分布的高密区状态.

需要注意的是，Metropolis算法中唯一限制是提议分布必须是对称的，满足条件的这类分布有正态分布、柯西分布、t分布和均匀分布.但是如果目标分布的支持集是不对称的，例如x∈(0，+∞)，我们将考虑使用Metropolis-Hastings算法.

3.3 Metropolis-Hastings算法

为了能够使用非对称的提议分布，Metropolis-Hastings算法定义了附加校正因子c，

校正因子c调整转移算子，以确保x(t)→x*的概率等于x*→x(t)的概率，而不管提议分布是否对称.此时候选状态的x*接受率为

(6)

下面证明构造的马氏链中状态x(t)→x*的概率等于状态x*→x(t)的概率，即满足细节平衡条件.为行文方便，引入以下符号：x≜x(t)，提议分布q(x,x*)≜q(x*|x(t))，状态转移算子T(x,x*)≜T(x(t)→x*)，则接受率

首先有T(x,x*)=q(x,x*)α(x,x*)，T(x*,x)=q(x*,x)α(x*,x)，只需证明p(x)T(x,x*)=p(x*)T(x*,x)，即满足细节平衡条件.

(1)当p(x*)q(x*,x)=p(x)q(x,x*)时，α(x,x*)=α(x*,x)=1，

(7)

比较公式(7)两组等式的最后项，可知p(x)T(x,x*)=p(x*)T(x*,x)，满足细节平衡条件.

(2)当p(x*)q(x*,x)>p(x)q(x,x*)时，α(x,x*)=1，α(x*,x)<1，

(8)

比较公式(8)两组等式的最后项，可知p(x)T(x,x*)=p(x*)T(x*,x)，同样满足细节平衡条件.

(3)证明过程同(2).

注意到Metropolis算法是Metropolis-Hastings算法的特例，因为提议分布是对称的，所以q(x(t)|x*)=q(x*|x(t))，校正因子c=1.

算法描述如下：

(a)从提议分布抽样候选样本点x*～q(x*|x(t))；

(d)从均匀分布抽取u～Unif(0,1)，如果u≤α接受x*作为新状态，否则x*=x(t)，转(a).

经过充分一段时间的迭代(也称burned in，煅烧成型时间)，马氏链最终进入并收敛到平稳状态，只是不同的初始状态，煅烧成型经历的时间不同.

3.4 Gibbs采样法(Gibbs Sampler)

生成一维随机变量并不困难，对于高维分布p(x)，x=(x1,x2,…,xn)，生成各分量非独立的随机向量就非常困难.Gibbs采样法把一次一个n维样本变为“n次一维”样本，但前提是我们事先知道完全条件概率分布p(xj|x1,…,xj-1,xj+1…,xn)(数学上称为半共轭，为方便常简写为p(xj|p(x-j)).不同于MH算法，该算法不引入新的提议分布q(·)，仅通过条件分布得到以给定分布p(x)为不变分布的马氏链的转移概率.

令提议分布q(·)为

算法描述从略[10].

需要说明的是，MCMC类采样技术获取的样本点违反样本独立性假设，因为样本点间存在自相关(autocorrelation)，这是由马氏链的固有性质决定的——采样点依赖于前一个.为了减少这种轻微依赖性，第一种方法采用删减(thinning)技术，对于生成的马氏链每隔K个保留一个，其余的全丢弃；第二种块采样(Blocked Gibbs Sampler)技术，若xi,xj在给定xk的前提下并不条件独立，可将xi,xj两个变量组合在一起，并根据它们的联合分布对所有其他变量进行采样，而不单独采样每个变量；第三种折叠采样(CollapsedGibbsSampler)技术[11]，若在给定xk的前提下，xi,xj不条件独立，但在xk未知的情况下，xi,xj条件独立，此时可积分掉xk，使得xi,xj条件独立，相关内容请参照概率图模型[4].

尽管我们采取上述技术可以减少样本依赖性，但不能完全消除.自相关性的存在将降低采样精度，因为高自相关同样意味着我们采样得到的样本平均值将高于从后验分布中抽取独立样本的均值.

4 实验

假设叶斯后验分布有如下形式：

联合分布密度函数等高线见图3.

采用Gibbs 采样技术，先求条件概率分布：

(9)

这里A(y),B(x)是关于y,x函数，当指定y,x时，它们是常数.很明显，公式(9)中随机变量x,y都符合正态分布，即

图3 联合分布密度函数的等高线图

图4 Bibbs采样点轨迹图

编写并运行Gibbs采样法，初值设定在(1，6)，采样4*105次.为消除自相关，采用删减技术，K=20.观察采样点轨迹图(图4)，采样点大多集中在高密区，它们是较具代表性的样本点.

5 结语

采样技术提供给统计学家一个不同的贝叶斯推断方法，把他们从因数学计算的方便性而被限制到那些可以从理论分析中找到后验分布的模型(例如共轭分布)中解放出来，从而可以选择使用基于观测模型得到的更现实的先验分布.MCMC类采样法是了解其参数的后验分布的、对复杂模型进行推理的有效方法，这使得统计学家只需关注模型的设计而不必担心计算能力，该方法的灵活运用大大提高了应用统计学、模式识别、机器学习理论领域处理来自现实世界复杂模型的能力.