成立花，張俊敏

（1.西安工程大學(xué) 理學(xué)院，陜西 西安710048；2.西安建筑科技大學(xué) 理學(xué)院，陜西 西安710055）

最小二乘法（LS）和遞推的最小二乘法（RLS）是系統(tǒng)辨識(shí)中重要的算法.然而，在實(shí)際應(yīng)用中，由于各種干擾因素的存在，辨識(shí)所需的數(shù)據(jù)往往會(huì)被污染，包含數(shù)量和類型均未知的異常點(diǎn).這些異常點(diǎn)使估計(jì)精度降低，甚至完全失效.解決這一問題通常有兩種方式：從數(shù)據(jù)中找出異常點(diǎn)并剔除掉，用剩余的數(shù)據(jù)進(jìn)行辨識(shí)；發(fā)展替代LS和RLS的魯棒回歸算法.由于后者簡(jiǎn)便實(shí)用，備受學(xué)者關(guān)注，已發(fā)展出了許多算法［1－9］.在自回歸模型的參數(shù)估計(jì)中，數(shù)據(jù)中的異常點(diǎn)通常被分為兩個(gè)基本類［10］：加性異常點(diǎn)（第一型的異常點(diǎn)）與革新異常點(diǎn)（第二型的異常點(diǎn)）.革新異常點(diǎn)在系統(tǒng)運(yùn)行過程中遵循了真實(shí)系統(tǒng)的關(guān)系，因此，就系統(tǒng)的辨識(shí)而言，其不利影響較小.然而，當(dāng)數(shù)據(jù)中包含較大量的加性異常點(diǎn)型杠桿點(diǎn)時(shí)，這些算法的估計(jì)性能會(huì)嚴(yán)重下降，甚至失效.針對(duì)廣義極大似然類（GM）估計(jì)器中存在的問題，本文提出新的遞推有界GM 回歸估計(jì)算法.

1 問題描述

考慮模型yi＝xiβ＋vi，i＝1，2，…，n.其中：β＝［b1，b2，…，bp］T；yi是第i時(shí)刻響應(yīng)變量觀測(cè)值；x＝［xi，1，xi，2，…，xi，p］；vi是相互獨(dú)立，且尺度參數(shù)σ未知的干擾項(xiàng).所用的數(shù)據(jù)記作（Xn，Yn）.其中：Xn＝［xT2，xT2，…，xTn］T，Y＝［y1，y2，…，yn］T.獲得第n時(shí)刻的參數(shù)估計(jì)后，第i時(shí)刻相應(yīng)的殘差記作rn，i＝y(tǒng)i－.而在迭代法中，相應(yīng)于起始估計(jì)在第i時(shí)刻的殘差記作

采用文獻(xiàn)［9］中的一般框架，并加以改進(jìn)的風(fēng)險(xiǎn)函數(shù)為

式（1）中：λ為遺忘因子，0＜λ≤1；d（xi），ξ（xi）為待定函數(shù)；ρ（·）是一個(gè)改進(jìn)于Huber函數(shù)的有界M－估計(jì)函數(shù)，即

上式中：M為常值參數(shù).新估計(jì)器是式（1）的解.

2 新的GM 估計(jì)器

相應(yīng)于式（1）的“正規(guī)方程”為

式（2）中：ψ＝ρ′.

定義權(quán)函數(shù)為

上式中：sign（·）是符號(hào)函數(shù).因此，式（2）可寫為

用矩陣形式描述，有

式（4）中：Λn＝diag（λn－1，λn－2，…，λn－i，…，1）；Wn＝diag（w1，w2，…，wi，…，wn），wi＝

定義

由式（4）可得

式（5），（6）是新算法的基礎(chǔ)公式.

為確定d（xi）和ξ（xi）的具體形式，定義3個(gè)矩陣，有

