基于蒙特卡洛方法的線性混合效應(yīng)狀態(tài)空間模型

唐愛萍

(西安電子科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，陜西西安 710071)

狀態(tài)空間模型［1］是一個(gè)重要的分析模型，在經(jīng)濟(jì)、工程、醫(yī)學(xué)等領(lǐng)域發(fā)揮著重要作用。而線性狀態(tài)空間模型是應(yīng)用廣泛的狀態(tài)空間模型，尤其是在時(shí)間序列及控制數(shù)據(jù)中，較多數(shù)據(jù)模型如ARMA模型等均可轉(zhuǎn)化成線性狀態(tài)空間模型。對(duì)于參數(shù)已知的高斯線性狀態(tài)空間模型，卡爾曼濾波［2］是一種非常有效的狀態(tài)估計(jì)算法。除此之外，文獻(xiàn)［3］首先提出了運(yùn)用序貫蒙特卡洛的方法解決狀態(tài)的估計(jì)問題。

縱向數(shù)據(jù)［4－5］的應(yīng)用廣泛，將其引入到線性狀態(tài)空間模型中具有較大應(yīng)用價(jià)值，但需克服的困難是如何處理高維數(shù)據(jù)以及同時(shí)實(shí)現(xiàn)每個(gè)個(gè)體的狀態(tài)估計(jì)問題。文獻(xiàn)［6］將個(gè)體參數(shù)當(dāng)做隱變量，引入到混合效應(yīng)方程［7］中，從而形成對(duì)混合效應(yīng)狀態(tài)空間模型(Mixed － Effects State Space Model，MESSM)的研究［7］，加入隱變量到線性狀態(tài)空間模型后，此時(shí)模型變?yōu)榉蔷€性狀態(tài)空間模型。對(duì)于非線性狀態(tài)空間模型［8－9］，目前國內(nèi)外已提出的方法有擴(kuò)展卡爾曼濾波等算法，但卻均不是最優(yōu)的解決非線性狀態(tài)空間模型的最優(yōu)算法，且始終沒有新的研究進(jìn)展。而對(duì)于非線性混合效應(yīng)狀態(tài)模型的研究則更少。

文獻(xiàn)［6］分別將狀態(tài)變量分為總體和個(gè)體的來研究，但前提是要應(yīng)用預(yù)測(cè)的個(gè)體參數(shù)，這增加了狀態(tài)估計(jì)的成本。因此對(duì)于縱向數(shù)據(jù)的狀態(tài)空間模型，研究出一種無需直接估計(jì)個(gè)體參數(shù)而進(jìn)行狀態(tài)估計(jì)的方法，將進(jìn)一步提高M(jìn)EMSE的應(yīng)用價(jià)值。

1 線性混合效應(yīng)狀態(tài)空間模型

線性混合效應(yīng)狀態(tài)空間模型是狀態(tài)空間模型的擴(kuò)展，其是在狀態(tài)空間模型式(1)和式(2)基礎(chǔ)上增加了隨機(jī)效應(yīng)方程(3)?；诒疚囊芯康目v向數(shù)據(jù)，給出第i(i=1，2，…，n)個(gè)個(gè)體的線性MESSM

