第三节　贝叶斯统计的常见算法

一、MCMC算法

贝叶斯统计基本理论和方法虽然简单易懂，但是在实际应用中由于后验分布涉及多重积分，无法给出解析解及数值解，因此在20世纪90年代前贝叶斯统计多数模型停留在理论研究阶段，只有少数模型可以采用共轭先验分布及近似解获得结果。20世纪90年代后随着随机模拟技术及计算机技术的发展，贝叶斯统计领域被广泛应用，使得贝叶斯统计计算得到迅猛发展，尤其是马尔可夫链蒙特卡罗模拟方法的引进。

马尔可夫链蒙特卡罗模拟（Markov Chain Monte Carlo Simulation，MCMC），由两部分内容组成，即采用马尔可夫链（Markov Chain，MC）获得目标分布的抽样样本，然后通过蒙特卡罗（Monte Carlo，MC）积分获得目标分布的相应结果。

蒙特卡罗方法由美国在第二次世界大战中研制原子弹的“曼哈顿计划”的成员塔尼斯拉夫·乌拉姆和约翰·冯·诺依曼于20世纪40年代首先提出。数学家约翰·冯·诺伊曼用驰名世界的赌城——摩纳哥的Monte Carlo的名字来命名这种方法。蒙特卡罗方法是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。蒙特卡罗方法的解题过程可以归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。具体如下：

若随机变量X的密度函数为f（x），假设g为一个可积函数，则随机变量Y=g（X）的期望为：。

如果从X的分布中获得一些随机数，则E［g（X）］的无偏估计即为相应的样本均值。通俗地理解即当一个积分求解困难时，我们可以通过构建一个函数的期望，从随机变量所服从分布中抽样，计算每个样本值的函数值，求得均值即为积分值。如为估计，若X₁…X_N为均匀分布（0，1）的抽样样本，则由强大数定律可知，以概率1收敛到期望E［g（X）］。

［例1.7］采用蒙特卡罗方法计算。提示该积分解析解为。

［解］（0，1）的均匀分布的密度函数为f（x）=1，从（0，1）的均匀分布中抽取n=10 000样本，则有，可由R语言编写简单代码实现，如：

计算结果为theta.hat=0.631 417，即，而积分值，两者计算结果相近。

蒙特卡罗方法将积分问题转化为从目标概率密度f（x）中产生随机样本，在实际计算中，后验分布较为复杂，无法直接抽样，必须解决如何方便得到各种复杂概率分布的对应采样样本集的问题，否则仍无法计算蒙特卡罗积分。马尔可夫链为解决抽样问题提供了很好的办法。在MCMC方法中，首先建立一个马尔可夫链，运行足够长时间后，使其收敛于平稳分布，产生的随机样本就是从达到平稳状态的马尔可夫链中产生的样本。

假设序列状态是…，X_t-2，X_t-1，X_t，X_t+1，…，那么在时刻X_t+1的状态的条件概率仅仅依赖于时刻X_t，即：。

具有此性质的序列状态称为马尔可夫链性质。不可约、非周期、正常返的马尔可夫链在运行充分长时间后将收敛于平稳分布，与初始分布无关。马尔可夫链的实值函数的遍历均值依概率1收敛于极限分布下的均值。遍历均值做合适的变化将依分布收敛于标准正态分布。这三条性质从理论上保证了MCMC算法的有效性。马尔可夫链的理论基础较为复杂，读者有兴趣可参考相关著作。MCMC算法通过构建马尔可夫链，涉及一大类抽样方法，包括Metropolis抽样、Metropolis-Hastings抽样、Gibbs抽样、随机游走抽样等。最为常见的是Metrapolis-Hastings抽样和Gibbs抽样。

Metropolis-Hastings抽样简称M-H抽样，这个算法首先由Metropolis提出，后被Hastings改进，因此被称之为Metropolis-Hastings抽样，主要是解决Metrapolis抽样效率低的问题。具体抽样算法如下：

1.构造合适的建议分布g（·X_t ），该分布满足不可约、正常返、非周期的正则化条件；

2.从分布g（·X_t ）中产生X₀；

3.重复下列过程，直至马尔可夫链达到平稳状态：

（1）从g（·X_t ）产生Y；

（2）从均匀分布U（0，1）中产生U；

（3）若，则接受Y并令X_t+1=Y，否则X_t+1=X；接受概率为α（X_t，Y）=

（4）增加t，返回到（1）。

