基于gAIC和滑動(dòng)窗的自回歸模型參數(shù)估計(jì)算法

楊淳沨 吳國(guó)成 伍家松 姜龍玉 孔佑勇 舒華忠

(東南大學(xué)計(jì)算機(jī)網(wǎng)絡(luò)和信息集成教育部重點(diǎn)實(shí)驗(yàn)室， 南京 210096)(東南大學(xué)計(jì)算機(jī)科學(xué)與工程學(xué)院， 南京 210096)(東南大學(xué)中法生物醫(yī)學(xué)信息研究中心, 南京 210096)

近年來(lái)，通過(guò)分析和處理大腦電信號(hào)來(lái)研究大腦連通性已成為腦工作原理和腦部重大疾病發(fā)生機(jī)制研究領(lǐng)域的熱點(diǎn)之一[1-2].WGCI (Wiener-Granger causality index)算法[3]、DTF (directed transfer function)算法[4]、PDC (partial directed coherence)算法[5]等都可有效應(yīng)用于研究大腦功能連通性和效應(yīng)連通性[6].這些算法都是基于隨機(jī)過(guò)程的自回歸(autoregressive, AR)模型的算法，模型參數(shù)(階數(shù)和系數(shù))的估計(jì)起著舉足輕重的作用.傳統(tǒng)的AR模型階數(shù)估計(jì)算法包括赤池信息準(zhǔn)則(Akaike information criterion, AIC)[7]、最后預(yù)測(cè)誤差法(final prediction error, FPE)和最小描述長(zhǎng)度法(minimum description length, MDL)[8]等.這些算法均假設(shè)不同信號(hào)具有相同階數(shù).然而，在實(shí)際情況中，不同信號(hào)各自的階數(shù)往往各不相同，因此這些算法會(huì)使模型項(xiàng)數(shù)過(guò)估計(jì)，從而導(dǎo)致后續(xù)的系數(shù)估計(jì)運(yùn)算量加大且精準(zhǔn)度不高.學(xué)者們繼而開(kāi)發(fā)了更加有效的AR模型參數(shù)估計(jì)算法[9-11].其中，最佳參數(shù)搜索(optimal parameter search, OPS)算法[12]是一種魯棒性較高的算法，即使在存在噪聲且模型初始階數(shù)選取不當(dāng)?shù)那闆r下，OPS算法仍能篩選出正確的模型項(xiàng).但是，在干擾項(xiàng)較多的情況下，OPS算法的性能下降較快，從而會(huì)產(chǎn)生較大誤差，導(dǎo)致難以篩選出所有正確的模型項(xiàng).

本文提出了一種基于gAIC (generalized AIC)[13]和滑動(dòng)窗的自回歸模型參數(shù)估計(jì)算法.其中，模型的階數(shù)估計(jì)采用gAIC算法，與AIC相比，它能更有效地減少干擾項(xiàng)數(shù)目.對(duì)原始信號(hào)進(jìn)行加窗處理，可以消除部分信號(hào)的不穩(wěn)定性.實(shí)驗(yàn)結(jié)果表明，相比傳統(tǒng)算法，所提算法在模型參數(shù)估計(jì)方面具有更好的魯棒性.

1 模型參數(shù)估計(jì)算法

對(duì)于一個(gè)平穩(wěn)時(shí)間序列y,y(n)表示其在n時(shí)刻的值.如果已知時(shí)間序列y在某一時(shí)刻n之前p個(gè)時(shí)刻的值，可以建立一個(gè)基于該時(shí)間序列過(guò)去p個(gè)值的自回歸模型，并且可用該模型來(lái)預(yù)測(cè)該時(shí)間序列將來(lái)的值，即

(1)

式中，p為該線性自回歸模型的階數(shù)；ai為y(n-i)項(xiàng)的系數(shù)；e(n)為真實(shí)值和預(yù)測(cè)值之間的誤差項(xiàng).

若在單變量的基礎(chǔ)上引入另一個(gè)平穩(wěn)時(shí)間序列x，可建立二維自回歸模型，即

(2)

1.1 自回歸模型的模型階數(shù)估計(jì)算法

