唐愛萍

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

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

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

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

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)個個體的線性MESSM

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

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

對于線性 MESSM，設(shè) Yit?{yi1，yi2，…，yit}表示第i個個體截止時刻t的所有觀測數(shù)據(jù)的集合，此部分討論參數(shù)已知，給定觀測數(shù)據(jù){Yit}mi=1時，通過卡爾曼濾波與蒙特卡洛方法的結(jié)合實(shí)現(xiàn)對MESSM中狀態(tài)的估計問題，在此重點(diǎn)討論狀態(tài)的一步向前預(yù)測估計問題，其他情況均可由濾波算法推導(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］計算得到。若從 f(θi中進(jìn)行抽樣得到樣本集{θ(ij)}Nj=1，由

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

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

(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)，對于每個參數(shù) θ0(j)一步向前預(yù)測

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

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

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

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

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

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

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

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

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

［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.