新的递推有界GM 回归估计算法

成立花，张俊敏

（1.西安工程大学 理学院，陕西 西安710048；2.西安建筑科技大学 理学院，陕西 西安710055）

最小二乘法（LS）和递推的最小二乘法（RLS）是系统辨识中重要的算法.然而，在实际应用中，由于各种干扰因素的存在，辨识所需的数据往往会被污染，包含数量和类型均未知的异常点.这些异常点使估计精度降低，甚至完全失效.解决这一问题通常有两种方式：从数据中找出异常点并剔除掉，用剩余的数据进行辨识；发展替代LS和RLS的鲁棒回归算法.由于后者简便实用，备受学者关注，已发展出了许多算法［1－9］.在自回归模型的参数估计中，数据中的异常点通常被分为两个基本类［10］：加性异常点（第一型的异常点）与革新异常点（第二型的异常点）.革新异常点在系统运行过程中遵循了真实系统的关系，因此，就系统的辨识而言，其不利影响较小.然而，当数据中包含较大量的加性异常点型杠杆点时，这些算法的估计性能会严重下降，甚至失效.针对广义极大似然类（GM）估计器中存在的问题，本文提出新的递推有界GM 回归估计算法.

1 问题描述

考虑模型yi＝xiβ＋vi，i＝1，2，…，n.其中：β＝［b1，b2，…，bp］T；yi是第i时刻响应变量观测值；x＝［xi，1，xi，2，…，xi，p］；vi是相互独立，且尺度参数σ未知的干扰项.所用的数据记作（Xn，Yn）.其中：Xn＝［xT2，xT2，…，xTn］T，Y＝［y1，y2，…，yn］T.获得第n时刻的参数估计后，第i时刻相应的残差记作rn，i＝yi－.而在迭代法中，相应于起始估计在第i时刻的残差记作

采用文献［9］中的一般框架，并加以改进的风险函数为

式（1）中：λ为遗忘因子，0＜λ≤1；d（xi），ξ（xi）为待定函数；ρ（·）是一个改进于Huber函数的有界M－估计函数，即

上式中：M为常值参数.新估计器是式（1）的解.

2 新的GM 估计器

相应于式（1）的“正规方程”为

式（2）中：ψ＝ρ′.

定义权函数为

上式中：sign（·）是符号函数.因此，式（2）可写为

用矩阵形式描述，有

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

定义

由式（4）可得

式（5），（6）是新算法的基础公式.

为确定d（xi）和ξ（xi）的具体形式，定义3个矩阵，有

上式中：（Wn）1／2＝一般地，这3个矩阵具有如下4个性质：1）是对称，但非幂等的；2）是幂等，但非对称的；3）HWn是幂等对称的，且其对角元素hwi的范围为的对角元素的范围为0≤≤1.

引理［11］设XTnXn可逆，则矩阵Hn＝Xn（XTnXn）－1XTn为幂等对称矩阵，且其对角元素hi的范围为0≤hi≤1.矩阵Hn的对角元素通常被用来检测杠杆点.

证明 性质1），2）的结论很容易得到，所以只需证明性质3）即可.

假设ˉWn是Wn中非零的相应部分，，分别是Λn，Xn中相应于的部分，且类似于普通最小二乘法中XTnXn是可逆的假设.新估计器中假定是可逆的，式（5）可简化为WPn＝.而，的乘积是非奇异的，令在HWn中相应于ˉWn的部分可以表示为

根据性质1），HWn是幂等对称的，且其对角元素的范围为0≤≤1.又因为HWn中相应于权重为0的对角元素也是0，所以对于HWn的所有对角元素，有0≤hwj≤1.

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

定理 残差rn，i的方差（σ（i）λ）2与之间满足

证明 根据式（6），残差向量的形式有

式（7）中：In为n阶单位矩阵.

在观测噪声为独立同分布的假设之下，残差向量rn的协方差矩阵为

于是残差rn，i的方差为

式（9）中：是的第i行第j列的元素；表示其对角线上的相应元素.注意到（λn－j）≤λn－j，j＝1，2，…，n，且0≤wj≤1，易得

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