上式中：（Wn）1／2＝一般地，這3個(gè)矩陣具有如下4個(gè)性質(zhì)：1）是對(duì)稱，但非冪等的；2）是冪等，但非對(duì)稱的；3）HWn是冪等對(duì)稱的，且其對(duì)角元素hwi的范圍為的對(duì)角元素的范圍為0≤≤1.

引理［11］設(shè)XTnXn可逆，則矩陣Hn＝Xn（XTnXn）－1XTn為冪等對(duì)稱矩陣，且其對(duì)角元素hi的范圍為0≤hi≤1.矩陣Hn的對(duì)角元素通常被用來(lái)檢測(cè)杠桿點(diǎn).

證明 性質(zhì)1），2）的結(jié)論很容易得到，所以只需證明性質(zhì)3）即可.

假設(shè)ˉWn是Wn中非零的相應(yīng)部分，，分別是Λn，Xn中相應(yīng)于的部分，且類似于普通最小二乘法中XTnXn是可逆的假設(shè).新估計(jì)器中假定是可逆的，式（5）可簡(jiǎn)化為WPn＝.而，的乘積是非奇異的，令在HWn中相應(yīng)于ˉWn的部分可以表示為

根據(jù)性質(zhì)1），HWn是冪等對(duì)稱的，且其對(duì)角元素的范圍為0≤≤1.又因?yàn)镠Wn中相應(yīng)于權(quán)重為0的對(duì)角元素也是0，所以對(duì)于HWn的所有對(duì)角元素，有0≤hwj≤1.

下面分析基于式（5），（6）的殘差特性.

定理 殘差rn，i的方差（σ（i）λ）2與之間滿足

證明 根據(jù)式（6），殘差向量的形式有

式（7）中：In為n階單位矩陣.

在觀測(cè)噪聲為獨(dú)立同分布的假設(shè)之下，殘差向量rn的協(xié)方差矩陣為

于是殘差rn，i的方差為

式（9）中：是的第i行第j列的元素；表示其對(duì)角線上的相應(yīng)元素.注意到（λn－j）≤λn－j，j＝1，2，…，n，且0≤wj≤1，易得

因而，把式（11）代入式（10），可得

由性質(zhì)4），又因?yàn)槭剑?2）右邊是非負(fù)的，考慮到在權(quán)重矩陣和遺忘因子確定的情況下是關(guān)于xi的函數(shù)，取并定義一個(gè)統(tǒng)計(jì)量考慮到Hn的對(duì)角元素hi是的特殊形式，且對(duì)普通的最小二乘法而言，hi有檢測(cè)回歸變量中異常點(diǎn)的功能.為了避免的值接近1時(shí)，造成計(jì)算上的誤差，定義ξ（·）函數(shù)為

上式中：K為可調(diào)的參數(shù).該參數(shù)確定數(shù)據(jù)xi是否為異常點(diǎn).如果＜K，認(rèn)為xi是正常的；否則，認(rèn)為xi是異常的.為方便，將ξ（xi），d（xi）分別記作

新遞推估計(jì)器采用一步迭代法推導(dǎo).假定在n－1時(shí)刻的權(quán)矩陣Wn－1和相關(guān)的估計(jì)，Pn－1已經(jīng)得到，那么相應(yīng)的式（5），（6）可以表示為

我八歲的時(shí)候，我父親就去世了，我母親一個(gè)人帶大我們哥倆。我們?cè)趦?nèi)蒙古偏遠(yuǎn)的地方帶大，我在北京沒有一個(gè)親戚，我沒有因?yàn)樽约旱墓ぷ魉瓦^一回禮，我不也走到了今天嗎？我知道社會(huì)上有很多不良的現(xiàn)象，我告訴你，信那些該信的東西，因?yàn)樗芨淖兡?。因?yàn)槿绻阋拍切┠銢]法不憤怒的事情，它只能害了你。

此時(shí)，初始權(quán)重矩陣取W（0）n＝diag（diag（Wn－1），1），而，Pn的初始估計(jì)為