模型階數(shù)估計(jì)是自回歸模型參數(shù)估計(jì)過(guò)程中的重要步驟.赤池信息準(zhǔn)則[7]是常用的模型階數(shù)估計(jì)算法，它是一種基于信息論的算法，主要思想是在模型與數(shù)據(jù)的擬合精準(zhǔn)度和模型復(fù)雜度之間進(jìn)行權(quán)衡，得出模型的最優(yōu)階數(shù).

給定一個(gè)零均值的m維向量自回歸過(guò)程zn={zn(1),zn(2),…,zn(m)}T，即

(3)

式中，Φi(i=1,2,…,p)為m×m的系數(shù)矩陣；wn為零均值獨(dú)立同分布的向量.假定對(duì)于任意的變量i>0，wn和zn-i是相互獨(dú)立的.

假設(shè)觀測(cè)值z(mì)1,z2,…,zN都是由式(3)中的自回歸過(guò)程所產(chǎn)生，其中N為序列的長(zhǎng)度.為了求得模型階數(shù)p的估計(jì)值q，設(shè)定一個(gè)估計(jì)階數(shù)的最大上限K，令q取1～K中所有整數(shù)值.對(duì)于每一個(gè)q值，均采用最小二乘法對(duì)原模型進(jìn)行估計(jì)，得到如下模型：

(4)

ZAj=zq+1,jj=1,2,…,m

(5)

式中

式中，Φi(j)表示Φi的第j行.

利用最小二乘法求解式(5)可得

(6)

原模型的預(yù)測(cè)誤差協(xié)方差矩陣可表示為

(7)

采用AIC定義式選擇最優(yōu)的階數(shù)估計(jì)q，即

(8)

令q取遍1～K中所有整數(shù)值，當(dāng)AIC(q)值最小時(shí)，對(duì)應(yīng)的q值即為該模型的最優(yōu)階數(shù).

對(duì)于一個(gè)二維自回歸模型(見(jiàn)式(2))，由式(8)可知，AIC算法是在假定自回歸模型中不同信號(hào)具有相同階數(shù)的基礎(chǔ)上進(jìn)行估計(jì)的，即q=qxx=qyx=qxy=qyy.然而，在實(shí)際情況中，自回歸模型中不同信號(hào)的階數(shù)往往各不相同，AIC算法會(huì)在估計(jì)模型階數(shù)時(shí)產(chǎn)生過(guò)估計(jì)，導(dǎo)致后續(xù)系數(shù)估計(jì)運(yùn)算量加大且精準(zhǔn)度不高.針對(duì)這一問(wèn)題，本文采用另一種改進(jìn)的模型階數(shù)估計(jì)算法gAIC[13],其定義如下：

(9)

1.2 OPS算法

以式(2)中的時(shí)間序列y為例，可得

(10)

采用OPS算法時(shí)，首先從式(10)模型項(xiàng)所組成的候選向量中篩選出線性無(wú)關(guān)項(xiàng).這些候選向量所構(gòu)成的矩陣表示為

[V1V2…Vqyy+qxy]T

(11)

式(11)中每一行代表一個(gè)候選向量.選擇V1作為第1個(gè)選出的向量，判定下一個(gè)候選向量V2和V1是否線性無(wú)關(guān)，即利用V1來(lái)擬合V2，并計(jì)算V2和估計(jì)向量的誤差值.一般情況下，在無(wú)噪聲的信號(hào)中，向量之間應(yīng)該是線性無(wú)關(guān)的.但是在有噪聲的情況 (即誤差不為零) 下，可以設(shè)置一個(gè)最大誤差值，當(dāng)誤差值小于該最大值時(shí)，向量V2便可被選為線性無(wú)關(guān)的向量，反之則丟棄向量V2.選出V2作為極大線性無(wú)關(guān)組中的另一個(gè)候選向量后，采用向量V1和V2來(lái)估計(jì)V3的線性無(wú)關(guān)性.重復(fù)以上步驟直到選出一個(gè)新的候選向量組合φ=[ω1ω2…ωR]，其中R為所選出的線性無(wú)關(guān)向量總數(shù).

篩選出新的線性無(wú)關(guān)的候選向量后，對(duì)其進(jìn)行最小二乘法處理，即將式(10)改寫(xiě)為

(12)

式中

θg={g1,g2,…,gR}

(13)

式中，gi(i=1,2,…,R)為模型系數(shù).為了最小化誤差ey(t)，令下式達(dá)到最小值：

(14)

然后利用最小二乘法可得