其中，xit，yit分別表示時(shí)刻t第i個(gè)個(gè)體的p×1狀態(tài)變量和q×1觀測(cè)變量;H(θi)，G(θi)分別是 p×p狀態(tài)轉(zhuǎn)移矩陣和q×p觀測(cè)矩陣，且含有未知參數(shù)θi;vit，wit分別是維數(shù)為p×1，q×1的狀態(tài)噪聲和觀測(cè)擾動(dòng)，當(dāng)i=j且 t=t'時(shí)，Cov(vit，vjt')=Q，Cov(wit，wjt')=R，其余情況均等于0。式(3)中eθ是固定效應(yīng)，bi是隨機(jī)誤差。對(duì)于正態(tài)的線性狀態(tài)空間模型，卡爾曼濾波算法可實(shí)現(xiàn)狀態(tài)xt最小無偏的在線估計(jì)。

2 基于序貫蒙特卡洛算法的狀態(tài)估計(jì)

對(duì)于線性 MESSM，設(shè) Yit?{yi1，yi2，…，yit}表示第i個(gè)個(gè)體截止時(shí)刻t的所有觀測(cè)數(shù)據(jù)的集合，此部分討論參數(shù)已知，給定觀測(cè)數(shù)據(jù){Yit}mi=1時(shí)，通過卡爾曼濾波與蒙特卡洛方法的結(jié)合實(shí)現(xiàn)對(duì)MESSM中狀態(tài)的估計(jì)問題，在此重點(diǎn)討論狀態(tài)的一步向前預(yù)測(cè)估計(jì)問題，其他情況均可由濾波算法推導(dǎo)。

f(xi，t+1是狀態(tài)變量xi，t+1的條件概率密度函數(shù)，式(4)可進(jìn)一步的寫為

其中，f(θiYit)是θi的后驗(yàn)分布。進(jìn)一步，由條件方差公式可得 xi，t+1的方差滿足

則由式(5)和式(6)可得出，只要給定θi，條件期望和條件方差可由卡爾曼濾波公式［1］計(jì)算得到。若從 f(θi中進(jìn)行抽樣得到樣本集{θ(ij)}Nj=1，由

對(duì)于每一個(gè) θ(j)i(j=1，2，…，N)，由卡爾曼濾波迭代公式可得

由以上分析可得出，對(duì)于線性MESSM，在隨機(jī)效應(yīng)θi未知時(shí)，將序貫蒙特卡洛重要性抽樣方法與卡爾曼濾波算法相結(jié)合，可實(shí)現(xiàn)對(duì)狀態(tài)的在線估計(jì)問題。具體做法如下:

(1)初始 t=0，假設(shè)初始狀態(tài)x0服從正態(tài)分布N(a0，P0)，其中 a0，P0已知。從先驗(yàn)分布 g(θi)中產(chǎn)生隨機(jī)樣本集{θ0(j)}Nj=1，且參數(shù) θ0(j)有權(quán)重 ω0(j)=，(j=1，2，…，N)，對(duì)于每個(gè)參數(shù) θ0(j)一步向前預(yù)測(cè)

(2)時(shí)刻 t=1，對(duì)于每個(gè)參數(shù) θ0(j)，用 Ω(tj)=表示相應(yīng)參數(shù)θ(tj)的卡爾曼濾波估計(jì)值。假設(shè)在 t－1 時(shí)刻，有觀測(cè)數(shù)據(jù) Yi，t－1，有粒子集{θ(t－j)1}Nj=1服從后驗(yàn)分布，則對(duì)于每個(gè) θ(t－j)1，對(duì)應(yīng)時(shí)刻的狀態(tài)估計(jì)值 Ω(t－j)1也可得到，從而由式(6)和式(7)可得

(3)在時(shí)刻 t=2，3，…，當(dāng)?shù)玫叫碌挠^測(cè)數(shù)據(jù) Yi，t=(yi1，…，yit)，對(duì)于 θ(t－j)1更 新 Ω(tj)。更新權(quán)重 w(tj)=然后從集合以概率產(chǎn)生隨機(jī)樣本

3 仿真數(shù)據(jù)的研究

為進(jìn)一步驗(yàn)證上文提到算法的有效性，對(duì)線性MSE模型進(jìn)行仿真實(shí)驗(yàn)。對(duì)于式(1)～式(3)，H(θi)=θi，G(θi)=1 設(shè)參數(shù)已知，取 eθ=0.8，D=0.006 4，Q=1，R=0.8，m=50，時(shí)間分別取T=20 和 T=40，從而產(chǎn)生了兩組模擬數(shù)據(jù)。為便于比較，在個(gè)體隨機(jī)效應(yīng)θi已知情況下運(yùn)用卡爾曼濾波算法進(jìn)行狀態(tài)的估計(jì)，而在其未知時(shí)運(yùn)用序貫蒙特卡洛濾波算法做狀態(tài)的估計(jì)。在運(yùn)用蒙特卡洛濾波算法進(jìn)行狀態(tài)估計(jì)時(shí)，對(duì)于每組模擬數(shù)據(jù)，分別取3組抽樣樣本 N=100，200，500。最終仿真結(jié)果如表1、圖1和圖2所示。

在表1中分別發(fā)給了在不同時(shí)間及不同采樣數(shù)下，兩種算法的均方誤差的比較數(shù)據(jù);在圖1和圖2中給出了兩種算法狀態(tài)估計(jì)值以及估計(jì)根均方誤差(RMSE)的比較結(jié)果，圖中標(biāo)注的ss(t)表示模擬出的真實(shí)狀態(tài)值，X1(t)，kf1(t)分別表示用蒙特卡洛濾波方法和卡爾曼濾波算法得到的狀態(tài)估計(jì)值，RMSE1，RMSE2分別表示兩種方法的RMSE。

表1 蒙特卡洛濾波和卡爾曼濾波算法均方誤差的比較

圖1 T=20時(shí)蒙特卡洛濾波與卡爾曼濾波結(jié)果比較

圖2 T=20時(shí)蒙特卡洛濾波與卡爾曼濾波RMSE

從表1、圖1和圖2的結(jié)果可看出，在同樣的數(shù)據(jù)下，兩種方法的估計(jì)結(jié)果相差較小。因此，對(duì)于混合效應(yīng)的狀態(tài)空間模型，在隨機(jī)效應(yīng)未知的情況下，本文提到的蒙特卡洛濾波算法是一種有效的狀態(tài)估計(jì)方法。本例中用到的線性模型雖較為簡單，但也容易拓展到高維線性模型中。同時(shí)，對(duì)于非線性狀態(tài)空間模型，只要可運(yùn)用線性化技術(shù)轉(zhuǎn)化為MESSM，在個(gè)體參數(shù)未知的情況下，本文中提到的方法也可以做狀態(tài)的在線估計(jì)問題。

［1］James Durbin.Times series analysis by state space methods［M］.US:Oxford University Press，2001.

［2］Harold W Sorenson.Kalman filtering:theory and application［M］.USA:IEEE Press，1985.

［3］Jun S Liu，Rong Chen.Sequential monte carlo methods for dynamic systems［J］.Journal of the American Association，1998，93(443):1031 －1044

［4］Mun J，Lindstrom M J.Diagnostics for repeated measurements in linear mixed effects models［J］.Statistics in Medicine，2013，32(8):1361 －1375.

［5］Liu D，Lu T，Niu X，et al.Mixed － effects state － space models for analysis of longitudinal dynamic systems［J］.Biometrics，2011，67(2):476 －485.

［6］Jun S Liu.Monte carlo strategies in scientific computing［M］.Germany:Springer，2002.

［7］Shayler Searle，George Casella.Variance components［M］.America:Wiley － Interscience，2011.

［8］Pedersen M W，Berg C W，Thygesen U H.Estimation methods for nonlinear state － space models in ecology［J］.Ecological Modelling，2011，222(8):1394 －1400.

［9］Zhou J，Han L，Liu SY.Nonlinear mixed － effects state space models with applications to HIV dynamic［J］.Statistics and Probability Letters，2013，83(5):1448 －1456.

［10］Huang Y X，Liu D C，Wu H L.Hierarchical bayesian methods for estimation of parameters in a longitudinal HIV dynamic system［J］.Biometrics，2006，62(2):413 －423.

［11］Markus Hurzeler，Hans R Kunschk，Monte Carlo.Approximations for general state － space models［J］.Journal of Computational and Graphical Statistics，1998，7(2):175 －193.

［12］Jun SLiu.Metropolized independent sampling with comparisons to rejection sampling and important sampling［J］.Statistics and Computing，1996，6(2):113 －119.