回归分析第16讲——稳健回归模型

通过前面的学习，我们已经了解了各种经典回归模型的架构。然而，无论哪种回归模型，都存在一个问题：对数据集中的异常值过于敏感。如果数据集中出现了异常值，就会使得回归平面向异常值偏移，造成拟合偏差。你或许想把这些异常值剔除掉，但是在实际应用中，有些异常值是真实观测值而非仪器、人为等原因造成，这些点不能简单剔除，否则会影响模型的泛化性。因此，我们需要在不改变模型结构的情况下使得模型能够“容忍”异常值，削弱异常值对回归平面的影响。这就是本讲要介绍的稳健回归(Robust Regression)模型。

泛化性(Generalization)，是指已完成训练的模型在遇到没有见过的新数据时也能够表现出良好的能力，这是所有模型都希望做到的一点。如果模型对异常值的容忍度很差，那么在遇到新数据时表现就会很差，即出现很大的预测误差。

本讲中，我们先补充稳健统计的预备知识，之后再介绍一些常见的稳健回归模型。

一、稳健统计简介

1.1 稳健统计量的概念

稳健(Robustness)，是指统计量在遇到异常值时仍然能够保持精度(Precision)和准度(Accuracy)，而不会因为异常值的出现使得统计量发生大幅变化。例如，我们熟知的样本均值和方差就是典型的非稳健的统计量：

$$\overline x=\frac{1}{n}\sum_{i=1}^nx_i\ \ \ s^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\overline x)^2$$

如果在样本中出现了一个极端值，会显著影响两个统计量的值，从而使得我们无法把握绝大部分“正常值”的特征，因此我们说这两个统计量是非稳健的。那么什么样的统计量是稳健的呢？

例如，样本的中位数就是稳健的。如果我们称一个数$t$是一组样本$x=(x_1,\cdots,x_n)$的中位数，那么$t$可以记作：

$$t=Med(x),if\ \ \#(x_i>t)=\#(x_i<t)$$

其中$\#()$代表满足括号内条件的样本的个数。也就是说，中位数的左右两边的样本量是一致的。

如果我们把样本重新编为顺序统计量$x=(x_{(1)},\cdots,x_{(n)})$，那么有：

$$Med(x)= \begin{cases} x_{(m)},&n\ is\ odd,\ n=2m-1\\ (x_{(m)}+x_{(m+1)})/2,&n\ is\ even,\ n=2m \end{cases}$$

可以看出，样本中位数是样本均值的稳健版本，因为如果出现了极端值，中位数只会发生微弱变化（中位数只与相邻的数有关）。

理解了什么是稳健统计量后，接下来我们看一些常见的稳健统计量。一般地，我们把稳健统计量分为两类：一类是描述数据位置的统计量，另一类是描述数据波动的统计量。

1.2 数据位置的稳健统计量

数据位置，是指一组样本数据的分布中心(Location)。严谨来说，我们认为每一个样本的观测值$x_i$由两部分组成：一部分是测量对象的真实值$\mu$，另一部分是随机误差$u_i$，即：

$$x_i=\mu+u_i\ \ \ i=1,\cdots,n$$

如果误差独立同分布$F_0$，那么观测值也独立同分布：

$$F(x)=F_0(x-\mu)$$

数据位置的稳健统计量，就是要找到一个稳健估计子$\hat\mu$，使得$\hat\mu\approx\mu$。一个好的位置稳健估计子应当具有两个性质，第一是平移不变性(Shift Equivariance)，即：

$$\hat\mu(x_1+c,\cdots,x_n+c)=\hat\mu(x_1,\cdots,x_n)+c$$

第二是尺度不变性(Scale Equivariance)，即：

$$\hat\mu(cx_1,\cdots,cx_n)=c\hat\mu(x_1,\cdots,x_n)$$

很多时候，我们会假设随机误差精确服从均值为0的正态分布，那么观测值就服从均值为$\mu$的正态分布，此时样本均值就是最优的估计子，并且具有平移不变性和尺度不变性。然而，实际测量数据是渐进正态的，所以为了构造稳健统计量，我们需要对观测值的分布进行新的假定。一般有两种方式。

