第七章-稳健回归---浙江大学数学科学学院

第八章稳健回归§8.1异常值例8.1设)1,0(~N，利用随机数产生6个独立观测值，,,,621对于给定的ix值，按如下模型产生y的观测值：6,1226,2ixyixyiiiiii显然第6个点为异常值。我们给出如下的SAS模拟计算的程序：libnametj’d:user’;datatj.a1;doi=1to10;settj.p183;r=normal(1234567);ifi=6theny=-2-x+12+r;elsey=-2-x+r;output;end;run;procregdata=tj.a1;modely=x;outputout=b1p=y1hat;run;datatj.a2;settj.a1;ifx5thendelete;run;procregdata=tj.a2;modely=x;outputout=b2p=y2hat;run;procsortdata=b1;byx;run;procsortdata=b2;byx;run;datab;mergeb1b2;run;procgplotdata=b;symbol1v=starc=greenI=none;symbol2v=plusc=redi=l;symbol3v=plusc=bluei=l;ploty*xhaty1*xhaty2*x/overlay;run;第一个方程为：xy32.095.0ˆ第二个方程为：xy01.14127.1ˆ§7.2M估计在最小二乘估计中，回归系数ˆ是使niie12达到最小，其中tjijjiixye0，若令2)(~xx，则ˆ使niiniiee121)(~达到最小。显然函数2)(~xx随x的增大而迅速增加，为使上式达到最小，就必须照顾某些点，特别是一些异常点如上例中的第6个点。从上面的图可以看出，第6点去掉后，回归直线的位置发生了较大的变化。为了减小异常点对估计的影响，就要设法找到一个函数)(x去替代2x，该函数除了随x的增长速度应比2x慢外，其它性质与2x类似。这便是M估计的想法。定义7.2.1设函数)(x满足下列条件：（1）定义在,上的连续函数；（2）)()(,0)0(xx，当0yx时有)()(yx；（3）0)(lim2xxx记nixyemjijjii,,2,1,1，又令)(,,,)(21ime若*满足)(min)(*则称*是的M估计。当)(x关于的导数)()(xx存在时，则*应满足1,,,2,1,0)(1tmmjxeniiji在一定条件下，方程组必有解，且其解即为M估计。定理7.2.2设在上面定义中，)(x是定义在,上的一个凸函数，)()(xx处处存在，且)(limxx，mXR)(，则方程组1,,,2,1,0)(1tmmjxeniiji必有解满足)(min)(*。反之，任何满足条件)(min)(*的*必为1,,,2,1,0)(1tmmjxeniiji的解。若进一步假定)(x为严格凸函数，则解为唯一。在上面定中主要涉及了)(x的导数，所以也有人直接从这个角度去定义M估计。定义7.2.3设函数)(x满足下列条件：（1）定义在,上的连续函数；（2）(0)0,()()xx，0x时有0,0)(xx时0)(x；（3）0)(limxxx则方程组1,,,2,1,0)(1tmmjxeniiji的解*是的M估计。常用的)(x有：（1）Huber提出的,,0,,,,)(1kxkkkxxkxkx（2）Hampel提出的cxcxbbcxcxsignabxaxsignaaxxsignxx,0,)()()()(2其中cba0，.0,1,0,0,0,1)(xxxxsign(3)Andrews提出的,,0,),sin()(2cxcxcxcxx由于上述函数所决定的方程组是一个非线性方程组，因而常常只能用迭代法来求解，为此需要给出初值、迭代公式及收敛准则。（1）初值的确定为了能与最小二乘估计作比较，因而在M估计中常将最小二乘估计作为初值，即取初值ˆ)0(。（2）迭代公式设)(l是第l次迭代值，其迭代公式为：,,2,1,0,)()()1(llll这里)(l是增量。为求)(l常假定)(x在)(l处可作台劳展开，并取下列近似式：nixeeelkikmklilii,,2,1,0,)()()()(1其中nixyeijmjljili,,2,1,1)(将展开式代入方程组，可得：mjxxeexeniijnilkikmililiiji,,2,1,0)()()(11)(1上述方程组是关于)()(2)(1,,,lmll的线性方程组。因而很容易解出)(l（3）收敛准则这类的收敛准则很多，如：事先给定一个0，要求)(maxljj；事先给定一个0，要求mjlj1:2)(；§7.2R估计R估计是将原始数据转换成相应的秩后，利用最小二乘建立y的秩关于1x的秩、2x的秩、…tx的秩间的回归方程，先求出y的秩的预测值后再求出y的预测值的一种方法。它是一种非参数估计方法，又由于秩受异常值的影响小，所以它也是一种稳健回归的方法。R估计的具体步骤如下：（1）设给出的原始数据为：niyxxxitii,,2,1,,,,,21记：ijx在itiixxx,,,21中的秩为tjnixRij,,2,1,,,2,1,，记iy在nyyy,,,21中的秩为iyR，则对应的秩的数据记为：niyRxRxRxRitii,,2,1,,,,,21(2)利用最小二乘估计建立yR关于,,,,21itiixRxRxR的回归方程：ittiixRxRxRyRˆˆˆˆˆ22110(3)在给定点txxx00201,,,作预测：先将jx0转换成为jxR0，方法如下：设itiixxx,,,21排成次序统计量为jnjjxxx)()2()1(，则：00110000111,,,,,1,2,1,1,2jnjnjjjjjjijijjijjijijijijijijijRxxxRxxxRxRxxxxxRxRxRxxxxxxinjt将txRxRxR00201,,,代入秩回归方程，求出0ˆyR.由0ˆyR求出0ˆy，方法如下：设nyyy21，则1,,2,1,1)(ˆ,)(ˆ,)(ˆ,,)(ˆ,,1)(ˆ,ˆ00100010niiyRiiyRyyyiyRynyRyyRyyiiiin例在把氨氧化成硝酸的生产过程中，希望求得氨的损失率y关于空气流速x1，冷却水温度x2，吸收液中硝酸的浓度x3的回归方程。收集到21天的数据。经分析其中21天的用最小二乘法估计后所得残差特别大，并根据生产经验知第1，3，4天的生产条件也不正常。下面是用不同的估计方法所得的21天，及去掉不正常的天后的数据，即17天的数据所作回归方程后的残差平方和最小二乘法M估计R估计21点178.83245.03109.1717点（除去1，3，4，和21点）16.4420.6019.21