式（15），（16）中：P（－0）n為P（0）n的逆矩陣.

令A(yù)＝，C＝Ip，B＝DT＝xTn，并將矩陣的逆的公式（A＋BCD）－1＝A－1－A－1B（DA－1B＋C－1）－1DA－1應(yīng)用于式（15），則的遞推計(jì)算公式變?yōu)?/p>

迭代開始后，需要計(jì)算再加權(quán)矩陣.此時(shí)，只對(duì)wn進(jìn)行更新，保持W（0）n中的相應(yīng)元素不變，即?。絛iag（diag（Wn－1），wn）.由于λ0＝1，＝1，易知因此，.一個(gè)更好的做法是將W（0）n的對(duì)角元素都進(jìn)行更新.然而，為了得到遞推的計(jì)算公式，這種做法不得不放棄.但是，引入的遺忘因子和所用M－估計(jì)函數(shù)的有界性可以降低這種影響.

式（16）～（19）形成了遞推有界廣義極大似然類（RBGM）回歸估計(jì)器.為了降低計(jì)算量，新估計(jì)器還可以描述為

輸入：，Pn－1，（xn，yn）；

輸出：，Pn；

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

在觀測(cè)數(shù)據(jù)中含有較大量加性異常點(diǎn)的情況下進(jìn)行AR 模型參數(shù)估計(jì)時(shí)，RBGM 也是有偏的.為此，需要對(duì)其進(jìn)行必要的改造，把面向AR 模型參數(shù)估計(jì)的算法記作AR－RBGM.

首先，用，替換RBGM 中的xn，βn.＝［b1，b2，…，bp，η1］，＝［yn－1，yn－2，…，yn－p，εn－1］（η1為增廣參數(shù)，εn－1為增廣變量）.在初始化時(shí)，εn－1和η1 都賦0值.在運(yùn)行期間，隨著算法自動(dòng)更新.

其次，增加一個(gè)加性異常點(diǎn)的檢測(cè)過程.在此過程中涉及兩種殘差的計(jì)算，一是不包括增廣項(xiàng)的rn，n，二是包括增廣項(xiàng)的殘差＝y(tǒng)n－.在當(dāng)前n時(shí)刻的估計(jì)已經(jīng)得到后，假設(shè)以前的加性異常點(diǎn)的影響已經(jīng)體現(xiàn)在增廣項(xiàng)里，且當(dāng)前的觀測(cè)數(shù)據(jù)不是加性異常點(diǎn)，那么相應(yīng)的接受域?yàn)椋鄳?yīng)的拒絕域?yàn)槠渲校害檬莻€(gè)可調(diào)常數(shù)，一般取1.5左右即可.據(jù)檢驗(yàn)結(jié)果可知：如果拒絕了假設(shè)，則認(rèn)為yn是加性異常點(diǎn)，其影響應(yīng)該體現(xiàn)在增廣項(xiàng)中，相應(yīng)地取εn＝rn，n；否則，取εn＝0.此外，是用替換rn，j后，通過式（20）計(jì)算得到的，而在用RLS進(jìn)行啟動(dòng)的階段，可直接用med估計(jì)，也可以借用先驗(yàn)估計(jì)結(jié)果.

在AR 模型參數(shù)估計(jì)時(shí)，無(wú)論數(shù)據(jù)中包含或不包含異常點(diǎn)，異常點(diǎn)是革新或加性，AR－RBGM 中的參數(shù)都無(wú)須重置.另外，RBGM 和AR－RBGM 追蹤系統(tǒng)突變的能力主要依賴于遺忘因子的大小，當(dāng)需要追蹤突變時(shí)，遺忘因子相對(duì)取小一點(diǎn)，如0.99等；不需要追蹤突變時(shí)，遺忘因子則取得相對(duì)大一點(diǎn)，如0.999或0.999 9等.在估計(jì)精度和對(duì)突變的追蹤能力方面需要折中取舍.

4 仿真實(shí)驗(yàn)