M-H抽样虽然很好地解决了蒙特卡罗方法需要的任意概率分布的样本集问题。但是M-H抽样有两个缺点：一是需要计算接受率，在高维时计算量大，计算接受率导致算法收敛时间变长。二是有些高维数据，特征的条件概率分布好求，但是特征的联合分布不好求。Gibbs抽样很好地解决了M-H抽样的缺点。Gibbs抽样是M-H抽样的一个特例，重新构建了细致平衡条件，使得抽样的接受率等于1，即保留了所有的抽样样本。Gibbs抽样是基于条件概率分布进行的，其基本做法是从联合概率分布定义条件概率分布，依次对条件概率分布进行抽样得到样本的序列。可以证明这样的抽样过程是在一个马尔可夫链上的随机游走，每一个样本对应着马尔可夫链的状态，平稳分布就是目标的联合分布。预烧期之后的样本就是联合分布的随机样本。具体抽样算法如下：

1.随机初始化｛x_i：i=1，…，n｝；

2.对t=0，1，2…循环抽样：

（1）

（2）

（3）…

（4）

（5）…

（6）

Gibbs抽样大大提高了抽样的效率，成为众多软件的算法基础，如WinBUGS、JAGS等。下面简单演示正态分布参数的Gibbs抽样过程。

［例1.8］假设X=（X₁，X₂，…X_n）为从N（θ，σ²）总体中抽取的随机样本。取θ的先验分布为N（μ，τ²），σ²先验分布为IG（a，b），其中μ，σ²，a，b已知，采用Gibbs抽样对θ，σ²进行抽样。

［解］（θ，σ²）的后验分布核为

化简得到参数的条件分布核如下：

继而化简可知：

由上述结果可知，参数θ的条件分布为含有σ²的正态分布，而σ²的条件分布为含有θ的IG分布；在给定初始值后，Gibbs抽样过程示意图如图1-2：

［例1.9］以Mackowiak等（1992）的数据为例，该数据记录了100个人的体温（华氏）、性别和每分钟的心跳数。实验的目的是检验Carl Wunderlich的观点——健康成年人的平均体温为37℃（98.6℉），试采用Gibbs抽样进行贝叶斯估计。

［解］采用N（0，100）及IG（0.01，0.01）弱信息先验分布。R语言代码如下：

体温均数、标准差抽样分布及迭代100路径分布图如图1-3：

二、MCMC抽样的收敛性诊断

MCMC方法通过构造一个概率转移矩阵，建立一个以分布π（x）为平稳分布的马尔可夫链来得到π（x）的样本，通过随机抽样得到的这些样本进行各种统计推断。马尔可夫链的极限分布是目标分布，但是如何判断马尔可夫链收敛于目标分布？这就涉及MCMC的收敛，目前尚无准确无误判断收敛的方法，但是有可以辅助判断是否收敛的手段。所谓算法的收敛性指的是所得到的马尔可夫链是否达到了平稳状态。如果达到平稳状态，则认为所获得的样本来自目标分布。一般情况下，并不知道算法需要运行多长时间才能到达平稳状态，因此对马尔可夫链收敛性进行监视。

MCMC产生的马尔可夫链的收敛性受到三个因素的影响：①初始值：用来初始化马尔可夫链。需要利用一切可得到的信息尽量使得初始值离参数真值更近。如果初始值远离后验密度的最高区域，且算法的迭代次数不足以消除初始值的影响，则会对后验推断产生影响。②后验分布密度：如果后验密度的形状是单峰的，问题相对简单。如果后验密度的形状是多峰的，则要小心伪收敛现象，也就是可能会陷入非最大值的其他极值点。③工具分布：好的工具分布的中心区域应当与后验密度的中心区域有较多的重合，尾部比较厚，支撑包含了后验密度的支撑。

收敛涉及相关的几个术语：①预烧期（burn-in period）：在MCMC机制中，用以保证马尔可夫链达到平稳状态所运行的时间称为预烧期。一般只要链运行时间足够长，预烧期对后验推断几乎没有影响。②筛选间隔或抽样步长（thinning interval or sampling lag）：马尔可夫链产生的样本并不是相互对立的，因此采用每间隔L个样本来抽取一个（近似）独立的样本。③蒙特卡罗误差（MC error）：其度量了估计因随机波动而导致的波动性。在计算感兴趣的参数时，其精度随着样本容量递增，与样本容量成反比。有组平均方法及窗口估计量计算方法。

收敛性的判断主要包括两大类方法，一是图示法，包括踪迹图、遍历均值图、自相关图、核密度估计图；二是收敛性指标，包括Gelman-Rubin法、Geweke法、Raftery-Lewis法、Heidelberg-Welch法等。

