基于MCMC方法的貝葉斯統(tǒng)計模型應(yīng)用研究

摘"要：為驗證貝葉斯統(tǒng)計后驗分布計算的收斂性，文章研究并構(gòu)建了線性回歸和時間序列兩類經(jīng)典的貝葉斯統(tǒng)計模型。基于馬爾可夫鏈蒙特卡羅方法，運用birdextinct數(shù)據(jù)集驗證了貝葉斯線性回歸模型，并分析鳥類的平均滅絕時間與平均筑巢數(shù)、種群規(guī)模、棲息狀態(tài)三個量之間的關(guān)系。此外，根據(jù)19822021年居民消費指數(shù)（CPI）的序列，結(jié)合差分處理數(shù)據(jù)，建立了AR（p）時間序列預(yù)測模型。結(jié)果表明，貝葉斯線性回歸模型和貝葉斯時間序列預(yù)測模型的馬爾可夫鏈均收斂，且貝葉斯時間序列預(yù)測模型的預(yù)測誤差僅為0.78%。MCMC方法能有效應(yīng)用于貝葉斯統(tǒng)計模型分析。

關(guān)鍵詞：MCMC 方法；貝葉斯統(tǒng)計；收斂性診斷；時間序列預(yù)測

中圖分類號：O212.8"""""文獻標(biāo)識碼：A"""""""文章編號：20959699（2024）03007005

貝葉斯統(tǒng)計方法憑借其在先驗分布上的數(shù)據(jù)統(tǒng)計優(yōu)勢，廣泛應(yīng)用于經(jīng)濟學(xué)、社會科學(xué)、自然科學(xué)、醫(yī)學(xué)和心理學(xué)等領(lǐng)域[1]。貝葉斯統(tǒng)計的關(guān)鍵計算包括兩個方面，其一是先驗分布的確定；其二是后驗分布的計算[2]。目前，先驗分布的確定方法主要包括：無信息先驗分布、共軛先驗分布、概率匹配先驗分布、參照先驗分布、專家經(jīng)驗確定法、最大熵方法以及隨機加權(quán)法等[34]。而在貝葉斯統(tǒng)計方法上，以往一般采用步驟繁多、難度大的分析方法獲取信息，已經(jīng)無法滿足現(xiàn)代信息計算技術(shù)和貝葉斯統(tǒng)計發(fā)展的需要[5]。為了解決這一問題，Metropolis等[6]學(xué)者另辟蹊徑，充分結(jié)合現(xiàn)代計算技術(shù)，基于馬爾可夫理論和蒙特卡羅模擬方法，使用馬爾可夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）方法，回避后驗分布中的計算難題，實現(xiàn)了對后驗分布的直接模擬和分析，以獲得相關(guān)統(tǒng)計特征值信息。王緒民等[7]采用BayesianMCMC算法，忽視先驗分布與似然函數(shù)的共軛性，從而減少計算冗余量，提高了水利工程報價數(shù)據(jù)分布的準(zhǔn)確性。王若嬋等[8]通過MCMC算法與顫振裕度法結(jié)合，進行樣本采樣與概率分布獲取，實現(xiàn)了紊流激勵下顫振邊界結(jié)構(gòu)數(shù)據(jù)的預(yù)測。胡紅波[9]研究了概率因素對不確定度概率的影響問題，利用MCMC算法提高抽樣樣本的相關(guān)性，從而實現(xiàn)復(fù)雜場合下不確定度樣本的預(yù)測。然而，已有的研究多集中于改進方法與聯(lián)合應(yīng)用，尚未對貝葉斯統(tǒng)計方法在回歸分析與預(yù)測方面的應(yīng)用進行深入研究[10]。

文中分別研究了貝葉斯線性回歸模型和貝葉斯時間序列預(yù)測模型。通過模型參數(shù)的收斂性與預(yù)測能力，為基于MCMC方法的貝葉斯統(tǒng)計模型的構(gòu)建與應(yīng)用研究提供了理論支撐。

1"相關(guān)理論

1.1"MCMC算法

MCMC算法由兩部分組成：蒙特卡羅模擬方法和馬爾可夫鏈。MCMC算法是通過計算假設(shè)概率來探討貝葉斯統(tǒng)計中的積分計算問題[9]。在給定概率密度函數(shù)f（x|θ）的情況下，需要求解g（θ）的期望值，可以用后驗分布表示為：

g∧θ=∫Θgθπθ|xdθ（1）

其中，g∧θ即為后驗分布的均值。顯然，此均值可由式（2）近似求得：

g-=1m∑mi=1g（θi）（2）

其中，θ1，θ2，…，θm取自后驗分布πθ|x的容量為m的樣本。

1.2"貝葉斯AR（p）模型

貝葉斯AR（p）模型稱為自回歸模型，描述的是當(dāng)前值與歷史值之間的關(guān)系，并通過使用變量自身的歷史數(shù)據(jù)來預(yù)測模型。通常，一個隨機時間序列用一組按時間順序排列的隨機變量y1，y2，…，yt來表示，記為{yt}，設(shè)θ1，θ2，…，θp為模型的自回歸系數(shù)，p為模型的階數(shù)，則其AR（p）模型[10]為：

yt=θ1yt-1+θ2yt-2+…+θpyt-p+εt（3）

其中，t=1，2，…，n。即在t時刻的隨機變量yt的取值是前p期yt-1，yt-2，…，yt-p的多元線性函數(shù)。誤差項是當(dāng)期的隨機干擾εt，誤差項服從正態(tài)分布N（0，τ-1）且相互獨立，其中τ為εt的精度。在文中，假定AR（p）模型是一個平衡過程，并采用基于MCMC算法的貝葉斯隨機搜索模型選擇方法，建立特定數(shù)據(jù)序列的自回歸AR（ p）模型[10]。設(shè)xt為原始數(shù)據(jù)，yt是對xt進行差分后的數(shù)據(jù)序列，將新得到的數(shù)據(jù)序列作為初始值y0，y1，…，yt-p，設(shè)θ=θ1，θ2，…，θpT，此時隨機變量yt服從正態(tài)分布，即：

yt|yt-1，yt-2，…，yt-p=Nθ1yt-1+θ2yt-2+…+θpyt-p，τ-1（4）

則模型的似然函數(shù)為：

Lθ，τ=fy1，y2，…，yn|θ，τ∝τ2πn2exp-τ2∑nt=1yt-θ1yt-1+θ2yt-2+…+θpyt-p。

模型的未知參數(shù)的先驗分布采用正態(tài)Gamma分布族，即對于給定的τ，參數(shù)θ的先驗分布為正態(tài)分布，精度參數(shù)τ的先驗分布為正態(tài)Gamma分布，即：

πθ，τ=πθ|τπτ（5）

其中，πθ，τ∝τ2exp-τ2θ-μTQθ-μ，πτ∝τα-1e-τβ。此處，超參數(shù)μ∈Rp，αgt;0，βgt;0。Q為p階正定矩陣。定義Y=y1，y2，…，ynT，根據(jù)貝葉斯定理，可得參數(shù)θ，τ的聯(lián)合后驗分布密度函數(shù)與參數(shù)先驗分布密度函數(shù)式（4）、模型擬合函數(shù)式（5）二者的乘積成正比，即：

πθ，τ|Y∝Lθ，τπθ，τ（6）

通過在區(qū)間0，SymboleB@上πθ，τ|Y對τ的積分，求出自回歸系數(shù)θ的后驗分布函數(shù)為：

πθ|Y=∫τgt;0πθ，τ|Ydτ（7）

類似地，將πθ，τ|Y對θ進行積分，則可得到參數(shù)τ的后驗分布密度函數(shù)為：

πτ|Y=∫θ∈Rpπθ，τ|Ydθ（8）

1.3"MCMC收斂性診斷

1.3.1"收斂性軌跡圖

軌跡圖（trace plot） [11]：根據(jù)迭代次數(shù)繪制生成的樣本路徑圖。在實際應(yīng)用中，為避免出現(xiàn)鏈條落入目標(biāo)局部分布區(qū)域，通常可以從不同的初始點出發(fā)同時生成多個馬爾可夫鏈。在經(jīng)過若干次迭代后，如果它們的樣本軌跡圖基本穩(wěn)定并重合，則可以認(rèn)為抽樣已經(jīng)收斂。

1.3.2"蒙特卡羅誤差

馬爾可夫鏈的收斂從初始點開始，隨著迭代次數(shù)的增加，波動逐漸減小。因此，通過參數(shù)或者其函數(shù)估計（即后驗樣本均值）的方差或標(biāo)準(zhǔn)差（稱為MC誤差）來評估對應(yīng)馬爾可夫鏈的收斂性[12]。較小的MC誤差可以提高參數(shù)估計的準(zhǔn)確性，通常情況下，模擬次數(shù)越多，蒙特卡羅誤差越小，模型的估計結(jié)果也就越準(zhǔn)確。

2"基于MCMC方法的貝葉斯統(tǒng)計模型應(yīng)用

2.1"基于MCMC方法的線性回歸模型

2.1.1"數(shù)據(jù)集

LearnBayes軟件包中的Birdextinct數(shù)據(jù)集包含了過去幾十年中在英國附近16個島嶼上62種鳥的四類數(shù)據(jù)：（1）在島上的平均滅絕時間（TIME）；（2）平均筑巢數(shù)（NESTING）；（3）種群規(guī)模（SIZE），分為“大”（用1表示）與“小”（用0表示）兩類；（4）棲息狀態(tài)（STATUS），分為“遷徙”（用1表示）與“久居”（用0表示）兩類。數(shù)據(jù)如表1所示，主要研究目的是找出該地區(qū)鳥類的滅絕時間與另外三個變量之間的關(guān)系。預(yù)測變量的顯著性分析如下：用y表示響應(yīng)變量平均滅絕時間；用x1表示預(yù)測變量平均筑巢數(shù)；用x2表示預(yù)測變量種群規(guī)模；用x3表示預(yù)測變量棲息狀態(tài)。由于前期分析中發(fā)現(xiàn)變量y嚴(yán)重右偏，因此對其進行對數(shù)變換。將此問題歸為線性回歸模型：

Log（y）=β0+β1x1+β2x2+β3x3+ε（9）

其中，β0、β1、β2、β3為回歸系數(shù)常數(shù)，ε為誤差項。

2.1.2"建模與收斂性

對Birdextinct數(shù)據(jù)集進行線性回歸建模，結(jié)果如表2所示。

根據(jù)模型結(jié)果，可得到的回歸方程為y=04306+0.2651x1-0.6526x2+0.5047x3。從表2可知，種群的滅絕時間與筑巢數(shù)的關(guān)系高度顯著，表明筑巢數(shù)越多，種群滅絕時間越長；而種群規(guī)模和棲息狀態(tài)與種群滅絕時間的相關(guān)顯著性較弱，表明規(guī)模越大的鳥類其滅絕時間越短，而遷徙的鳥類其滅絕時間較長，其模型的模擬軌跡如圖1所示。

從圖1中可知，平均筑巢數(shù)的軌跡圖中的路徑?jīng)]有明顯的周期性和趨勢性。通過計算各個迭代次數(shù)的平均筑巢數(shù)，得到均值約為0.25＜1.2[13]，因此認(rèn)為該馬爾可夫鏈?zhǔn)鞘諗康摹?/p>

2.2"基于MCMC的貝葉斯時間序列預(yù)測模型

2.2.1"CPI數(shù)據(jù)集來源

本研究從國家統(tǒng)計局網(wǎng)站上獲取1982-2021年居民消費價格指數(shù)（CPI），部分?jǐn)?shù)據(jù)如表3所示。在建立模型過程中，使用1982-2020年數(shù)據(jù)用作訓(xùn)練集，2021年數(shù)據(jù)作為測試數(shù)據(jù)。

從表3可知，居民消費指數(shù)從1982年到2020年有明顯的上升趨勢，原始序列是非平穩(wěn)的。

2.2.2"性能預(yù)測

對數(shù)據(jù)進行增項DF單位根檢驗（Augmented DickeyFuller test，ADF），ADF檢驗結(jié)果的p值為0999 0，顯著大于0.05。顯然，單位根檢驗沒有通過，證實原CPI序列是非平穩(wěn)序列。為了進行時間序列分析，需要對數(shù)據(jù)進行一階差分，即Yt= Zt-Zt-1，得到序列{Yt1} ，對其進行ADF單位根檢驗時，檢驗結(jié)果的p值為0.974 7，此時序列仍然非平穩(wěn)，因此需要進行二次差分以得到平穩(wěn)序列{Yt}。通過二次差分變換后時間數(shù)據(jù)序列間的自相關(guān)系數(shù)（ACF）和偏自相關(guān)系數(shù)（PACF）結(jié)果如圖2所示。

從圖2可知，中國居民消費指數(shù)數(shù)據(jù)在經(jīng)過二次差分變換后的序列具有1階自相關(guān)特征，得到適應(yīng)的模型為AR（1），即yt=β1yt-1+εt。選擇正態(tài)-Gamma分布族作為未知參數(shù)初始值的先驗分布，即β1～N（10-6，0.1），τ～Gamma（10-5，10-3），其中τ=σ-1。在吉布斯采樣模擬過程中，為了確保參數(shù)的收斂性，首先進行1 000次預(yù)迭代，然后丟棄預(yù)迭代的值，重新執(zhí)行99 000次迭代，可以得到β1的均值為-0.3423，95%置信區(qū)間為（-0.663 8，-0.021 0），結(jié)果如表4所示。

從表4可知，參數(shù)的蒙特卡羅誤差值均相對較小，因此可以認(rèn)為馬爾可夫鏈?zhǔn)諗俊８鶕?jù)參數(shù)估計，2021年CPI指數(shù)為2 334.4，與真實值2 352.8的誤差為0.78%。

2.2.3"收斂性診斷

將兩組不同初始值分別進行100 000次迭代分析，其中1 000次作為預(yù)迭代分析。如圖3所示，在迭代次數(shù)達(dá)20 000次后，β1和τ收斂性判斷圖中，兩者都接近為一條水平直線并重合，說明馬爾可夫鏈抽樣收斂。

3"結(jié)語

為研究MCMC 算法在貝葉斯統(tǒng)計模型中的應(yīng)用，文章運用birdextinct數(shù)據(jù)集構(gòu)建貝葉斯線性回歸模型，并對數(shù)據(jù)進行預(yù)處理，通過R語言程序?qū)?shù)據(jù)集進行收斂性診斷，研究鳥類的滅絕時間與其他三個量之間的關(guān)系。同時，選取1982-2021年居民消費價格指數(shù)（CPI）作為時間序列，結(jié)合差分處理數(shù)據(jù)建立AR（p）時間序列預(yù)測模型。根據(jù)參數(shù)估計得到2021年CPI數(shù)據(jù)為2 334.4，與真實值2 352.8的誤差為0.78%。通過貝葉斯線性回歸模型與AR（p）時間序列預(yù)測模型的案例分析，驗證了MCMC方法能有效解決貝葉斯統(tǒng)計問題。然而，文中使用的MCMC算法估計模型參數(shù)時主要采用了Gibbs抽樣方法，盡管該方法高效，但當(dāng)某些條件后驗分布不是標(biāo)準(zhǔn)分布時不適用。因此，后續(xù)需要對數(shù)據(jù)集進行進一步分析，選取合適的先驗信息分布和初始值，以便更好地應(yīng)用MCMC算法解決統(tǒng)計模型問題。

參考文獻：

[1]李貴玉，顧昕.貝葉斯統(tǒng)計方法的應(yīng)用與現(xiàn)狀[J].心理學(xué)探新，2021，41（05）：466473.

[2]孟光磊，叢澤林，宋彬，等.貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)綜述[J/OL].北京航空航天大學(xué)學(xué)報，124[20231129].https：//link.cnki.net/doi/10.13700/j.bh.10015965.2023.0445.

[3]谷恒明.經(jīng)典統(tǒng)計學(xué)與貝葉斯統(tǒng)計學(xué)在回歸模型中的比較研究[D].北京：軍事科學(xué)院，2018.

[4]梁瑛，季憲軍.小樣本合格率估計的隨機加權(quán)貝葉斯方法[J].統(tǒng)計與決策，2016（12）：46.

[5]黨紅.貝葉斯統(tǒng)計中單參數(shù)后驗分布的精確計算方法[J].長治學(xué)院學(xué)報，2019，36（02）：67.

[6]Metropolis N，Ulam S.The Monte Carlo method[J].Journal of the American Statistical Association，1949，44：335341.

[7]王緒民，鄭順超.基于BayesianMCMC算法的水利工程投標(biāo)報價分布預(yù)測[J].水電能源科學(xué)，2023，41（09）：155158，206.

[8]王若嬋，嚴(yán)剛，李揚，等.紊流激勵下基于貝葉斯統(tǒng)計推斷的顫振邊界預(yù)測方法研究[J].振動與沖擊， 2023，42（12）：283289.

[9]胡紅波.MCMC方法在測量不確定度評估中的應(yīng)用[J].計量技術(shù)，2020（05）：8994，88.

[10]王婧，鮑貴.貝葉斯統(tǒng)計與傳統(tǒng)統(tǒng)計方法的比較[J].統(tǒng)計與決策，2021，37（01）：2429.

[11]張航綺.利用蒙特卡羅法計算投資組合VaR：以嘉實新興產(chǎn)業(yè)股票（000751）為例[J].現(xiàn)代營銷（經(jīng)營版），2022（03）：9193.

[12]靳海強，文俊，王彤彤，等.基于改進方差特性二維混合抽樣蒙特卡羅法的多諧波源疊加研究[J].電網(wǎng)與清潔能源，2021，37（01）：1623，31.

[13]Geman A， Rubin DB.Inference for iterative simulation using multiple sequences[J].Statistics Science，1992，7（04）：457511.

責(zé)任編輯：肖祖銘

Research on the Application of Bayesian Statistical Model Based on MCMC Method

ZHANG Zongyu， XU Jun， JIANG Kui， CHEN Shichao

（School of Computer and Data Engineering， Bengbu College of Business and Industry， Bengbu 233000， China）

Abstract：To verify the convergence of Bayesian statistical posterior distribution calculation， two classic Bayesian statistical models， namely linear regression and time series， are investigated and constructed in this paper. Based on the Markov Chain Monte Carlo （MCMC） method， the Bayesian linear regression model is validated using the Birdextinct dataset to analyze the relationship between the average extinction time of birds and three variables： average number of nests， population size， and habitat status. Furthermore， based on the series of Consumer Price Index （CPI） from 1982 to 2021， an AR （p） time series forecasting model is established， incorporating differential processing data. The results indicate that both the Bayesian linear regression model and the Bayesian time series forecasting model have converged Markov chains， with the error of only 0.78% for the Bayesian time series forecasting model. The MCMC method can be effectively applied to analysis of Bayesian statistical models.

Keywords： MCMC method; Bayesian statistics; convergence diagnosis; time series forecasting