Robustness regression

稳健回归（Robustness regression）

 https://www.jianshu.com/p/cdbfedcaceaf

最小二乘法的弊端

之前文章里的关于线性回归的模型，都是基于最小二乘法来实现的。但是，当数据样本点出现很多的异常点（outliers），这些异常点对回归模型的影响会非常的大，传统的基于最小二乘的回归方法将不适用。

比如下图中所示，数据中存在一个异常点，如果不剔除改点，适用OLS方法来做回归的话，那么就会得到途中红色的那条线；如果将这个异常点剔除掉的话，那么就可以得到图中蓝色的那条线。显然，蓝色的线比红色的线对数据有更强的解释性，这就是OLS在做回归分析时候的弊端。

当然，可以考虑在做回归分析之前，对数据做预处理，剔除掉那些异常点。但是，在实际的数据中，存在两个问题：

1. 异常点并不能很好的确定，并没有一个很好的标准用于确定哪些点是异常点
2. 即便确定了异常点，但这些被确定为异常的点，真的是错误的数据吗？很有可能这看似异常的点，就是原始模型的数据，如果是这样的话，那么这些异常的点就会带有大量的原始模型的信息，剔除之后就会丢失大量的信息。

再比如下面这幅图，其中红色的都是异常点，但是很难从数据中剔除出去。

稳健回归

稳健回归（Robust regression），就是当最小二乘法遇到上述的，数据样本点存在异常点的时候，用于代替最小二乘法的一个算法。当然，稳健回归还可以用于异常点检测，或者是找出那些对模型影响最大的样本点。

Breakdown point

关于稳健回归，有一个名词需要做解释：Breakdown point，这个名词我并不想翻译，我也没找到一个很好的中文翻译。对于一个估计器而言，原始数据中混入了脏数据，那么，Breakdown point 指的就是在这个估计器给出错误的模型估计之前，脏数据最大的比例 α ，Breakdown point 代表的是一个估计器对脏数据的最大容忍度。

举个简单的例子：有 n 个随机变量， (X1,X2,…,Xn) ， 其对应的数据为 (x1,x2,…,xn) ，那么，我么可以求出这 n 个随机变量的均值：

 

X¯¯¯=X1+X2+⋯+Xnn

 

这个均值估计器的Breakdown point 为0，因为使任意一个 xi 变成足够大的脏数据之后，上面估计出来的均值，就不再正确了。

毫无疑问，Breakdown point越大，估计器就越稳健。

Breakdown point 是不可能达到 50% 的，因为如果总体样本中超过一半的数据是脏数据了，那么从统计上来说，就无法将样本中的隐藏分布和脏数据的分布给区分开来。

本文主要介绍两种稳健回归模型：RANSAC（RANdom SAmple Consensus 随机采样一致性）和Theil-Sen estimator。

RANSAC随机采样一致性算法

RANSAC算法的输入是一组观测数据（往往含有较大的噪声或无效点），它是一种重采样技术（resampling technique），通过估计模型参数所需的最小的样本点数，来得到备选模型集合，然后在不断的对集合进行扩充，其算法步骤为：

1. 随机的选择估计模型参数所需的最少的样本点。
2. 估计出模型的参数。
3. 找出在误差 ϵ 内，有多少点适合当前这个模型，并将这些点标记为模型内点
4. 如果内点的数目占总样本点的比例达到了事先设定的阈值 τ ，那么基于这些内点重新估计模型的参数，并以此为最终模型， 终止程序。
5. 否则重复执行1到4步。

RANSAC算法是从输入样本集合的内点的随机子集中学习模型。

RANSAC算法是一个非确定性算法（non-deterministic algorithm），这个算法只能得以一定的概率得到一个还不错的结果，在基本模型已定的情况下，结果的好坏程度主要取决于算法最大的迭代次数。

RANSAC算法在线性和非线性回归中都得到了广泛的应用，而其最典型也是最成功的应用，莫过于在图像处理中处理图像拼接问题，这部分在Opencv中有相关的实现。