踪迹图本质是一个时间序列图，收敛的马尔可夫链应当是没有趋势与周期的，它会在一个水平线上下小幅度波动（形状上类似白噪声序列）。当马尔可夫链收敛时遍历均值将接近一条直线，只有小幅波动。马尔可夫链收敛后，抽样样本量越大，核密度估计图越光滑。下面继续以体温数据为例，绘制上述图形，如图1-4和图1-5，并判断马尔可夫链的收敛性。R语言代码如下：

由踪迹图可知，theta及sigmasqr围绕均值两端均匀分布；核密度估计图提示二者曲线光滑，theta呈现对称分布，sigmasqr呈左偏态分布；自相关图提示相关性低，最大自相关系数为0.01左右；遍历均值图提示，在迭代早期，遍历均值很快达到稳定状态。以上结果提示抽样分布收敛性好。

Gelman-Rubin法由Gelman及Rubin提出来，又称为方差比法。基本思想是当多条马尔可夫链平稳收敛后，马尔可夫链的统计特征应该是一样的，利用类似方差分析的方法构造出一个统计量诊断马尔可夫链的收敛性，分别计算链间的方差和链内的方差，然后计算二者的加权组合来估计方差，最后由加权方差估计量及链内方差估计量计算获得潜在尺度缩减因子（potential scale reduction factor，PSRF），简写为R。一般而言，R<1.1提示马尔可夫链收敛。具体计算过程可参考相关著作。R语言coda包可以给出Gelman-Rubin方法的结果。首先生成两组2 000个标准正态分布样本作为两条链；生成一组2 000个N（2，1）正态分布样本为一条链，进行Gelman-Rubin检验，如图1-6、图1-7。具体如下：

X1由两组2 000个标准正态分布样本组成，可视为同一目标分布产生的两条链；X2由两组不同正态分布样本，可视为没有收敛到目标分布。由图1-7可知X1的两条马尔可夫链混杂在一起，均匀分布在0两侧，同时核密度估计图示两条链的样本呈现一条光滑的曲线。而X2的两条马尔可夫链分别来自N（0，1）和N（3，1），因此踪迹图并未混在一起；核密度估计图也呈现出双峰表现。由Gelman-Rubin图可知，X1在迭代500次后，R≈1；而X2的R≈4；说明X2两条链并未收敛。

Geweke检验法，包括了数量诊断与图形诊断两个部分。这种方法的基本思想是：如果马尔可夫链是收敛的，那么这个马尔可夫链前一部分与后一部分的均值应该是相等的。在马尔可夫链前一部分与后一部分渐进独立的假设条件下，Geweke构造的检验统计量渐进服从标准正态分布。数量诊断中，Geweke检验对每个分量马尔可夫链算出检验统计量的值，当值小于1.96或2的时候，马尔可夫链被认为是收敛的。图形诊断中，采用同样的方法算出一个统计量的值，然后，把最前边的一小段马尔可夫链切除，再用同样的方法算出第二个统计量的值，以此类推。一共算出了若干个统计量值之后，最后在坐标系上画出他们的位置，同时用虚线画出某一置信度的置信区间（一般认为是0.95）形成的置信带。如果绝大多数统计量值在可信带内，那么认为马尔可夫链是收敛的。R语言coda包的geweke.diag命令可实现，具体略。

三、贝叶斯统计的其他算法

贝叶斯后验分布因涉及多个参数，计算时面临多重积分的挑战，在早期限制了其在多元统计建模中的使用。MCMC方法的出现很好地解决了计算困难的问题，但仍存在一些问题，比如M-H算法接受率低下、高维参数空间Gibbs抽样效率低下及所有基于马尔可夫链抽样样本收敛性的问题。随着大型数据集和复杂模型的出现，出现了更先进的算法，如期望最大化（EM）法、积分嵌套拉普拉斯逼近（integrated nested Laplace approximation，INLA）、哈密尔顿蒙特卡罗（Hamiltonian Monte Carlo，HMC）、变分逼近（variational approximation，VA），这些算法有的可以避免收敛性问题，有的可以提高计算效率。这些先进的算法往往涉及复杂的数学及统计原理，本小节只简单介绍INLA及HMC的基本思想，读者感兴趣可参考相关著作。

INLA是综合贝叶斯推断、潜变量高斯模型、高斯马尔可夫随机场及拉普拉斯近似的方法。首先按条件概率重新编写后验分布，采用拉普拉斯近似值替代边际分布p（xθ，y ），再计算分母边际分布的拉普拉斯近似值，采用高斯马尔可夫随机场的方法。其中拉普拉斯近似是基于泰勒级数展开式，将概率密度函数对数化后，按泰勒级数展开，分别计算二阶导数后，可以与正态分布关联，将求解的概率密度函数逼近特定均数及方差的正态分布。Finn Lindgren等人在2011年进行的研究提出了一种针对INLA的随机偏微分方程（SPDE）方法，该方法可以在较大的地理区域内进行更实际的空间建模，这些地理区域的曲面具有点数据，从而具有连续的协方差函数。因此，该方法在处理空间数据上具有极大的优势。R语言的inla包可以实现贝叶斯统计的非迭代抽样计算。