(15)

1.3 自適應(yīng)OPS算法

對(duì)于1.2節(jié)中的OPS算法而言，確定候選向量集合，即式(11)中行向量所形成的矩陣，是一個(gè)非常重要的過(guò)程.AIC準(zhǔn)則存在過(guò)度擬合的問(wèn)題，且對(duì)二維自回歸模型的每個(gè)維度得到相同的階數(shù)估計(jì)值，因此，本文采用gAIC算法對(duì)自回歸模型進(jìn)行階數(shù)估計(jì)，得到每個(gè)維度的各自階數(shù)，繼而進(jìn)行后續(xù)計(jì)算.另一方面，考慮到OPS算法對(duì)于噪聲敏感，為了減少噪聲的影響，對(duì)原始信號(hào)進(jìn)行重疊加窗處理.以式(10)為例，對(duì)于一個(gè)長(zhǎng)度為N的信號(hào)，選擇一個(gè)寬度為W的滑動(dòng)窗，其滑動(dòng)步長(zhǎng)L，具體步驟如下：

①?gòu)男盘?hào)中提取長(zhǎng)度為W的信號(hào)，記為Wi；

②對(duì)每段窗信號(hào)Wi，利用gAIC估計(jì)出模型階數(shù)，得到式(11)中的矩陣V，進(jìn)而使用OPS算法計(jì)算矩陣V中每個(gè)候選項(xiàng)的權(quán)重值；

③將窗向后滑動(dòng)一段距離L；

⑤計(jì)算式(4)中T個(gè)窗信號(hào)每個(gè)候選項(xiàng)權(quán)重的均值，并以此作為該候選項(xiàng)的權(quán)重值，對(duì)候選項(xiàng)進(jìn)行最終篩選.

2 實(shí)驗(yàn)結(jié)果及分析

實(shí)驗(yàn)中采用了如下的二維自回歸模型：

(16)

式中，ex(n)和ey(n)表示均值為0、方差為1的相互獨(dú)立的高斯白噪聲序列.

為比較OPS算法和自適應(yīng)OPS算法的效果，進(jìn)行了1 000組實(shí)驗(yàn)，并將結(jié)果應(yīng)用于式(16)y(n)的參數(shù)估計(jì)中.在OPS算法中，為了考察模型階數(shù)過(guò)大情況下算法的性能，需預(yù)設(shè)一個(gè)比實(shí)際情況大的階數(shù)上限.針對(duì)本文實(shí)驗(yàn)中所用模型，該階數(shù)上限設(shè)為10.由1.2節(jié)可知，模型階數(shù)越大，候選向量越多，從而導(dǎo)致計(jì)算量越大.針對(duì)這一問(wèn)題，分別采用AIC和gAIC算法進(jìn)行階數(shù)估計(jì)，以減少不必要的計(jì)算.由AIC得到的估計(jì)結(jié)果為6，即候選項(xiàng)為Y(n-1),Y(n-2), …,Y(n-6),X(n-1),X(n-2), …,X(n-6)，從而將候選向量個(gè)數(shù)由20降至12.采用gAIC算法對(duì)模型階數(shù)進(jìn)行估計(jì)，得到的階數(shù)qyy=6，qxy=3，即候選項(xiàng)為Y(n-1),Y(n-2), …,Y(n-6),X(n-1),X(n-2),X(n-3).然后，考慮加窗處理對(duì)OPS算法估計(jì)模型階數(shù)正確率的影響.由此設(shè)計(jì)了6種組合實(shí)驗(yàn)方案，見(jiàn)表1.

表1 6種組合實(shí)驗(yàn)方案

下面以表1中第6種組合實(shí)驗(yàn)方案為例，說(shuō)明模型階數(shù)的選取過(guò)程.首先，由gAIC算法得到信號(hào)y(n)中的模型階數(shù)qyy=6，qxy=3.然后，由第1節(jié)中所述算法求得每個(gè)候選項(xiàng)的權(quán)重值，結(jié)果見(jiàn)圖1.圖中，候選項(xiàng)的權(quán)重值表示每個(gè)候選項(xiàng)的重要程度，按遞減順序排列.由圖可知，X(n-3)項(xiàng)的權(quán)重值最大.按照1.2節(jié)所述篩選候選項(xiàng)的原則，權(quán)重值前四大項(xiàng)(即X(n-3)，Y(n-6)，X(n-1)，Y(n-5))相加之和已超出預(yù)先設(shè)定的閾值0.92，這4項(xiàng)即為最終需保留的候選項(xiàng).其余5項(xiàng)的權(quán)重值較小，接近于0，表明其對(duì)信號(hào)Y的貢獻(xiàn)可以忽略不計(jì)，該結(jié)果與式(16)中模型一致.