由性质4），又因为式（12）右边是非负的，考虑到在权重矩阵和遗忘因子确定的情况下是关于xi的函数，取并定义一个统计量考虑到Hn的对角元素hi是的特殊形式，且对普通的最小二乘法而言，hi有检测回归变量中异常点的功能.为了避免的值接近1时，造成计算上的误差，定义ξ（·）函数为

上式中：K为可调的参数.该参数确定数据xi是否为异常点.如果＜K，认为xi是正常的；否则，认为xi是异常的.为方便，将ξ（xi），d（xi）分别记作

新递推估计器采用一步迭代法推导.假定在n－1时刻的权矩阵Wn－1和相关的估计，Pn－1已经得到，那么相应的式（5），（6）可以表示为

我八岁的时候，我父亲就去世了，我母亲一个人带大我们哥俩。我们在内蒙古偏远的地方带大，我在北京没有一个亲戚，我没有因为自己的工作送过一回礼，我不也走到了今天吗？我知道社会上有很多不良的现象，我告诉你，信那些该信的东西，因为它能改变你。因为如果你要信那些你没法不愤怒的事情，它只能害了你。

此时，初始权重矩阵取W（0）n＝diag（diag（Wn－1），1），而，Pn的初始估计为

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

令A＝，C＝Ip，B＝DT＝xTn，并将矩阵的逆的公式（A＋BCD）－1＝A－1－A－1B（DA－1B＋C－1）－1DA－1应用于式（15），则的递推计算公式变为

迭代开始后，需要计算再加权矩阵.此时，只对wn进行更新，保持W（0）n中的相应元素不变，即取＝diag（diag（Wn－1），wn）.由于λ0＝1，＝1，易知因此，.一个更好的做法是将W（0）n的对角元素都进行更新.然而，为了得到递推的计算公式，这种做法不得不放弃.但是，引入的遗忘因子和所用M－估计函数的有界性可以降低这种影响.

式（16）～（19）形成了递推有界广义极大似然类（RBGM）回归估计器.为了降低计算量，新估计器还可以描述为

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

输出：，Pn；

3 AR 模型参数估计算法

在观测数据中含有较大量加性异常点的情况下进行AR 模型参数估计时，RBGM 也是有偏的.为此，需要对其进行必要的改造，把面向AR 模型参数估计的算法记作AR－RBGM.

首先，用，替换RBGM 中的xn，βn.＝［b1，b2，…，bp，η1］，＝［yn－1，yn－2，…，yn－p，εn－1］（η1为增广参数，εn－1为增广变量）.在初始化时，εn－1和η1 都赋0值.在运行期间，随着算法自动更新.

其次，增加一个加性异常点的检测过程.在此过程中涉及两种残差的计算，一是不包括增广项的rn，n，二是包括增广项的残差＝yn－.在当前n时刻的估计已经得到后，假设以前的加性异常点的影响已经体现在增广项里，且当前的观测数据不是加性异常点，那么相应的接受域为，而相应的拒绝域为其中：γ是个可调常数，一般取1.5左右即可.据检验结果可知：如果拒绝了假设，则认为yn是加性异常点，其影响应该体现在增广项中，相应地取εn＝rn，n；否则，取εn＝0.此外，是用替换rn，j后，通过式（20）计算得到的，而在用RLS进行启动的阶段，可直接用med估计，也可以借用先验估计结果.

在AR 模型参数估计时，无论数据中包含或不包含异常点，异常点是革新或加性，AR－RBGM 中的参数都无须重置.另外，RBGM 和AR－RBGM 追踪系统突变的能力主要依赖于遗忘因子的大小，当需要追踪突变时，遗忘因子相对取小一点，如0.99等；不需要追踪突变时，遗忘因子则取得相对大一点，如0.999或0.999 9等.在估计精度和对突变的追踪能力方面需要折中取舍.

4 仿真实验

在实验中，假定需要辨识的真实系统为