第一种称为污染正态分布(Contaminated Normal Distribution)，又称为Tukey-Huber模型，它假设有$\epsilon$比例的观测值不是由正态分布G产生的，而是由其他未知分布H产生，因此观测值的分布可以写成分布的混合叠加，即：

$$F=(1-\epsilon)G+\epsilon H$$

第二种是使用重尾分布(Heavy-tailed Distribution)，即概率密度趋近于0的速度更慢，例如$Student$分布家族：

$$f(x,v)=\frac{\Gamma((v+1)/2)}{\sqrt{v\pi}\ \Gamma(v/2)}\left(1+\frac{x^2}{v}\right)^{-(v+1)/2}$$

常见可取$v=1$，这样就变成了柯西分布(Cauchy Distribution)。

有了估计模型和分布之后，我们来考虑数据位置的稳健估计量。数据位置的稳健估计量可以分为两大类：M估计和L估计。

M估计。M估计是Huber在1964年提出的一类估计，其基本思想是，最小化测量值和真实值之间的随机误差，即：

$$\hat\mu=\arg\min_\mu\sum_{i=1}^n\rho(x_i-\mu)$$

其中$\rho=-\log f_0$，称为$\rho$函数；$f_0$是随机误差$u_i$的概率密度函数。

进一步地，如果$f_0$可微分，那么令$\psi=\rho^{'}$，那么上述优化目标相当于求解：

$$\sum_{i=1}^n\psi(x_i-\hat\mu)=0$$

如果随机误差服从标准正态分布（白噪声），那么可以求得$\rho(x)=x^2/2$，于是就类似于最小二乘法，相当于最小化平方损失。通过上述方程可以解得$\hat\mu=\overline x$，此时M估计就是样本均值。

如果随机误差服从双指数分布，即$f_0(x)=e^{- |x|}/2$，那么可以求得$\rho(x)= |x|$。这种估计被称为最小绝对偏差(Least Absolute Deviation, LAD)。可以解得$\hat\mu=Med(x)$，此时M估计就是样本中位数。

当然，除了这些，Huber本人还提出了自己的估计方法，即取：

$$\rho(x,k)=\begin{cases} \frac{1}{2}x^2,|x|\le k\\\\ k|x|-\frac{1}{2}k^2,|x|>k \end{cases}$$

还有一些别的常用方法，例如Hampel方法是取：

$$\psi(x,a,b,c)= \begin{cases} x,|x|\le a\\ a\times sign(x),a< |x|\le b\\ a\times\frac{c- |x| }{c-b}\times sign(x),b< |x|\le c\\ 0,|x| >c \end{cases}$$

Bisquare方法（或称Tukey’s Biweight方法）是取：

$$\psi(x,c)= \begin{cases} x\left(1-\left(\frac{x}{c}\right)^2\right)^2,|x|\le c\\ 0,|x| >c \end{cases}$$

可以证明，所有的M估计子都具有平移不变性；均值和中位数估计子具有放缩不变性，但并不是所有的M估计子都具有放缩不变性。

L估计。L估计又称为线性组合估计，也就是说，真实值的估计值表达为各个样本测量值的线性组合。一般地，我们使用顺序统计量的权重之和来估计，即：

$$\hat\mu=\sum_{i=1}^n a_ix_{(i)}$$

这里我们介绍两种L估计。第一种称为截断均值，即去掉前$m$个最小值和后$m$个最大值之后，再取平均：

$$\overline x_\alpha=\frac{1}{n-2m}\sum_{i=m+1}^{n-m}x_{(i)}$$

其中$\alpha\in[0,1/2)$，$m=[n\alpha]$，$[]$表示取整。

如果写成L估计的标准形式就是取：

$$a_i=\frac{1}{n-2m}I(m+1\le i\le n-m)$$

第二种是Winsor均值，其实就是在截断均值的基础上加入两个截断端点的值，但是依据样本总量来平均，即：

$$\overline x_\alpha=\frac{1}{n}\left(mx_{(m)}+mx_{(n-m+1)}+\sum_{i=m+1}^{n-m}x_{(i)}\right)$$

其中$m,\alpha$同前。

如果写成L估计的标准形式就是取：

$$a_i=\frac{1}{n}\left[mI(i=m)+mI(i=n-m+1)+I(m+1\le i\le n-m)\right]$$

可以证明，当满足下列条件时，L估计子都满足平移不变性和放缩不变性：

$$a_i\ge0,\sum_{i=1}^n a_i=1,a_i=a_{n-i+1}$$

1.3 数据波动的稳健统计量

除了考虑数据分布的中心，还要考虑数据分布的波动性(Dispersion)。描述数据分布波动的方法大致有两类：一类是带尺度参数的形式，另一类是非参数的形式。

带尺度参数的形式。带尺度参数，是指样本观测值可以表示为随机误差的倍数，即：

$$x_i=\sigma u_i\ \ i=1,\cdots,n$$

其中$\sigma$被称为尺度参数。

要估计尺度参数，可以使用1.2小节中提到的M估计：

$$\frac{1}{n}\sum_{i=1}^n\rho\left(\frac{x_i}{\hat\sigma}\right)=\delta$$

其中$\rho()$代表M估计中的各种$\rho$函数，$\delta$是大于0的常数。

要使得上述方程有解，需要规定$0<\delta<\rho(\infty)$。例如，如果$\rho$函数在无穷处的上界为1，那么$\delta\in(0,1)$。

非参数的形式。有时候我们无法假定数据波动的形式，就需要一些不含参数的统计量来描述数据的波动性。我们最开始接触的统计量就是标准差（或方差），即：

$$SD(x)=\left[\frac{1}{n-1}\sum_{i=1}^n(x_i-\overline x)^2\right]^{1/2}$$

不难看出，标准差具有平移不变性和尺度不变性。

然而，标准差是不稳健的，因为它和样本均值紧密联系。因此，我们应该考虑构建一个稳定版本的标准差/方差。我们首先定义中位绝对偏差(median absolute deviation about the median, MAD)：

$$MAD(x)=Med(|x-Med(x)|)$$

MAD就是稳健版本的样本方差，并且具有平移不变性和尺度不变性。

有了MAD就可以考虑构造稳健版本的样本标准差，不过我们肯定不能直接开平方根。我们这样构造：

$$MADN(x)=\frac{MAD(x)}{0.6745}=1.4826\times MAD(x)$$

这个被称为标准化MAD(Normalized MAD, MADN)，相当于是稳健版本的样本标准差。

你或许对0.6745感到很莫名其妙。取这个数是因为一个标准正态分布的MAD大约是其标准差的0.6745倍，因此我们采用这个值来进行标准化。

除了MAD，我们再介绍两个统计量用于描述数据的波动性。

极差(Range)。极差就是样本最大值和最小值的差，即：

$$max(x)-min(x)=x_{(n)}-x_{(1)}$$

极差是不稳健的，因为异常值往往发生在最小值和最大值处。

四分位距(Inter-Quartile Range)。四分位距就是第三个四分位数和第一个四分位数的差，即：

$$IQR(x)=x_{(n-m+1)}-x_{(m)},m=[n/4]$$

四分位距是稳健的，因为异常值只会微弱影响分位数的取值。另外，IQR是可以标准化的：

$$IQRN=\frac{IQR}{2\times 0.675}$$

1.4 数据位置和波动的同时估计

前面我们介绍了多种不同的稳健统计量，但是存在一些问题。首先，位置模型的M估计很高效，但是无法估计数据波动性；其次，数据波动的稳健估计需要对全体数据进行排序，在数据量大时效率很低。因此，我们考虑把各方的优缺点结合起来，定义新的模型：

$$x_i=\mu+\sigma u_i$$

这个模型被称为位置-分散模型(Location-Disperstion Model)，这里的$\sigma$是随机变量$u_i$的一个尺度参数，同时也是关于$x_i$的一个分散参数。通过这个模型，我们可以同时估计数据中心和数据波动。

设误差$u_i$的分布密度是$f_0$，取$\rho_0=-\log f_0$，那么上述模型的优化目标是：

$$(\hat\mu,\hat\sigma)=\arg\min_{\mu,\sigma}\left\{\frac{1}{n}\sum_{i=1}^n\rho_0\left(\frac{x_i-\mu}{\sigma}\right)+\log\sigma\right\}$$

分别对两个参数求偏导，等价于求解：

$$\begin{aligned} &\sum_{i=1}^n\psi\left(\frac{x_i-\hat\mu}{\hat\sigma}\right)=0\\ &\frac{1}{n}\sum_{i=1}^n\rho_{scale}\left(\frac{x_i-\hat\mu}{\hat\sigma}\right)=\delta \end{aligned}$$

其中：

$$\psi(x)=-\rho_0^{'},\rho_{scale}(x)=x\psi(x),\delta=1$$

这个方程组的结构很重要，在下一节稳健回归模型中的很多估计子就是这样的结构。

1.5 稳健性的测量

在稳健统计的最后，我们探讨一个问题：如何衡量稳健性？

一个比较直接的想法是，对于已有的估计$\hat\mu(x_1,\cdots,x_n)$，我们额外加入一个新的样本点$x_0$，重新计算$\hat\mu(x_1,\cdots,x_n,x_0)$，再对二者进行比较：

$$SC=\hat\mu(x_1,\cdots,x_n,x_0)-\hat\mu(x_1,\cdots,x_n)$$

这里的SC就是敏感性曲线(Sensitivity Curve)，我们可以绘制这样一张敏感性曲线图，即以新样本$x_0$的值为横坐标，SC值为纵坐标，就能查看在加入新样本后统计量的SC的变化：

上图的原有样本来自标准正态分布，Huber的M估计取$k=1.37$，并且分别用SD和MAD估计了分散参数。

SC的方法虽然直观，但是从数学上来看不太严谨。接下来，我们把这种“扰动”再数学化一些。

假设一组样本$x=(x_1,\cdots,x_n)$由$n$个独立同分布$F$的随机变量构成，设估计子$\hat\theta_n=\hat\theta_n(x)$，当样本量趋于无穷时，该估计子依概率收敛于一个渐近值：

$$\hat\theta_n\stackrel{p}{\longrightarrow}\hat\theta_{\infty}(F)$$

这个值被称为估计子在分布F下的渐近值。

• 例如，如果估计子是样本均值，即$\hat\theta_n=\overline x$，那么$\hat\theta_{\infty}(F)=E_F(x)$，也就是分布F的期望；
• 再比如，如果估计子是样本中位数，即$\hat\theta_n=Med(x)$，那么$\hat\theta_{\infty}(F)=F^{-1}(0.5)$，也就是分布F的中位数；
• 再进一步，如果$\hat\theta_n$是位置模型的M估计子，且$\psi$是单调的，那么$\hat\theta_{\infty}(F)$是方程$E_F[\psi(x-\theta)]=0$的解。

有了渐近值的概念后，Hampel在1974年提出了敏感性曲线的一个渐近版本，称之为影响函数(Influence Function, IF)，其定义为：

$$\begin{aligned} IF_{\hat\theta}(x_0,F)&=\lim_{\epsilon\rightarrow0_+}\frac{\hat\theta_{\infty}((1-\epsilon)F+\epsilon\delta_{x_0})-\hat\theta_\infty(F)}{\epsilon}\\ &=\left.\frac{\partial}{\partial\epsilon}\hat\theta_\infty((1-\epsilon)F+\epsilon\delta_0)\right|_{\epsilon\rightarrow0_+} \end{aligned}$$

其中$\delta_{x_0}$是$x_0$处的点质量(point mass)，也就是$P(x=x_0)=1$的分布。

影响函数刻画了加入新样本$x_0$时对样本分布$F$下估计子的扰动，也就是稳健性的衡量。如果影响函数是有界的，说明这个估计子是稳健的。

对于样本均值估计子$\mu$，影响函数为$IF=x_0-\mu$，无界，说明是不稳健的；

对于样本中位数估计子$\theta$，影响函数为$IF=sign(x_0-\theta)/(2f(\theta))$，有界，说明是稳健的；

对于M估计，其估计子就是下面方程的解：

$$\sum_{i=1}^n\Psi(x_i,\hat\theta)=0$$

因此，M估计的影响函数为：

$$IF_{\hat\theta}(x_0,F)=-\frac{\Psi(x_0,\hat\theta_\infty)}{B(\hat\theta_\infty,\Psi)},B(\hat\theta_\infty,\Psi)=\frac{\partial}{\partial\theta}E\Psi(x,\theta)$$

对于截断均值估计，其影响函数为：

$$IF_{\hat\mu}(x_0,F)=\frac{\psi_k(x-\hat\mu_\infty)}{1-2\alpha}$$

其中$\psi$是Huber估计，$k=F^{-1}(1-\alpha)$。

除了影响函数，我们再介绍一个衡量指标：崩溃点(Breaking Point, BP)。崩溃点用于描述一个估计子最多能忍受多少比例的异常值，也就是说，在保证估计子估计不失效的情况下所能容忍的最大污染比例。

严谨来说，我们设$\Theta$是参数$\theta$的可取范围，如果分布为$F$的样本的一个估计子$\hat\theta$的崩溃点是$\epsilon^*$，$\epsilon^*\in(0,1)$，那么对于任意$\epsilon<\epsilon^*$和任意污染分布$G$，存在一个有界闭集$K\subset \Theta$，且$K$远离$\Theta$的边界，使得：

$$\hat\theta_\infty((1-\epsilon)F+\epsilon G)\in K$$

对于位置参数的M估计，我们令$k_1=-\psi(-\infty)$，$k_2=\psi(\infty)$，且$k_1,k_2$有界，那么崩溃点为：

$$\epsilon^*=\frac{\min(k_1,k_2)}{k_1+k_2}$$

例如，均值就是无界的，崩溃点是0；中位数的$\psi(x)=sgn(x)$，$k_1=k_2=1$，那么崩溃点就是$1/2$，这也很容易理解，即当一半数据都被污染时才会影响中位数；截断均值的崩溃点与$\alpha$的选择有关，崩溃点就等于$\alpha$。

对于分散参数的估计，标准差的崩溃点是0，MAD是$1/2$，IQR是$1/4$。

对于位置-分散模型，其崩溃点更复杂一些，但一般来说小于0.5。如果使用Huber的$\psi_k$，那么崩溃点是：

$$\epsilon^*=\min\left(0.5,\frac{0.675}{k+0.675}\right)$$

也就是说，Huber估计的崩溃点是可以调整的。

二、稳健线性回归模型

在第一节中，我们花费了大量篇幅讲解稳健统计，相信你应该对该领域的一些基础知识有一定了解。在这一节，我们要把稳健统计方法应用到线性回归模型当中。

2.1 M-线性回归

我们已经知道，线性回归的框架是：

$$y_i=x_i^T\beta+e_i$$

如果迁移到位置-分散模型框架下，那么$y_i$就是观测值，$x_i^T\beta$就是要拟合的真实值，$e_i$是随机误差。由于M估计是最小化随机误差，于是对于线性回归来说，其M估计为：

$$\hat\beta=\arg\min_\beta\sum_{i=1}^n\rho\left(\frac{r_i(\beta)}{\sigma}\right)$$

其中$r_i(\hat\beta)=y_i-x_i^T\beta$。

假定$\rho$可微分，要求解回归系数，只需对$\beta$求偏导并令其为0，于是得到：

$$\sum_{i=1}^n\psi\left(\frac{r_i(\hat\beta)}{\sigma}\right)x_i=0$$

其中$\psi=\rho_0^{'}$。

如果误差服从正态分布，那么$\psi(x)=x$，求得的回归系数就是最小二乘法估计；如果误差服从双指数分布，那么$\psi(x)= sgn(x)$，求得的回归系数就是最小绝对偏差估计（又叫做L1估计）；如果$\psi$取Huber函数，得到的就是Huber估计。

一般来说，直接求上述方程组的解比较困难，并且要假定$\sigma$已知，所以不太实用。我们将方程组改写为：

$$\sum_{i=1}^n\frac{\psi\left(\frac{r_i(\hat\beta)}{\sigma}\right)}{\frac{r_i(\hat\beta)}{\sigma}}r_i(\hat\beta)x_i=0$$

可以看到，$\psi(x)/x$部分在每一个样本中是一个常数，$r_i(\hat\beta)$是残差，你会发现这个方程组的解其实就是加权最小二乘解，$\psi(x)/x$部分就是权重。因此，我们可以把M估计转化为加权最小二乘解来进行，这样就非常高效。于是我们可以做如下迭代：

1. 初始化。先通过最小二乘法求得$\beta$的初始解$\beta_0$，然后计算残差$r_i=y_i-x_i^T\beta_0$，根据残差的MAD值求得分散参数$\sigma$的初始解$\sigma=MAD(r_i)/0.6745$。
2. 计算标准化残差，即$u_i=r_i/\sigma$。
3. 计算权重，即$w_i=\psi(u_i)/u_i$。$\psi()$可取各种$\rho$函数的导数。
4. 计算加权最小二乘解，即$\beta^{new}=(X^TWX)^{-1}X^TWy$。
5. 更新尺度参数，即$\sigma^{new}=MAD(r_i)/0.6745$。
6. 使用新的$\beta^{new}$计算新的残差$r_i=y_i-x_i^T\beta^{new}$，结合新的$\sigma^{new}$重复2~5步，直至收敛。

另外，尺度参数的更新除了使用标准MAD，还可以使用1.4小节中的：

$$\frac{1}{n}\sum_{i=1}^n\rho_{scale}\left(\frac{r_i}{\sigma}\right)=\delta$$

其中$\delta$一般取标准正态分布随机变量在$\rho$函数中的期望，即$\delta=E\rho_{scale}(Z),Z\sim N(0,1)$。

最后我们讨论M-线性回归的崩溃点。遗憾的是，由于M-线性回归仍然以残差来建模，所以导致M-线性回归的崩溃点很低。如果完全允许高杠杆点的出现（自变量很大并且响应变量很坏），那么崩溃点只有

$$\epsilon^*=\frac{1}{n}$$

如果自变量正常，但响应变量（目标观测值）异常，那么崩溃点会提高：

$$\epsilon^*=\frac{1}{p}$$

综合来说，M-线性回归不算特别稳健的线性回归模型，我们需要寻找崩溃点更高的估计。

2.2 LTS线性回归

最小截断二乘(Least Trimmed Squares, LTS)是Rousseeuw提出的方法，其基本思想还是利用残差，但是只使用残差最小的前$h$个观测点进行回归。也就是说，先计算残差$r_i=y_i-x_i^T\beta$，再对残差进行排序：

$$r_{(1)}\le \cdots\le r_{(n)}$$

取前$h$个最小的残差，求平方和，构成优化目标：

$$\hat\beta_{LTS}=\arg\min_\beta\sum_{i=1}^hr_{(i)}^2$$

一般地，$h$取

$$h=\left[\frac{n+p+1}{2}\right]$$

其中$[\ ]$表示取整。

可以证明，LTS估计的崩溃点可以达到50%，非常稳健。

然而，LTS的优化目标是非光滑、不连续的，所以不能使用梯度法来求解。这里我们介绍一个经典的算法：Fast-LTS。Fast-LTS的基本步骤为：

1. 随机抽取$p+1$个样本，计算最小二乘解$\beta_0=(X_0^TX_0)^{-1}X_0^Ty_0$作为初始解。
2. 计算所有残差，即$r_i=y_i-x_i^T\beta_0$，排序。
3. 取残差最小的前$h$个样本，计算最小二乘解$\beta_H=(X_H^TX_H)^{-1}X_H^Ty_H$，这一步被称为C步(Concentration Step)
4. 依照$\beta_H$计算残差，重复第三步，直至目标函数不再下降，得到局部最优解。
5. 重复1~4步，从而实现多次随机抽取并优化，选择目标函数最小的解作为LTS回归的最终解。

LTS回归的稳健性很好，适合混有大量异常值和污染数据的情形。然而，LTS的搜索效率较低、需要进行数据排序、计算难度很大，所以一般用于进行初步的异常值检测，或者为其他算法提供初始解。

2.3 S-线性回归

无论是M-线性回归还是LTS回归，本质上都在最小化残差的大小。然而，Rousseeuw和Yohai指出，还可以通过最小化残差的尺度参数来进行估计，这样能获得更稳健的结果。S-线性回归就是这样一种策略，通过求解以下方程先估计最优的$\hat\sigma$：

$$\frac{1}{n}\sum_{i=1}^n\rho\left(\frac{r_i(\beta)}{\sigma}\right)=\delta$$

一般情况下，这里的$\rho$函数要取有界的，例如Bisquare函数。

这样做的好处是，通过先优化残差尺度参数来提高稳健性。S-估计的崩溃点也能达到50%。

S-回归的拟合步骤如下：

1. 使用LTS回归获得初始解$\beta_0$。
2. 计算残差$r_i=y_i-x_i^T\beta_0$。
3. 通过求解方程$(1/n)\sum_{i=1}^n\rho(r_i/\sigma)=\delta$得到尺度参数估计值$\hat\sigma$。
4. 计算权重$w_i=\psi(r_i/\hat\sigma)/(r_i/\hat\sigma)$。
5. 使用加权最小二乘进行更新：$\beta_{new}=(X^TWX)^{-1}X^TWy$。
6. 更新残差$r_i=y_i-x_i^T\beta_{new}$，重复3~5步，直至收敛。

S-回归的稳健性很好，效率略优于LTS回归，常常用于生成一个高度稳健的初始值。

2.4 MM-线性回归

既然M-线性回归的计算效率很高，S-线性回归的稳健性非常好，那我们是否可以结合一下，形成一种稳健且高效的回归方法呢？当然可以，这就是本小节要介绍的MM-线性回归。

MM-线性回归的基本思想是，先使用S-估计（或者LTS估计）保证稳健性，再使用M-估计提升求解效率，最后调整尺度参数。这是目前比较实用的稳健线性回归算法。

下面介绍MM-线性回归的拟合步骤。

第一步，使用S-估计或LTS-估计得到初始解$\beta_0,\sigma_0$。

第二步，固定$\sigma_0$，进行M-估计，即：

$$\hat\beta_{new}=\arg\min_\beta\sum_{i=1}^n\rho_1\left(\frac{r_i(\beta)}{\sigma_0}\right)$$

这一步仍然可以使用前面的思路，计算权重$u_i=\psi(r_i/\sigma_0)/(r_i/\sigma_0)$，然后求加权最小二乘解$\hat\beta_{new}=(X^TWX)^{-1}X^TWy$，反复直至收敛。这个过程不要更新σ。

第三步，得到收敛的回归系数$\hat\beta_{MM}$后，再利用尺度估计方程更新尺度，即求解方程：

$$\frac{1}{n}\sum_{i=1}^n\rho_2\left(\frac{r_i(\hat\beta)}{\sigma}\right)=\delta$$

最终得到尺度估计$\hat\sigma_{MM}$。

对于MM-线性回归来说，S-估计保证了其稳健性，而M-估计又保证了运算效率（往往可以达到90%左右），所以成为了很多稳健回归的首选方法。

三、稳健广义线性回归模型简介

既然有稳健的线性回归模型，那肯定少不了稳健的广义线性回归模型（Robust GLM）。然而，稳健的广义线性回归的实现其实较为困难。一方面，广义线性回归的似然结构让经典的M-估计等方法变得不平滑，需要对损失函数进行复杂变换；另一方面，在求解过程中往往遇到的是非凸目标，求解和调参的难度很大。因此，在这一节我们只介绍几种常见的稳健估计方法，对于其他细节问题不再过多展开。

3.1 M-估计

我们在稳健线性回归中了解过M-估计，现在我们把它迁移到广义线性回归当中。对于GLM来说，由于响应变量服从指数分布族，因此其对数似然函数的影响函数是无界的，如果数据出现极端污染，就会导致估计崩溃。因此，我们考虑把$\rho$函数迁移到GLM中，也就是：

$$\hat\beta=\arg\min_\beta\sum_{i=1}^n\rho(r_i),r_i=\frac{y_i-\mu_i}{\sqrt{Var(y_i)}}$$

也就是说，把Pearson残差放进$\rho$函数来优化，这样就能获得稳健的估计。本质上就是求解：

$$\sum_{i=1}^n\psi(r_i)\frac{\partial\mu_i}{\partial\beta}=0$$

M-估计的优点是，简单直观，通用性较强；缺点是崩溃点不会很高，并且依赖$\rho$函数的调参。

GLM的M-估计仍然可以通过IRWLS算法来完成。首先计算工作变量：

$$z_i=g(\mu_i)+g^{'}(\mu_i)(y_i-\mu_i),\mu_i=g^{-1}(x_i^T\beta^{(t)})$$

之后计算权重：

$$w_i=\frac{\psi(r_i)}{r_i}$$

最后更新回归系数：

$$\beta^{(t+1)}=(X^TWX)^{-1}X^TWz$$

3.2 加权似然估计

M-估计是为损失函数添加稳健性，而加权似然(Weighted Likelihood, WL)则是直接改变似然函数的权重，也就是：

$$l=\sum_{i=1}^nw_il_i$$

也就是求解方程：

$$\sum_{i=1}^nw_iu_i(\beta)=0$$

其中$u_i(\beta)$是第$i$个样本的得分统计量。

权重函数的选择同上，一般是：

$$w_i=\frac{\psi(r_i)}{r_i}$$

也就是说，WL方法是对得分统计量加权，也就是对似然函数加权，而3.1小节的M-估计是对残差平方和加权。这样做的好处是，稳健性加强，并且适合计数回归，例如泊松回归或者Gamma回归中响应值异常大的情况。然而，WL方法对高杠杆点的处理仍然较弱，权重设计和理论推导上较难。

3.3 散度估计

散度(Divergence)是描述矢量场在空间中某点的 “发散” 或 “汇聚” 的程度，是一个标量，如今在信息论中有非常重要的应用。这里我们不展开，我们简单介绍两种运用散度思想实现稳健GLM的方法。

第一种称为密度幂散度(Density Power Divergence, DPD)。DPD是Basu等人于1998年提出的方法，其定义是，对于模型概率密度分布$f$和真实分布$g$来说，构造散度：

$$D_\alpha(g,f)=\int f^{1+\alpha}(y)dy-\frac{1+\alpha}{\alpha}\int g(y)f^\alpha(y)dy+\frac{1}{\alpha}\int g^{1+\alpha}(y)dy$$

如果是GLM模型，那么目标函数可以写作：

$$H_\alpha(\beta)=\int f_\beta^{1+\alpha}(y)dy-\frac{1+\alpha}{\alpha n}f_\beta^\alpha(y_i)$$

把指数分布族带入$f_\beta(y)$求解优化即可。

第二种称为γ散度(γ/Gamma-Divergence)。γ散度由Fujisawa和Eguchi分别于2008年和2013年提出并优化，其定义是，对于模型概率密度分布$f$和真实分布$g$来说，构造散度：

$$D_\gamma(g,f)=\frac{1}{\gamma(\gamma+1)}\left[\int f^{1+\gamma}(y)dy-(1+\gamma)\int g(y)f^\gamma(y)dy+\int g^{1+\gamma}(y)dy\right]$$

γ散度和DPD在形式上类似，但是γ散度具有额外的归一化性质，对极端响应值的抑制更强，更适用于重尾分布。

如果是GLM模型，那么γ散度估计的目标函数为：

$$H_\gamma(\beta)=-\frac{1}{\gamma}\log\left(\sum_{i=1}^n f_\beta^\gamma(y_i)\right)+\frac{1}{\gamma+1}\log\left(\int f_\beta^{\gamma+1}(y)dy\right)$$

以上就是本讲的全部内容，相信你已经对稳健估计和稳健回归有一定的了解了。然而，稳健估计理论和稳健回归方法远远不止这些，也不可能把这些内容全部浓缩在这一篇讲义当中。因此，本讲的目的是，能够介绍一些基本的思想和方法，当后续的新内容再提及时也会让你感到熟悉。