在實(shí)驗(yàn)中，假定需要辨識(shí)的真實(shí)系統(tǒng)為

式（20）中：ei為相互獨(dú)立同分布的隨機(jī)變量，且ei～N（0，σ21）.該模型已經(jīng)在文獻(xiàn)［5，10］中被用來(lái)檢驗(yàn)M 和GM 估計(jì)器的魯棒性.

在觀測(cè)值的仿真實(shí)驗(yàn)中，如果第一型異常點(diǎn)出現(xiàn)，相應(yīng)的觀測(cè)值yi用yi＝zi＋vi進(jìn)行模擬.其中：vi～（1－κ1）ד0”＋κ1×N（0，σ22），如果第二型異常點(diǎn)出現(xiàn)，相應(yīng)的觀測(cè)值yi通過把式（21）中的ei替換為另一個(gè)隨機(jī)變量ni，并令yi＝zi進(jìn)行模擬.其中：ni～N（0，σ21）＋κ2×N（0，σ23）.在仿真實(shí)驗(yàn)中，設(shè)計(jì)了4種更加復(fù)雜的情形，觀測(cè)數(shù)據(jù)如表1所示.每一種情形都將整個(gè)過程分成3個(gè)階段，并且假定第一型異常點(diǎn)或者第二型異常點(diǎn)分別出現(xiàn)在某一個(gè)階段.

在所有情形中，取σ21＝25，κ1＝0.05，κ2＝0.05，σ22＝σ23＝400.采用對(duì)數(shù)化平均相對(duì)誤差（LMRE）來(lái)刻畫估計(jì)精度及收斂性，即

式（21）中：M為程序運(yùn)行的次數(shù)；β為系統(tǒng)的真實(shí)參數(shù)；為在第m次運(yùn)行時(shí)，n時(shí)刻β的估計(jì).

表1 4種情形的觀測(cè)數(shù)據(jù)Tab.1 Four types of observed data

實(shí)驗(yàn)中用經(jīng)典的RLS和另外兩種遞推的GM 估計(jì)器RGM［8］和RKW［4］做了比較.在沒有異常點(diǎn)的情況下這幾種估計(jì)器的表現(xiàn)，如圖1所示.圖1 中：LMRE 為對(duì)數(shù)化平均相對(duì)誤差；n為時(shí)序.在比較時(shí)，每種估計(jì)器都進(jìn)行了參數(shù)調(diào)整以使其表現(xiàn)盡可能達(dá)到最好.具體地，對(duì)RLS，λ＝0.999 9；對(duì)RGM，λ＝0.999 9，c＝2.8，P（0）＝100I2；對(duì)RKW，P（0）＝A－1（0）＝100I2，λ＝0.999 9，c＝2，a＝5；對(duì)RBGM 和AR－RBGM，λ＝0.999 9，NInitial＝10，λσ＝0.98，γ＝1.5，L＝20，K＝0.85，M1＝1.88，M2＝2.41.為了方便比較，所有的初始尺度估計(jì)中，σ0＝5.實(shí)驗(yàn)結(jié)果表明：在此情形下，除了AR－RBGM 的表現(xiàn)稍差以外，其他估計(jì)器表現(xiàn)相當(dāng).因?yàn)锳R－RBGM 中的增廣項(xiàng)改變了原始模型的結(jié)構(gòu)，對(duì)其性能產(chǎn)生了有限的影響.