式（20）中：ei为相互独立同分布的随机变量，且ei～N（0，σ21）.该模型已经在文献［5，10］中被用来检验M 和GM 估计器的鲁棒性.

在观测值的仿真实验中，如果第一型异常点出现，相应的观测值yi用yi＝zi＋vi进行模拟.其中：vi～（1－κ1）ד0”＋κ1×N（0，σ22），如果第二型异常点出现，相应的观测值yi通过把式（21）中的ei替换为另一个随机变量ni，并令yi＝zi进行模拟.其中：ni～N（0，σ21）＋κ2×N（0，σ23）.在仿真实验中，设计了4种更加复杂的情形，观测数据如表1所示.每一种情形都将整个过程分成3个阶段，并且假定第一型异常点或者第二型异常点分别出现在某一个阶段.

在所有情形中，取σ21＝25，κ1＝0.05，κ2＝0.05，σ22＝σ23＝400.采用对数化平均相对误差（LMRE）来刻画估计精度及收敛性，即

式（21）中：M为程序运行的次数；β为系统的真实参数；为在第m次运行时，n时刻β的估计.

表1 4种情形的观测数据Tab.1 Four types of observed data

实验中用经典的RLS和另外两种递推的GM 估计器RGM［8］和RKW［4］做了比较.在没有异常点的情况下这几种估计器的表现，如图1所示.图1 中：LMRE 为对数化平均相对误差；n为时序.在比较时，每种估计器都进行了参数调整以使其表现尽可能达到最好.具体地，对RLS，λ＝0.999 9；对RGM，λ＝0.999 9，c＝2.8，P（0）＝100I2；对RKW，P（0）＝A－1（0）＝100I2，λ＝0.999 9，c＝2，a＝5；对RBGM 和AR－RBGM，λ＝0.999 9，NInitial＝10，λσ＝0.98，γ＝1.5，L＝20，K＝0.85，M1＝1.88，M2＝2.41.为了方便比较，所有的初始尺度估计中，σ0＝5.实验结果表明：在此情形下，除了AR－RBGM 的表现稍差以外，其他估计器表现相当.因为AR－RBGM 中的增广项改变了原始模型的结构，对其性能产生了有限的影响.

非平稳环境（情形2～4）下估计性能比较，如图2～4所示.图2～4中：LMRE为对数化平均相对误差；n为时序.由图2，3可知：即使当加性异常点出现在算法已达到或者接近稳定状态时，RLS，RGM，RKW 完全失效，而RBGM 对这些影响抑制效果是明显的，但精度仍然不高；AR－RBGM 能在各种情形下都保持较高的精度和收敛性，其精度与没有异常点时的精度差别很小.由图4可知：对加性异常点出现在初始阶段的情形，仍然是AR－RBGM 保持较高精度和良好收敛性，其他估计器全部失；加性异常点出现在起始阶段时，AR－RBGM 的相应的指标会变大一点，这主要是受到非鲁棒启动算法RLS的影响.

图1 平稳环境（情形1）下估计性能比较Fig.1 Performance comparison of the algorithm for the stability environment（case 1）

图2 非平稳环境（情形2）下估计性能比较Fig.2 Performance comparison of the algorithm for the unstability environment（case 2）

图3 非平稳环境（情形3）下估计性能比较Fig.3 Performance comparison of the algorithm for the unstability environment（case 3）

图4 非平稳环境（情形4）下估计性能比较Fig.4 Performance comparison of the algorithm for the unstability environment（case 4）

5 结束语

改进GM 估计器的一般框架，提出一种新的递推GM 回归估计器（RBGM），并针对AR 模型参数估计提出了AR－RBGM.RBGM 和AR－RBGM 均能对回归变量中的加性异常点的影响起到抑制作用，特别是AR－RBGM 能在非平稳环境下实现自适应的估计，并保持良好的精度和收敛性.新估计器还可以进行改进，一方面，鲁棒的启动算法可以提高性能；另一方面，可进一步推广到ARMA 模型参数的估计中，从而获得相应的鲁棒算法.然而，由于AR－RBGM 引入了增广变量，增加了额外的计算量.

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