基于近似貝葉斯的模型選擇

索文莉，李長(zhǎng)國(guó)

(陸軍軍事交通學(xué)院 基礎(chǔ)部，天津 300161)

在工程學(xué)、生物學(xué)、社會(huì)學(xué)和經(jīng)濟(jì)學(xué)等多個(gè)研究領(lǐng)域，學(xué)者們都會(huì)遇到模型的選擇和比較問(wèn)題，尤其是有多個(gè)可選模型來(lái)解釋數(shù)據(jù)時(shí)，從中選出最合適的模型是具有挑戰(zhàn)性的，常常還需要對(duì)模型背景有深度了解.

傳統(tǒng)的模型選擇中常用的IC準(zhǔn)則是基于極大似然估計(jì)的，利用其中的邊際似然來(lái)計(jì)算每個(gè)模型的可信度，但是隨著模型復(fù)雜度的逐漸提高，邊際似然的計(jì)算是個(gè)難題；對(duì)于很多非線性問(wèn)題，其密度不一定是高斯的，而是二項(xiàng)、多項(xiàng)或者泊松等等，這種情況下，IC就不能用來(lái)比較備選模型，這限制了它的廣泛應(yīng)用.

與傳統(tǒng)方法相比，貝葉斯方法已經(jīng)成功解決了不同類型模型的模型選擇和參數(shù)估計(jì)問(wèn)題[1-4]；Green[5]利用基于貝葉斯理論的MCMC用于模型選擇，但處理高維模型時(shí)，需定義額外的算法；Raftery[6]將貝葉斯因子用作模型比較，但僅僅提供了可選模型的相對(duì)比較，并沒(méi)有給出后驗(yàn)概率的絕對(duì)取值.

近似貝葉斯計(jì)算(Approximate Bayesian Computation，簡(jiǎn)寫為ABC)依托于一般的貝葉斯理論，提出了新的想法，把對(duì)似然函數(shù)的度量處理轉(zhuǎn)變?yōu)橛^測(cè)數(shù)據(jù)與模擬數(shù)據(jù)之間的相似度，基于拒絕方法產(chǎn)生參數(shù)的近似后驗(yàn)分布樣本[7]，前提是研究清楚模型的生成機(jī)理.ABC算法的思路簡(jiǎn)潔，克服了似然函數(shù)難以表示和高斯假設(shè)不成立的問(wèn)題，而且不需要計(jì)算額外的標(biāo)準(zhǔn)來(lái)判別候選模型；可同時(shí)處理多個(gè)模型和大量數(shù)據(jù)，解決了MCMC的限制.因此，基因?qū)W、生物學(xué)和心理學(xué)領(lǐng)域都運(yùn)用ABC算法處理過(guò)模型選擇和參數(shù)估計(jì)問(wèn)題[8-11]，其高效率的計(jì)算激勵(lì)了更多學(xué)者來(lái)研究其更多可用性，同時(shí)也促進(jìn)了ABC算法的快速發(fā)展和持續(xù)改進(jìn)，衍生了ABC算法的其他變種算法：Beaumont[12]提出將回歸方法應(yīng)用于ABC的參數(shù)估計(jì)中，減弱了對(duì)閾值的要求；Toni[13]將序列蒙特卡洛(SMC)思想與ABC算法結(jié)合，Marjoram[14]提出將ABC與MCMC算法相結(jié)合.

本文首先介紹ABC-SMC算法在模型選擇和參數(shù)估計(jì)中的應(yīng)用思想和算法流程，然后結(jié)合實(shí)例來(lái)說(shuō)明該算法的有效性，最后指出有待改進(jìn)之處以及以后的工作.

1 ABC-SMC算法

近似貝葉斯算法的目標(biāo)是獲得一個(gè)計(jì)算強(qiáng)度可承受，又比較好的近似后驗(yàn)分布結(jié)果：

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

其中：m表示可選模型，π(θ|m)表示模型m條件下參數(shù)θ的先驗(yàn)分布，L(xobs|θ,m)是給定模型m、參數(shù)θ條件下觀測(cè)數(shù)據(jù)xobs的似然函數(shù).為克服似然函數(shù)難以精確表示或者計(jì)算冗繁的問(wèn)題，ABC算法主要依賴于比較模擬數(shù)據(jù)xobs與觀測(cè)數(shù)據(jù)x*的之間的距離d(xobs,x*)，如果閾值足夠小，則能很好的近似真實(shí)后驗(yàn)分布πε(θ|xobs,m)≈π(θ|xobs,m).

應(yīng)用ABC-SMC[13]算法進(jìn)行參數(shù)估計(jì)和模型選擇，核心思想是：通過(guò)帶權(quán)重的樣本表示后驗(yàn)密度，通過(guò)一系列中間過(guò)程的后驗(yàn)分布樣本，根據(jù)研究者定義的收斂標(biāo)準(zhǔn)，逐漸收斂得到最優(yōu)近似后驗(yàn)分布.