非平穩(wěn)環(huán)境（情形2～4）下估計(jì)性能比較，如圖2～4所示.圖2～4中：LMRE為對(duì)數(shù)化平均相對(duì)誤差；n為時(shí)序.由圖2，3可知：即使當(dāng)加性異常點(diǎn)出現(xiàn)在算法已達(dá)到或者接近穩(wěn)定狀態(tài)時(shí)，RLS，RGM，RKW 完全失效，而RBGM 對(duì)這些影響抑制效果是明顯的，但精度仍然不高；AR－RBGM 能在各種情形下都保持較高的精度和收斂性，其精度與沒有異常點(diǎn)時(shí)的精度差別很小.由圖4可知：對(duì)加性異常點(diǎn)出現(xiàn)在初始階段的情形，仍然是AR－RBGM 保持較高精度和良好收斂性，其他估計(jì)器全部失；加性異常點(diǎn)出現(xiàn)在起始階段時(shí)，AR－RBGM 的相應(yīng)的指標(biāo)會(huì)變大一點(diǎn)，這主要是受到非魯棒啟動(dòng)算法RLS的影響.

圖1 平穩(wěn)環(huán)境（情形1）下估計(jì)性能比較Fig.1 Performance comparison of the algorithm for the stability environment（case 1）

圖2 非平穩(wěn)環(huán)境（情形2）下估計(jì)性能比較Fig.2 Performance comparison of the algorithm for the unstability environment（case 2）

圖3 非平穩(wěn)環(huán)境（情形3）下估計(jì)性能比較Fig.3 Performance comparison of the algorithm for the unstability environment（case 3）

圖4 非平穩(wěn)環(huán)境（情形4）下估計(jì)性能比較Fig.4 Performance comparison of the algorithm for the unstability environment（case 4）

5 結(jié)束語(yǔ)

改進(jìn)GM 估計(jì)器的一般框架，提出一種新的遞推GM 回歸估計(jì)器（RBGM），并針對(duì)AR 模型參數(shù)估計(jì)提出了AR－RBGM.RBGM 和AR－RBGM 均能對(duì)回歸變量中的加性異常點(diǎn)的影響起到抑制作用，特別是AR－RBGM 能在非平穩(wěn)環(huán)境下實(shí)現(xiàn)自適應(yīng)的估計(jì)，并保持良好的精度和收斂性.新估計(jì)器還可以進(jìn)行改進(jìn)，一方面，魯棒的啟動(dòng)算法可以提高性能；另一方面，可進(jìn)一步推廣到ARMA 模型參數(shù)的估計(jì)中，從而獲得相應(yīng)的魯棒算法.然而，由于AR－RBGM 引入了增廣變量，增加了額外的計(jì)算量.

［1］HUBER P J.Robust regression：Asymptotics，conjectures and monte carlo［J］.Annals of Statistics，1973，1（5）：799－821.

［2］CAMPBEL K.Recursive computation of M－estimates for the parameters of a finite autoregressive process［J］.The Annals of Stat，1982，10（2）：442－453.

［3］ANTOCH J，EKBLOM H.Recursive robust regression computational aspects and comparison［J］.Computational Statistics and Data Analysis，1995，19（2）：115－128.

［4］SEJLING K，et al.Methods for recursive robust estimation of AR parameters［J］.Computational Atatistics and Data Analysis，1994，17（5）：509－536.

［5］PHAM D S，ZOUBIR A M.A sequential algorithm for robust parameter estimation［J］.IEEE Signal Processing Lett，2005，12（1）：21－24.

［6］VEGA L R，REY H，BENESTY J，et al.A robust recursive least squares algorithm［J］.IEEE Trans Signal Process，2009，57（3）：1209－1216.

［7］KRASKER W S，WELSCH R E.Efficient bounded－influence regression estimation［J］.Journal of the American Statistical Association，1982，77（379）：595－604.

［8］GRILLENZONI C.Recursive generalized M－estimators of system parameters［J］.Technometrics，1997，39（2）：211－224.

［9］ENGIUND J E.Recursive versions of the algorithm by Krasker and Welsch［J］.Sequential Analysis，1991，10（3／4）：211－234.

［10］MARONNA R A，MARTIN R D，YOHAI V J.Robust statistics：Theory and methods［M］.West Sussex：John Wiley＆Sons，2006：888－889.

［11］ROUSSEEUW P J，LEROY A M.Robust regression and outlier detection［M］.New York：Wiley，1987：340－347.