从总体上来讲，RANSAC算法将输入样本分成了两个大的子集：内点（inliers）和外点（outliers）。其中内点的数据分布会受到噪声的影响；而外点主要来自于错误的测量手段或者是对数据错误的假设。而RANSAC算法最终的结果是基于算法所确定的内点集合得到的。

下面这份代码是RANSAC的适用实例：

# -*- coding: utf-8 -*-

""" author : duanxxnj@163.com time : 2016-07-07-15-36 """

import numpy as np
import time
from sklearn import linear_model,datasets

import matplotlib.pyplot as plt

# 产生数据样本点集合
# 样本点的特征X维度为1维，输出y的维度也为1维
# 输出是在输入的基础上加入了高斯噪声N（0,10）
# 产生的样本点数目为1000个

n_samples = 1000
X, y, coef = datasets.make_regression(n_samples=n_samples,
                                      n_features=1,
                                      n_informative=1,
                                      noise=10,
                                      coef=True,
                                      random_state=0)

# 将上面产生的样本点中的前50个设为异常点（外点）
# 即：让前50个点偏离原来的位置，模拟错误的测量带来的误差
n_outliers = 50
np.random.seed(int(time.time()) % 100)
X[:n_outliers] = 3 + 0.5 * np.random.normal(size=(n_outliers, 1))
y[:n_outliers] = -3 + 0.5 * np.random.normal(size=n_outliers)

# 用普通线性模型拟合X，y
model = linear_model.LinearRegression()
model.fit(X, y)

# 使用RANSAC算法拟合X,y
model_ransac = linear_model.RANSACRegressor(linear_model.LinearRegression())
model_ransac.fit(X, y)
inlier_mask = model_ransac.inlier_mask_　　　　#返回不是异常点的值， 布尔类型数组
outlier_mask = np.logical_not(inlier_mask)

# 使用一般回归模型和RANSAC算法分别对测试数据做预测
line_X = np.arange(-5, 5)
line_y = model.predict(line_X[:, np.newaxis])
line_y_ransac = model_ransac.predict(line_X[:, np.newaxis])

print "真实数据参数：", coef
print "线性回归模型参数：", model.coef_
print "RANSAC算法参数： ", model_ransac.estimator_.coef_


plt.plot(X[inlier_mask], y[inlier_mask], '.g', label='Inliers')
plt.plot(X[outlier_mask], y[outlier_mask], '.r', label='Outliers')
plt.plot(line_X, line_y, '-k', label='Linear Regression')
plt.plot(line_X, line_y_ransac, '-b', label="RANSAC Regression")
plt.legend(loc='upper left')
plt.show()

　　

运行结果为：

真实数据参数： 82.1903908408
线性回归模型参数： [ 55.19291974]
RANSAC算法参数：  [ 82.08533159]

Theil-Sen Regression 泰尔森回归

Theil-Sen回归是一个参数中值估计器，它适用泛化中值，对多维数据进行估计，因此其对多维的异常点（outliers 外点）有很强的稳健性。

一般的回归模型为：

 

y=α+βx+ϵ

 

其中， α,β 模型的参数，而 ϵ 为模型的随机误差。

Theil-Sen回归则是这么处理的：

 

β=Median{yi−yjxi−xj:xi≠xj,i<j=1,…,n}

 

在实践中发现，随着数据特征维度的提升，Theil-Sen回归的效果不断的下降，在高维数据中，Theil-Sen回归的效果有时甚至还不如OLS（最小二乘）。

在之间的文章《线性回归》中讨论过，OLS方法是渐进无偏的，Theil-Sen方法在渐进无偏方面和OLS性能相似。和OLS方法不同的是，Theil-Sen方法是一种非参数方法，其对数据的潜在分布不做任何的假设。Theil-Sen方法是一种基于中值的估计其，所以其对异常点有更强的稳健性。

在单变量回归问题中，Theil-Sen方法的Breakdown point为29.3%，也就是说，Theil-Sen方法可以容忍29.3%的数据是outliers。

# -*- coding: utf-8 -*-

""" @author : duanxxnj@163.com @time ；2016-07-08_08-50 Theil-Sen 回归 本例生成一个数据集，然后在该数据集上测试Theil-Sen回归 """

print __doc__

