索文莉，李長國

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

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

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

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

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

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

1 ABC-SMC算法

近似貝葉斯算法的目標是獲得一個計算強度可承受，又比較好的近似后驗分布結(jié)果：

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

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

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

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

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

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

(1)

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

1)t=1

(S1)i=1；

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

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

2)t=2∶T

(S1)i=1；

2 實例分析

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

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

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

其中：f(t)是均值為0，標準差為10的激勵序列，參數(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ù)，加上隨機噪聲之后作為觀測數(shù)據(jù)yobs，據(jù)此從上述2個模型中選擇合適的模型，并得到參數(shù)估計結(jié)果.

設(shè)參數(shù)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，定義模擬數(shù)據(jù)y*與觀測數(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算法編程計算，可得到每一次迭代各模型的后驗概率分布圖，如圖1所示，每一行從左向右依次是第1～12次迭代得到的兩個模型的后驗概率結(jié)果.

圖1 可選模型的后驗概率 Figure 1 Posterior probabilities of alternative models

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

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

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

3 結(jié) 語

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