HMC方法又称为混合蒙特卡罗（Hybrid Monte Carlo），其采用系统动力学而不是概率分布来建立马尔可夫链的未来状态。哈密尔顿动力学是经典力学的重要理论，采用位置、动量及时间来刻画物体的状态。按能量守恒定律，物体总能量等于势能和动能之和，当物体在整个系统中运动时，势能和动能之间相互转化，可通过Hamiltonian方程的微分方程来描述。描述运动系统的当前状态需要有物体所处的位置式（position）和物体的动量。位置对应了统计中参数的值，而动量则是人为添加的一个变量，也被称为辅助变量。通过统计力学的正则分布与目标分布关联起来，通过欧拉法或蛙跳法进行模拟。HMC法的关键在于选取轨迹的长度及步长，当选择不当时会增加拒绝率，对于高维复杂模型可能耗时。No-U-Turn抽样（NUTS）是HMC的扩展，可解决上述问题。Stan是贝叶斯推断的C++库，基于No-U-Turn采样器（NUTS）进行抽样，进而对后验分布进行参数估计。不同统计软件有不同的Stan包，如rstan、pystan等。

四、贝叶斯统计常用软件简介

尽管贝叶斯统计涉及复杂的计算，但是有不同软件及程序包可以方便地进行计算。这些软件包括专用软件及通用软件。专用软件主要包括WinBUGS、OpenBUGS、JAGS及Stan等。通用软件包括R、SAS、Matlab、Python等。通用软件中可以自行编程进行计算，但多数是依靠相应的包，如R语言中的R2WinBUGS、rjags、mcmcglmm、mcmcpack等；SAS软件中的proc mcmc包；Stata 16.0软件中的bayes系列命令；Python软件中的pymcmc、pymc3等。读者只须根据自己的使用习惯进行选择。推荐读者采用R语言结合专用软件进行贝叶斯Meta分析。

BUGS是Bayesian Inference Using Gibbs Sampling的缩写，指基于Gibbs抽样的贝叶斯统计推断。BUGS软件最初于1989年由位于英国剑桥的生物统计学研究所（Biostatistics the Medical Research Council，Cambridge，United Kingdom）研制，现在由这个研究所和位于伦敦的S.t Mary’s皇家学院医学分院（the Imperial College School of Medicine）共同开发。BUGS所使用的编程语言类似R语言，通过设定贝叶斯模型的似然函数及先验分布，并将其自动转换成马尔可夫链的转移核，迭代后获得后验分布（目标分布）的抽样数据。BUGS主要包括WinBUGS和OpenBUGS，WinBUGS目前已暂停更新。OpenBUGS作为WinBUGS的开源版本，在持续更新中，目前的版本为OpenBUGS323，可于Windows、Linux等系统。尽管Open-BUGS在持续更新，但其基本代码、数据、初始值及操作与WinBUGS无明显差异。主要差异可参考：http：//www.openbugs.net/w/OpenVsWin。WinBUGS下载地址为https：//www.mrcbsu.cam.ac.uk/wp-content/uploads/2018/11/winbugs143_unrestricted.zip；OpenBUGS下载地址为http：//www.openbugs.net/w/OpenBUGS_3_2_3？action=AttachFile&do=get&target=Open-BUGS323setup.exe。

JAGS为Just Another Gibbs sampler的缩写，是由Martyn Plummer创作的基于Gibbs抽样的软件，可在不同的操作系统内使用。目前主要是与R语言相关的包联合应用，如R2jags、runjags、rjags等。JAGS语法命令与BUGS相似，多数代码是通用的；最新的版本是JAGS4.3.0，下载地址为https：//sourceforge.net/projects/mcmc-jags/。

Stan是用于贝叶斯统计建模和高性能统计计算的最新平台。可在社会学、生物和物理科学、工程和商业领域进行统计建模、数据分析和预测。用户使用Stan的概率编程语言指定对数密度函数，使用No-U-Turn抽样器进行贝叶斯统计推断、带变分推理（ADVI）的近似贝叶斯推理、带有优化的惩罚最大似然估计（L-BFGS）。更加方便地进行复杂的贝叶斯模型的计算。