import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, TheilSenRegressor,\
                                 RANSACRegressor

estimators = [('OLS', LinearRegression()),
              ('Theil-Sen', TheilSenRegressor())]

# 异常值仅仅出现在y轴
np.random.seed((int)(time.time() % 100))
n_samples = 200

# 线性模型的函数形式为： y = 3 * x + N(2, .1 ** 2)
x = np.random.randn(n_samples)
w = 3.
c = 2.
noise = c + 0.1 * np.random.randn(n_samples)
y = w * x + noise

# 加入10%的异常值,最后20个值称为异常值
y[-20:] += -20 * x[-20:]

X = x[:, np.newaxis]
plt.plot(X, y, 'k+', mew=2, ms=8)
line_x = np.array([-3, 3])

for name, estimator in estimators:
    t0 = time.time()
    estimator.fit(X, y)
    elapsed_time = time.time() - t0
    y_pred = estimator.predict(line_x.reshape(2, 1))
    plt.plot(line_x, y_pred, label='%s (fit time: %.2fs)'
             %(name, elapsed_time))

plt.axis('tight')
plt.legend(loc='upper left')

plt.show()

　　

稳健回归能够在有损坏的数据下拟合出一个回归模型：异常值或模型中的误差。

受损的y

1.1.15.1. 不同的场景和有用的概念

在处理由异常值损坏的数据时有几点我们需要注意：

• 误差值在 X 还是 y ?

误差值在y	误差值在X
误差值在y	误差值在X

• 异常值的比值与误差的幅度

小误差	大误差
小误差	大误差

稳健拟合的一个重要概念是抛锚点(breakdown point)：即模型会为了拟合外部的少量数据，而开始丢失内部的数据。
要时刻注意这一点，稳健回归在拟合高维度的数据（大量的 n_feature）中的设置是很难的。所以在这种设置下稳健模型可能不会进行工作。

权衡：哪种估量器？
Scikit-learn提供了三种稳健回归估量器：RANSAC, Theil Sen 与 HuberRegressor

• 除非样本的数据量非常大，否则一般情况下 HuberRegressor 会比 RANSAC 与 Theil Sen 要快得多，即 样本数量 >> 特征数量。这是因为 RANSAC 和 Theil Sen 是拟合数据集中的子集。
• RANSAC 比 Theil Sen 要快，且会随着样本的数量增加而表现的更好。
• RANSAC 在处理大量于 y 里的异常值（同时这也是异常值最普遍出现的形式）表现良好。
• Theil Sen 能够比较好的应对于 X 里存在中等数量异常值的情况。但是这一特性会在高维中失效。
  也就是说，如果不知道要用那种估量器时，用 RANSAC 就对了。

1.1.15.2. RANSAC: 随机抽样一致

RANSAC（随机抽样一致）从完整数据中的随机子集来拟合模型。
RANSAC是一个非确定性算法他只会在一定概率下产生合理的结果，这概率取决于迭代的次数（可以通过 max_trals 参数设置）。其典型的应用场景是用来解决线性与非线性回归问题，并且它也在摄影测量计算机视觉领域中非常流行。
这个算法会把完整的输入拆分成一个集合集，同时这也会使得在每个集合中多多少少受到噪音与异常值的影响。例如。由错误测量引起的噪音或由无效的假设所收集的数据。最终产生的模型也是根据集合来决定。

RANSAC回归

1.1.15.2.1. 算法的细节

该算法在每次迭代中都会执行以下步骤：

1. 从原始数据中选出 min_samples 个随机样本并且检查数据集是否有效（ is_data_valid ）。
2. 根据随机出来的子集进行拟合。（base_estimator.fit），并且对估量出的模型进行判断其是否有效（is_model_valid）。
3. 通过计算估量模型的残差(base_estimator.predict(X) - y)将所有数据进行分类（内点与离群点）。数据的残差比 residual_threshold 要小的就作为内点。
4. 如果内点样本已达到最大值，就保存拟合后的模型为合适模型。以防止当前估计出的模型拥有相同数量的内点。不过如果模型的分数更高也会被认为是更好的模型。