圖1 各候選項(xiàng)的權(quán)重值

為了驗(yàn)證算法對(duì)不同長(zhǎng)度信號(hào)的有效性，對(duì)于表1中的6種方案，每種方案選取2種不同長(zhǎng)度(1 000和500)的信號(hào)，分別進(jìn)行1 000次實(shí)驗(yàn).對(duì)于長(zhǎng)度為1 000的信號(hào)，滑動(dòng)窗長(zhǎng)度為800，步長(zhǎng)為40；對(duì)于長(zhǎng)度為500的信號(hào)，滑動(dòng)窗長(zhǎng)度為300，步長(zhǎng)為40.考慮到OPS算法需要閾值，而閾值選取沒(méi)有統(tǒng)一標(biāo)準(zhǔn)，因此，本文采用不同閾值進(jìn)行實(shí)驗(yàn).實(shí)驗(yàn)結(jié)果見(jiàn)圖2.實(shí)驗(yàn)結(jié)果中的正確定義為算法篩選出完全正確的4個(gè)候選項(xiàng)，多選和少選都?xì)w為錯(cuò)誤.由圖可知，大部分閾值范圍(除高閾值范圍)內(nèi)，有窗情況都明顯優(yōu)于無(wú)窗情況；在同樣有窗/無(wú)窗情況下，gAIC算法的正確率均高于其他2種算法.N=1 000時(shí)最佳閾值約為0.92，N=500時(shí)最佳閾值約為0.88，正確率均接近90%.當(dāng)信號(hào)長(zhǎng)度從1 000減少到500時(shí)，OPS算法仍具有較好的性能，如gAIC加窗方案，最佳閾值不同，但正確率并未明顯下降，均接近于90%.

(a) N=1 000

(b) N=500圖2 不同閾值下6種方案結(jié)果對(duì)比

圖1中，Y(n-5)項(xiàng)之后候選項(xiàng)的投影值較前4項(xiàng)大幅減少，可以認(rèn)為候選項(xiàng)集合被分成2個(gè)部分，即圖中前4項(xiàng)和后5項(xiàng).由圖2可知，閾值較大時(shí)，算法的正確率急劇降低至0，根據(jù)OPS算法，當(dāng)閾值慢慢增大直到1時(shí)，包含的候選項(xiàng)增加，即逐漸將剩余的權(quán)重值較小的候選項(xiàng)全部選入.由此可知，閾值的選取對(duì)算法的正確率起著至關(guān)重要的作用.因此，采用了1.3節(jié)中的自適應(yīng)OPS算法選取最優(yōu)閾值.對(duì)于圖1中按權(quán)重值降序排列的9個(gè)候選項(xiàng)，依次計(jì)算相鄰2項(xiàng)的比值，結(jié)果見(jiàn)圖3.由圖可知，比值為Y(n-5)/X(n-2)時(shí)權(quán)重值取得最大值，其后項(xiàng)的權(quán)重值顯著降低.因此，僅保留前面4項(xiàng).

自適應(yīng)OPS算法的實(shí)驗(yàn)結(jié)果見(jiàn)圖4.實(shí)驗(yàn)中信號(hào)長(zhǎng)度為1 000.由圖可知，gAIC有窗的實(shí)驗(yàn)方案取得了最高的正確率，這與圖2(a)中的結(jié)果一致.

圖3 相鄰2項(xiàng)權(quán)重值的比值

綜上可知，通過(guò)在算法中引入gAIC，可以有效減少實(shí)驗(yàn)過(guò)程中候選項(xiàng)集合里干擾項(xiàng)的個(gè)數(shù)；采用加窗的方式，信號(hào)能更加平滑，從而減少噪聲帶來(lái)的影響.這2種改進(jìn)均能使OPS算法的正確率明顯提升.自適應(yīng)OPS算法的正確率較gAIC有窗方案有所下降，但其計(jì)算過(guò)程中不需要閾值，尤其適用于沒(méi)有先驗(yàn)知識(shí)的數(shù)據(jù)或模型.

