基于近似贝叶斯的模型选择

索文莉，李长国

(陆军军事交通学院 基础部，天津 300161)

在工程学、生物学、社会学和经济学等多个研究领域，学者们都会遇到模型的选择和比较问题，尤其是有多个可选模型来解释数据时，从中选出最合适的模型是具有挑战性的，常常还需要对模型背景有深度了解.

传统的模型选择中常用的IC准则是基于极大似然估计的，利用其中的边际似然来计算每个模型的可信度，但是随着模型复杂度的逐渐提高，边际似然的计算是个难题；对于很多非线性问题，其密度不一定是高斯的，而是二项、多项或者泊松等等，这种情况下，IC就不能用来比较备选模型，这限制了它的广泛应用.

与传统方法相比，贝叶斯方法已经成功解决了不同类型模型的模型选择和参数估计问题[1-4]；Green[5]利用基于贝叶斯理论的MCMC用于模型选择，但处理高维模型时，需定义额外的算法；Raftery[6]将贝叶斯因子用作模型比较，但仅仅提供了可选模型的相对比较，并没有给出后验概率的绝对取值.

近似贝叶斯计算(Approximate Bayesian Computation，简写为ABC)依托于一般的贝叶斯理论，提出了新的想法，把对似然函数的度量处理转变为观测数据与模拟数据之间的相似度，基于拒绝方法产生参数的近似后验分布样本[7]，前提是研究清楚模型的生成机理.ABC算法的思路简洁，克服了似然函数难以表示和高斯假设不成立的问题，而且不需要计算额外的标准来判别候选模型；可同时处理多个模型和大量数据，解决了MCMC的限制.因此，基因学、生物学和心理学领域都运用ABC算法处理过模型选择和参数估计问题[8-11]，其高效率的计算激励了更多学者来研究其更多可用性，同时也促进了ABC算法的快速发展和持续改进，衍生了ABC算法的其他变种算法：Beaumont[12]提出将回归方法应用于ABC的参数估计中，减弱了对阈值的要求；Toni[13]将序列蒙特卡洛(SMC)思想与ABC算法结合，Marjoram[14]提出将ABC与MCMC算法相结合.

本文首先介绍ABC-SMC算法在模型选择和参数估计中的应用思想和算法流程，然后结合实例来说明该算法的有效性，最后指出有待改进之处以及以后的工作.

1 ABC-SMC算法

近似贝叶斯算法的目标是获得一个计算强度可承受，又比较好的近似后验分布结果：

π(θ|xobs,m)∝L(xobs|θ,m)π(θ|m)

其中：m表示可选模型，π(θ|m)表示模型m条件下参数θ的先验分布，L(xobs|θ,m)是给定模型m、参数θ条件下观测数据xobs的似然函数.为克服似然函数难以精确表示或者计算冗繁的问题，ABC算法主要依赖于比较模拟数据xobs与观测数据x*的之间的距离d(xobs,x*)，如果阈值足够小，则能很好的近似真实后验分布πε(θ|xobs,m)≈π(θ|xobs,m).

应用ABC-SMC[13]算法进行参数估计和模型选择，核心思想是：通过带权重的样本表示后验密度，通过一系列中间过程的后验分布样本，根据研究者定义的收敛标准，逐渐收敛得到最优近似后验分布.

算法的第1次迭代过程：从一个任意阈值ε开始，避免条件太苛刻导致接受率太低，降低计算效率；从先验分布π(m)与π(θ|m)中抽取随机样本，其中π(m)表示可选模型的先验分布，在模型m、参数θ的条件下，生成模拟数据x*；计算距离d(xobs,x*)，与阈值ε相比，决定接受还是拒绝；重复上述过程直到产生N个接受的参数粒子；对接受的参数粒子分别赋予权重1/N.

算法的第t(t=2∶T)次迭代过程：设置阈值序列εt，满足ε1>ε2…>εT，最终选择的阈值εT，依赖于专业人员的目的与要求.阈值序列εt的选择可手动也可动态调整，例如第二次迭代的阈值可设为前一次迭代距离的p%，这也是最常用的阈值序列选择法则，直观简单易定义，也可保持较合适的接受率.当然也可自己手动定义阈值序列，逐渐减小到研究者想要达到的阈值，本文选择这种方式.

然后计算d(xobs,x*)，如果d(xobs,x*)≤εt，接受新粒子，否则拒绝；重复这个过程直到生成N个粒子，更新迭代次数t；依据更新粒子的权重和阈值序列迭代，直至得到最后一次迭代得到的粒子总体，可得模型mk的近似边缘后验分布为

(1)

具体的算法流程如下所示.

1)t=1

(S1)i=1；

(S5)重复步骤S2～S4，直到接受N个参数粒子；

(S6)根据式(1)计算模型mk的后验分布，作为下一次迭代模型的先验分布；

2)t=2∶T

(S1)i=1；

2 实例分析

假设有两个可选模型，

ay″+by′+cy=f(t) (m1)

ay″+by′+cy+dy2=f(t) (m2)

其中：f(t)是均值为0，标准差为10的激励序列，参数为a,b,c,d,e.

给定初始值y(0)=0,y′(0)=0，设模型m2中参数a=1,b=0.05,c=50,d=100，采用Runge-Kutta方法生成数据，加上随机噪声之后作为观测数据yobs，据此从上述2个模型中选择合适的模型，并得到参数估计结果.

设参数a,b,c,d,e先验分布为U(0.1,10),U(0.01,1),U(5,500),U(10,500),U(10,1 000)，每次迭代生成N个粒子，考虑到计算时间此例取N=50，定义模拟数据y*与观测数据yobs的距离函数为

阈值序列为ε=(50 20 10 5 2 0.5 0.3 0.1 0.05 0.01 0.001 0.0001)，根据ABC-SMC算法编程计算，可得到每一次迭代各模型的后验概率分布图，如图1所示，每一行从左向右依次是第1～12次迭代得到的两个模型的后验概率结果.

图1 可选模型的后验概率 Figure 1 Posterior probabilities of alternative models

随着迭代次数的增加，可以看出模型2为最终选择的模型，这与数据的生成模型一致，说明模型选择正确.同时可以得到每次迭代后各模型对应参数的后验直方图，这里只输出模型2最后一次迭代得到的参数可得各参数的后验直方图，如图2所示，取后验样本的均值(中位数)作为参数估计值，与真实值比较可见算法的估计精度，如表1所示.

图2 模型参数的后验直方图Figure 2 posterior histogram of model parameters

表1 模型的参数估计结果Table 1 parameter estimation result of model

3 结 语

现实问题越来越趋于复杂，模型对应的似然函数也随之越来越冗繁，难以精确写出或者进行数值计算，贝叶斯统计方法逐渐称为研究者们常用的算法之一.ABC-SMC算法作为近似贝叶斯算法的改进算法，极大提高了计算效率，估计精度较高，同时可用于模型的选择问题，得到很好的结果.但是在算法过程中，粒子受核函数的影响，实施中需要选择一定数量的超参数，定义距离函数、阈值序列，这些都需要谨慎选择，因为算法的实施表现非常依赖它们，选取的不好会导致计算时间难以承受，效率较低.