这些步骤的最大执行次数（max_trials）或直到满足特殊停止条件之一（详见 stop_n_inliers 和 stop_score）后停止。最后会使用先前所有估计出的合适模型来验证所有内点样本来决定最合适的模型。

is_data_valid 和 is_model_valid 函数允许用来识别和拒绝随机子样本的退化组合。如果估计出的模型不再需要鉴定退化样本，那么就应该在拟合模型之前调用is_data_valid，使的模型能够获得一个更好的计算性能。

示例

• 使用RANSAC的稳健线性模型
• 稳健线性估计器的拟合过程

引用

• https://en.wikipedia.org/wiki/RANSAC
• “Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography” Martin A. Fischler and Robert C. Bolles - SRI International (1981)
• “Performance Evaluation of RANSAC Family” Sunglok Choi, Taemin Kim and Wonpil Yu - BMVC (2009)

1.1.15.3. Theil-Sen估量器：普通中值估量器

TheilSenRegressor 估量器使用了多维数据中的中值概括。因此对于有多变量异常值时仍旧是稳健的。但无论如何都要注意的是，估量器的鲁棒性会随着问题的维度的增加而快速减少。所以在高维度的情况其在失去鲁棒性的情况，并且效率可能还不如OLS高。

示例

• Theil-Sen回归
• 稳健线性回归估计器的拟合过程

引用

• https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator

1.1.15.3.1. 理论限制

TheilSenRegressor作为一个无偏差值的估计器而言，在拟合效率上与OLS相似。与OLS相反，Theil-Sen是一个非参数方法，这意味着它不会对地底层分布的数据作任何假设。因为 Theil-Sen 是一种基于中值的估量器，所以它对于损坏的数据更具鲁棒。在单变量设置里，在简单线性回归的情况下Theil-Sen的击破点率为29.3%，这意味着它可以容忍高达29.3％的任意损坏的数据。

受损y

在scikit-learn中，TheilSenRegressor 的实现是根据多元线性回归的空间中值这一多维中值的概念[7][8]。

在时间和空间复杂度方面，Theil-Sen的尺度根据

Tenil-Sen的尺度

使得在大量数据与特征的情况下，不可能将所有的子集样本应用到该问题中。因此子集的量级可从考虑所有的随机子集的集合来限制时间和空间复杂度。

示例

• Theil-Sen回归

引用

1. Xin Dang, Hanxiang Peng, Xueqin Wang and Heping Zhang: Theil-Sen Estimators in a Multiple Linear Regression Model.
2. T. Kärkkäinen and S. Äyrämö: On Computation of Spatial Median for Robust Data Mining.

1.1.15.4. Huber回归

HuberRegressor 与 Ridge的不同在于前者应用线性损失到样本时，是将其归类为异常值。样本的分类为内点是根据其的绝对错误是否低于一个确定的阈值。其跟 TheilSenRegressor 的 RANSACRegressor 不同之处在于它并不会忽略异常值，相反而是给予他们一个很小的权重值。

Huber与Ridge的对比

HuberRegressor 的最小化损失函数的公式为：

Huber的最小化损失公式

其中：

Huber的最小化损失公式中的H

建议设置参数 epsilon 的值为1.35以获得95%统计效率。

1.1.15.5. 注意事项

在损失设置上， HuberRegressor 与 SGDRegressor 的不同之处在于：

• HuberRegressor 是缩放不变的。当 epsilon 被设置了后，通过任意值缩放 X 与 y 都不会改成其异常值的鲁棒性。而 SGDRegressor 在缩放 X 与 y 后仍需要重新设置 epsilon 的值。
• HuberRegressor 能够充分使用少量样本的数据，而 SGDRegressor 为了产生一致的鲁棒性则需要多次传递数据。

示例

• 在强异常值下 Hube回归 vs Ridge回归

引用

• Peter J. Huber, Elvezio M. Ronchetti: Robust Statistics, Concomitant scale estimates, pg 172

同样的，该估计器与R语言实现的稳健回归(http://www.ats.ucla.edu/stat/r/dae/rreg.htm)也是有不同的地方。因为R语言的实现是基于加权最小二乘，根据残差大于特定阈值多少来给予每个样本的权重。