算法的第1次迭代過(guò)程：從一個(gè)任意閾值ε開(kāi)始，避免條件太苛刻導(dǎo)致接受率太低，降低計(jì)算效率；從先驗(yàn)分布π(m)與π(θ|m)中抽取隨機(jī)樣本，其中π(m)表示可選模型的先驗(yàn)分布，在模型m、參數(shù)θ的條件下，生成模擬數(shù)據(jù)x*；計(jì)算距離d(xobs,x*)，與閾值ε相比，決定接受還是拒絕；重復(fù)上述過(guò)程直到產(chǎn)生N個(gè)接受的參數(shù)粒子；對(duì)接受的參數(shù)粒子分別賦予權(quán)重1/N.

算法的第t(t=2∶T)次迭代過(guò)程：設(shè)置閾值序列εt，滿足ε1>ε2…>εT，最終選擇的閾值εT，依賴于專業(yè)人員的目的與要求.閾值序列εt的選擇可手動(dòng)也可動(dòng)態(tài)調(diào)整，例如第二次迭代的閾值可設(shè)為前一次迭代距離的p%，這也是最常用的閾值序列選擇法則，直觀簡(jiǎn)單易定義，也可保持較合適的接受率.當(dāng)然也可自己手動(dòng)定義閾值序列，逐漸減小到研究者想要達(dá)到的閾值，本文選擇這種方式.

然后計(jì)算d(xobs,x*)，如果d(xobs,x*)≤εt，接受新粒子，否則拒絕；重復(fù)這個(gè)過(guò)程直到生成N個(gè)粒子，更新迭代次數(shù)t；依據(jù)更新粒子的權(quán)重和閾值序列迭代，直至得到最后一次迭代得到的粒子總體，可得模型mk的近似邊緣后驗(yàn)分布為

(1)

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

1)t=1

(S1)i=1；

(S5)重復(fù)步驟S2～S4，直到接受N個(gè)參數(shù)粒子；

(S6)根據(jù)式(1)計(jì)算模型mk的后驗(yàn)分布，作為下一次迭代模型的先驗(yàn)分布；

2)t=2∶T

(S1)i=1；

2 實(shí)例分析

假設(shè)有兩個(gè)可選模型，

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

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

其中：f(t)是均值為0，標(biāo)準(zhǔn)差為10的激勵(lì)序列，參數(shù)為a,b,c,d,e.

給定初始值y(0)=0,y′(0)=0，設(shè)模型m2中參數(shù)a=1,b=0.05,c=50,d=100，采用Runge-Kutta方法生成數(shù)據(jù)，加上隨機(jī)噪聲之后作為觀測(cè)數(shù)據(jù)yobs，據(jù)此從上述2個(gè)模型中選擇合適的模型，并得到參數(shù)估計(jì)結(jié)果.

設(shè)參數(shù)a,b,c,d,e先驗(yàn)分布為U(0.1,10),U(0.01,1),U(5,500),U(10,500),U(10,1 000)，每次迭代生成N個(gè)粒子，考慮到計(jì)算時(shí)間此例取N=50，定義模擬數(shù)據(jù)y*與觀測(cè)數(shù)據(jù)yobs的距離函數(shù)為

閾值序列為ε=(50 20 10 5 2 0.5 0.3 0.1 0.05 0.01 0.001 0.0001)，根據(jù)ABC-SMC算法編程計(jì)算，可得到每一次迭代各模型的后驗(yàn)概率分布圖，如圖1所示，每一行從左向右依次是第1～12次迭代得到的兩個(gè)模型的后驗(yàn)概率結(jié)果.

圖1 可選模型的后驗(yàn)概率 Figure 1 Posterior probabilities of alternative models

隨著迭代次數(shù)的增加，可以看出模型2為最終選擇的模型，這與數(shù)據(jù)的生成模型一致，說(shuō)明模型選擇正確.同時(shí)可以得到每次迭代后各模型對(duì)應(yīng)參數(shù)的后驗(yàn)直方圖，這里只輸出模型2最后一次迭代得到的參數(shù)可得各參數(shù)的后驗(yàn)直方圖，如圖2所示，取后驗(yàn)樣本的均值(中位數(shù))作為參數(shù)估計(jì)值，與真實(shí)值比較可見(jiàn)算法的估計(jì)精度，如表1所示.

圖2 模型參數(shù)的后驗(yàn)直方圖Figure 2 posterior histogram of model parameters

表1 模型的參數(shù)估計(jì)結(jié)果Table 1 parameter estimation result of model

3 結(jié) 語(yǔ)

現(xiàn)實(shí)問(wèn)題越來(lái)越趨于復(fù)雜，模型對(duì)應(yīng)的似然函數(shù)也隨之越來(lái)越冗繁，難以精確寫出或者進(jìn)行數(shù)值計(jì)算，貝葉斯統(tǒng)計(jì)方法逐漸稱為研究者們常用的算法之一.ABC-SMC算法作為近似貝葉斯算法的改進(jìn)算法，極大提高了計(jì)算效率，估計(jì)精度較高，同時(shí)可用于模型的選擇問(wèn)題，得到很好的結(jié)果.但是在算法過(guò)程中，粒子受核函數(shù)的影響，實(shí)施中需要選擇一定數(shù)量的超參數(shù)，定義距離函數(shù)、閾值序列，這些都需要謹(jǐn)慎選擇，因?yàn)樗惴ǖ膶?shí)施表現(xiàn)非常依賴它們，選取的不好會(huì)導(dǎo)致計(jì)算時(shí)間難以承受，效率較低.