圖4 自適應(yīng)OPS算法的實(shí)驗(yàn)結(jié)果

3 結(jié)論

1) 針對(duì)自回歸模型參數(shù)估計(jì)問(wèn)題，本文提出了一種基于gAIC和滑動(dòng)窗的模型參數(shù)估計(jì)算法.

2) 在模型參數(shù)估計(jì)中，使用gAIC算法對(duì)模型階數(shù)進(jìn)行估計(jì)，從而有效減少模型中候選項(xiàng)集合里干擾項(xiàng)的個(gè)數(shù).然而，對(duì)信號(hào)進(jìn)行加窗處理，使得信號(hào)能更加平滑，減少噪聲帶來(lái)的影響.最后，采用自適應(yīng)OPS算法，進(jìn)一步篩選候選項(xiàng)集合中的候選項(xiàng)，得到最終的模型選項(xiàng)及相應(yīng)的模型系數(shù)值.

3) 實(shí)驗(yàn)結(jié)果表明，對(duì)于2種不同的實(shí)驗(yàn)信號(hào)長(zhǎng)度，6種組合實(shí)驗(yàn)方案下所提算法都表現(xiàn)出最優(yōu)的魯棒性，其正確率接近于90%.

參考文獻(xiàn)(References)

[1] 蒲慕明, 徐波, 譚鐵牛. 腦科學(xué)與類腦研究概述[J]. 中國(guó)科學(xué)院院刊, 2016, 31(7): 714,725-736. DOI: 10.16418/j.issn.1000-3045.2016.07.001.

[2] Sporns O. Brain connectivity [EB/OL]. (2007-10-27)[2017-05-10]. http://scholarpedia.org/article/Brain_connectivity.

[3] Granger C W J. Investigating causal relations by econometric models and cross-spectral methods [J].Econometrica, 1969,37(3): 424-438. DOI: 10.2307/1912791.

[4] Kamiński M J, Blinowska K J. A new method of the description of the information flow in the brain structures[J].BiologicalCybernetics, 1991,65(3): 203-210. DOI: 10.1007/BF00198091.

[5] Baccalá L A, Sameshima K. Partial directed coherence: A new concept in neural structure determination[J].BiologicalCybernetics, 2001,84(6): 463-474. DOI: 10.1007/PL00007990.

[6] Friston K J. Functional and effective connectivity:A review[J].BrainConnectivity, 2011,1(1): 13-36. DOI: 10.1089/brain.2011.0008.

[7] Akaike H. Power spectrum estimation through autoregressive model fitting[J].AnnalsoftheInstituteofStatisticalMathematics, 1969,21(1): 407-419. DOI: 10.1007/BF02532269.

[8] Rissanen J. A universal prior for integers and estimation by minimum description length[J].TheAnnalsofstatistics, 1983,11(2): 416-431. DOI:10.1214/aos/1176346150.

[9] Khorshidi S, Karimi M, Nematollahi A R. New autoregressive (AR) order selection criteria based on the prediction error estimation[J].SignalProcessing, 2011,91(10): 2359-2370. DOI:10.1016/j.sigpro.2011.04.021.

[10] Zhao Y, Billings S A, Wei H L, et al. A parametric method to measure time-varying linear and nonlinear causality with applications to EEG data[J].IEEETransactionsonBiomedicalEngineering, 2013,60(11): 3141-3148. DOI: 10.1109/TBME.2013.2269766.

[11] Stoica P, Babu P. Model order estimation via penalizing adaptively the likelihood (PAL)[J].SignalProcessing, 2013,93(11): 2865-2871. DOI: 10.1016/j.sigpro.2013.03.014.

[12] Lu S, Ju K H, Chon K H. A new algorithm for linear and nonlinear ARMA model parameter estimation using affine geometry[J].IEEETransactionsonBiomedicalEngineering, 2001,48(10): 1116-1124. DOI: 10.1109/ 10.951514.

[13] Yang C, Le Bouquin Jeannes R, Bellanger J J, et al. A new strategy for model order identification and its application to transfer entropy for EEG signals analysis[J].IEEETransactionsonBiomedicalEngineering, 2013,60(5): 1318-1327. DOI: 10.1109/TBME.2012.2234125.