做回归任务的人，多少都被异常值坑过。一组数据里混进几个离谱的点，最小二乘的拟合线就会整个歪掉。稳健回归（Robust Regression）就是专门解决这类问题的：Theil-Sen回归、RANSAC回归、Huber回归，是三种最常见也最实用的稳健回归方案。这篇文章会把这三种方法的原理、选型逻辑、完整代码和调参经验一次讲透，适合做数据建模、质量监控、指标清洗的朋友参考。
1. 先从最小二乘的“脆弱性”说起
1.1 一个异常点就能让拟合“崩掉”
最小二乘（Ordinary Least Squares, OLS）是我们机器学习路上最先接触的回归方法。它的优化目标是让“预测值与真实值的差值的平方和”最小。问题恰恰出在这个“平方”上：残差是2，平方后是4；残差是10，平方后是100；残差是50，平方后是2500。异常值只要离群一点，它的残差平方就会把整个损失函数的重心拽过去。
我经常用一句话解释这件事：最小二乘是“讨好所有人”，但异常值不是“人”，它是“熊”。模型为了减少一个熊的怒吼，宁愿牺牲一大片普通样本的拟合精度。举例来说，100个样本分布在一条斜率为3的直线附近，只要混进1个偏离特别大的点，拟合出来的斜率可能直接掉到1以下。这不是夸张，而是平方损失的特性决定的。
很多同学遇到这种情况，第一反应是肉眼找异常值、手动剔除。数据量小的时候可以这么干，但现实中的数据采集过程充满不确定性，传感器偶发故障、人工录入出错、指标口径调整，都会产生异常点。你不可能每来一批数据都手动筛选一遍。真正可落地的方式，是让算法本身对异常值“不敏感”，这就是稳健回归存在的意义。
1.2 稳健回归的共同出发点
既然问题的根源是“平方把异常值的影响放大了”，那么解决思路无非三条方向：
第一，不用均值，改用中位数。均值会被极端值拉偏，中位数对极端值免疫，Theil-Sen回归就建立在这个朴素但扎实的想法上。
第二，先筛选、后拟合。与其让异常值参与计算再从结果里纠正，不如先判断哪些样本是“正常样本”，只拿正常样本做拟合。RANSAC就是典型的“筛选派”。
第三，改损失函数。既然平方损失对残差大的点惩罚太重，那就设计一个折中的损失：残差小的时候用平方，残差大的时候改用线性增长，甚至封顶。Huber回归就是这个思路的经典代表。
这三种思路分别对应了“改统计量”、“改采样策略”、“改损失函数”三个层面的设计。理解了这一点，你再看后面的代码和参数，就不会觉得它们是一堆需要死记硬背的配置项，而是一套有明确意图的策略。
2. 三种稳健回归的原理与选型逻辑
2.1 Theil-Sen回归：把所有斜率的“中位数”当成答案
Theil-Sen回归在统计学里被称为“中位数斜率估计”。它的核心操作极简：对于二维数据，任意两个点可以确定一条直线，得到一个斜率。如果有n个样本，理论上可以得到 n(n-1)/2 个斜率，Theil-Sen回归取的就是这些斜率的“中位数”。
中位数为什么稳健？因为哪怕有一半的点都出了问题，只要正常点仍然占多数，斜率的“中位数”就不会被极端斜率带走。这比最小二乘用均值来估计参数要牢固得多。
但使用Theil-Sen回归有一个显眼的代价：当样本量较大时，两两配对的计算量接近O(n²)，如果样本量超过几千，计算成本会变得难以接受。sklearn中的实现提供了max_subpopulation和下采样机制来缓解这个问题，但本质上，它更适合小规模数据、低维特征的场景。
2.2 RANSAC回归：先找到“干净子集”再拟合
RANSAC的全称是“随机采样一致性”（Random Sample Consensus），它和我上面说的“筛选派”思路完全一致。它的工作流程分四步：
- 从样本中随机抽取一个最小子集（比如直线拟合至少需要2个点），用这个子集拟合一棵临时模型。
- 计算所有样本到该模型的残差，残差小于阈值
residual_threshold的样本被称为“内点”，其余是“外点”。 - 记录当前模型的内点数量，并只用内点重新拟合模型。
- 重复上述过程若干次，最终保留内点数量最多、且模型质量最好的一组结果。
说得直白一点，RANSAC的做法是“先假设一个模型，看谁能接受它；谁接受它，谁就是自己人；然后用自己人重新拟合模型”。这个流程非常直观，而且它对异常值比例很宽容，哪怕有一半数据是噪声，只要噪声不是严重偏离到大规模聚集，它依然有机会找到正确的模型。
RANSAC最大的特点其实是“随机性”。它依赖随机采样来寻找最优模型，意味着同样的数据和参数，换一个随机种子，结果可能不同。所以实际使用中，设定random_state是必须养成的习惯，否则你复现不出实验结果。
2.3 Huber回归：给损失函数加一个“封顶”
Huber回归属于损失函数派的代表，它试图在“最小二乘的精度”和“绝对损失的稳健”之间找平衡。它的损失函数长这样：
- 当残差绝对值小于等于阈值 ε（epsilon）时，按平方损失计算；
- 当残差绝对值大于 ε 时，按线性损失计算。
用生活化的类比理解：对Huber来说，误差小的时候，它像严格的老师，用平方惩罚，希望模型精雕细琢；误差大的时候，它变成宽容的管理者，只按线性比例处罚，不让极端残差通过平方放大来左右模型。
Huber回归非常适合那些“异常值几乎不可避免，但正常数据本身又服从类正态分布”的场景。它既保留了正常数据下的高精度，又对异常值有不错的耐受能力。其核心参数epsilon就是“平方区域”和“线性区域”之间的一个分界线，决定模型对异常值的容忍程度。
2.4 选型对照：该用哪一种
三种方法各有各的舒适区，不能一概而论。我曾经把三种方法放在同一组脏数据上进行对比，结果是：异常值比例不高时，Theil-Sen和RANSAC表现接近；异常值比例高时，RANSAC更稳；当数据本身噪声很大但分布连续时，Huber的预测偏差最小。
| 方法 | 核心思想 | 适合场景 | 缺点与注意点 |
|---|---|---|---|
| Theil-Sen | 配对斜率的单位数 | 小样本、低维、异常比例中等 | 样本量大时计算开销大 |
| RANSAC | 反复采样，内点子集拟合 | 异常比例高、数据点较多 | 结果依赖随机种子，需调阈值 |
| Huber | 分段损失，平方与线性结合 | 正常数据近似正态、伴随连续异常 | epsilon需要调，参数敏感 |
我自己在实践中常用一个简单原则：如果只想给线性回归加一道“保险杠”，选Huber；如果确定数据里混入了大量脏点，优先试RANSAC；如果计算资源充足但样本量不大，Theil-Sen往往能够给出非常稳健、可解释的回归结果。
3. 全代码实战：从模拟数据到模型效果对比
3.1 构造带异常值的数据集
实践出真知，我们直接构造一组带异常值的数据来对比。这里用Python和scikit-learn。
import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression, TheilSenRegressor, RANSACRegressor, HuberRegressor from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score # 固定随机种子，保证结果可复现 np.random.seed(42) n = 80 X = np.random.uniform(0, 10, n).reshape(-1, 1) # 真实关系：y = 6 + 3x + 噪声 true_slope = 3.0 true_intercept = 6.0 y = true_intercept + true_slope * X.ravel() + np.random.normal(0, 1.0, n) # 随机挑10个样本，注入大幅异常值 outlier_idx = np.random.choice(n, size=10, replace=False) y[outlier_idx] = np.random.uniform(-25, 50, size=10)我这里故意让真实关系非常简单：斜率为3，截距为6。这样做的好处是，后续你一眼就能看出哪种模型的参数更接近真实值。注入的异常值范围是-25到50，对于正常位于5~35区间内的数据来说，这些异常点足够“刺眼”。
如果你担心RANSAC的随机性导致每次结果漂移，可以把np.random.seed和下面每个模型的random_state都固定。我把所有模型的随机种子都设为42，方便你直接复现。
df = pd.DataFrame({'x': X.ravel(), 'y': y}) print(df.head()) print("异常值索引：", sorted(outlier_idx))3.2 最小二乘基线：看看它到底多容易翻车
先跑一个普通的线性回归做基线。这个对比很有必要，因为只有看到最小二乘“翻车”的程度，你才能直观感受稳健回归的价值。
ols = LinearRegression() ols.fit(X, y) print("OLS 斜率：", round(ols.coef_[0], 4)) print("OLS 截距：", round(ols.intercept_, 4)) print("OLS R2：", round(r2_score(y, ols.predict(X)), 4))你在自己机器上运行后会发现，斜率已经和真实的3.0差了很远，甚至可能变成负数。这个结果非常符合我对“平方损失被异常值直接支配”的判断。在10%的异常值比例下，最小二乘的参数估计就已经严重劣化。
3.3 三行代码跑通三种稳健回归
scikit-learn把三种稳健回归都封装成了和LinearRegression同接口的类，调用方式非常统一。这也是为什么我推荐先从sklearn入手，先把方法论跑通，再按需去读底层实现。
# Theil-Sen回归 theil_sen = TheilSenRegressor(random_state=42) theil_sen.fit(X, y) # RANSAC回归 ransac = RANSACRegressor(random_state=42) ransac.fit(X, y) # Huber回归 huber = HuberRegressor(epsilon=1.35, max_iter=1000, alpha=0.0001) huber.fit(X, y)需要注意的一点是，HuberRegressor在alpha=0时会给出警告，因为正则化项为0时模型行为可能不稳定，我建议保留非常小的alpha值，比如0.0001，这样既不影响稳健性，也让求解器走得稳。
3.4 效果评估与可视化输出
模型的优劣不能只靠肉眼看，我们用多个指标统一评估。因为数据里含异常值，我会重点关注 MAE（平均绝对误差），它对异常值的惩罚相对温和，更能体现稳健模型在“普通样本上的拟合质量”。
models = { 'OLS': ols, 'TheilSen': theil_sen, 'RANSAC': ransac, 'Huber': huber, } metrics = [] for name, model in models.items(): pred = model.predict(X) metrics.append({ 'Model': name, 'Slope': round(model.coef_[0], 4), 'Intercept': round(model.intercept_, 4), 'MAE': round(mean_absolute_error(y, pred), 4), 'RMSE': round(mean_squared_error(y, pred) ** 0.5, 4), 'R2': round(r2_score(y, pred), 4), }) result_df = pd.DataFrame(metrics) print(result_df)如果你运行上面这段代码，大概率会看到：OLS的斜率严重偏离3.0，而三种稳健回归的斜率和截距都更接近真实值。其中RANSAC因为显式剔除了离群点，拟合线往往最贴近真实直线。
再画一张诊断图，把散点、真实直线和各模型拟合线画到一起：
x_line = np.linspace(0, 10, 100).reshape(-1, 1) plt.figure(figsize=(10, 6)) plt.scatter(X.ravel(), y, alpha=0.6, label='样本点', color='gray') plt.plot(x_line, true_intercept + true_slope * x_line.ravel(), color='black', linestyle='--', linewidth=2, label='真实直线') for name, model in models.items(): plt.plot(x_line, model.predict(x_line), label=name, linewidth=2) plt.xlabel('x') plt.ylabel('y') plt.legend() plt.title('不同回归模型在含异常值数据上的表现') plt.show()图中最扎眼的一定是OLS那条线，它为了“照顾”右上角或下方的异常点，把整个斜率拉变形了。稳健回归的三条线则基本重叠在真实直线附近。这个可视化实验我做过很多次，每次都能让团队成员直观理解为什么不能直接用普通线性回归处理生产数据。
4. 参数调优与使用陷阱
4.1 Theil-Sen参数：当样本量大了，必须缩n_subsamples
Theil-Sen回归没有太多需要调参的地方，核心关注点是n_subsamples和max_subpopulation。
n_subsamples表示实际用于估计斜率的样本数。默认值是样本总量的一半。当样本量较大时，如果你想控制计算量，可以显式设置一个较小的值，比如全部样本的20%。但要注意，n_subsamples至少要为特征数加1，否则无法拟合成一条有效直线。
max_subpopulation则是对配对斜率数量的一种“上限截断”。sklearn文档中它默认是1e6，意思是如果可能的配对数量超过这个上限，算法会采用随机抽样的方式去找一部分配对，而不是枚举全部组合。如果你希望追求更高精度的估计，可以把这个值调大，但代价是计算时间明显上升。
theil_sen = TheilSenRegressor( n_subsamples=20, max_subpopulation=1e5, random_state=42 )4.2 RANSAC参数：residual_threshold是最影响结果的旋钮
RANSAC的参数相对多一些，但最重要的无疑是residual_threshold。它决定了“内点”的判定标准：残差小于这个阈值的样本才被认为是正常样本。阈值设小了，会把大量正常样本误判为外点，导致模型只由一小撮点拟合，方差变大；阈值设大了，异常点又会被当成内点混进拟合，稳健性下降。
一个实用的做法是先用普通线性回归做一个初步拟合，然后计算残差的标准差，把阈值设置为残差标准差的整数倍。例如：
import numpy as np from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_absolute_error # 先用OLS粗拟合，估计残差尺度 tmp_ols = LinearRegression() tmp_ols.fit(X, y) residuals = np.abs(y - tmp_ols.predict(X)) threshold = np.percentile(residuals, 95) # 或使用若干倍标准差 ransac = RANSACRegressor( min_samples=5, residual_threshold=threshold, max_trials=200, random_state=42 ) ransac.fit(X, y)此外，min_samples表示每次采样的最小拟合样本数。对于二维直线拟合，理论上设为2即可，但在实际使用中我建议设置为总样本的10%左右，让初始拟合段更稳定。max_trials是最大迭代次数，默认值是100，当异常值比例较高时可以提升到300甚至更多。
4.3 Huber参数：epsilon决定“容忍度”
Huber回归的关键参数就一个：epsilon。它控制平方损失向线性损失切换的阈值。sklearn 默认值是1.35。
epsilon越小，模型对异常值越“冷酷”，稍微偏离预期就会改用线性惩罚；epsilon越大，模型越接近普通最小二乘，因为大部分残差都会落在平方区域内。实际调参时，我通常从默认值1.35出发，用交叉验证或网格搜索在0.5~3.0之间扫一遍。
from sklearn.model_selection import GridSearchCV param_grid = {'epsilon': [1.0, 1.2, 1.35, 1.5, 2.0, 2.5, 3.0]} huber_cv = GridSearchCV( HuberRegressor(alpha=0.0001), param_grid=param_grid, cv=5, scoring='neg_mean_absolute_error' ) huber_cv.fit(X, y) print("最优 epsilon：", huber_cv.best_params_)注意，HuberRegressor默认没有random_state参数，因为它不依赖随机采样，而是通过迭代优化求解。如果你在可视化时发现Huber的效果比RANSAC差一点，别意外，Huber的定位本来就不是“剔除”异常值，而是“弱化”异常值的影响，二者的行为逻辑不同。
4.4 交叉验证与随机数种子：别让结果忽好忽坏
无论用哪种稳健回归，有两点建议都值得遵守：
第一，固定随机种子。RANSAC和Theil-Sen在实现中都可能用到随机抽样，如果不固定random_state，同一个脚本每次跑出来的模型都不一样，这对上线和回溯都是灾难。
第二，用交叉验证评估，不要只看单次拟合的R²。稳健回归在训练集上表现优异不代表泛化能力好，尤其RANSAC，它在训练集上可能因为只选了部分内点而“过拟合”了内点模式。我习惯把数据按时间或随机方式拆分多次，比较不同折上的MAE变化。
5. 常见报错与场景化选型建议
5.1 五种常见报错与排查速查表
| 报错或运行异常 | 常见原因 | 解决方案 |
|---|---|---|
| “n_subsamples must be greater than n_features” | n_subsamples设置过小 | 确保 n_subsamples >= 特征数 + 1 |
| RANSAC运行时间过长 | max_trials过大或每次拟合样本多 | 减少max_trials，设置stop_n_inliers提前终止 |
| HuberRegressor警告alpha=0不可靠 | 设置了alpha=0 | 改成0.0001这类极小值 |
| 结果每次都不同 | 未固定random_state | 为模型设置固定random_state |
| Theil-Sen在数据量大的时候特别慢 | n(n-1)/2配对斜率太大 | 设置max_subpopulation小一些，或改用RANSAC |
其中，RANSAC运行时间过长是新手最容易忽略的问题。stop_n_inliers参数可以帮助提前终止迭代——如果已经找到了一个内点数量达到某个比例的模型，就停止继续搜索。比如我设定stop_n_inliers=int(0.75 * len(X))，结合stop_score=9999，能明显缩短耗时。
5.2 现实场景：传感器噪声、业务指标、异常流量清理
理论说再多，最终还是要回到业务场景。这里分享几个我曾经处理过的典型情形：
某物流公司做运费预测时，数据里混入了大量“测试订单”，这些订单的运费和重量完全没有相关性。一开始用普通回归做预测，误差大得离谱。后来改用RANSAC回归，把测试订单当成外点筛掉，模型才稳定下来。
某质量监控平台上，传感器偶尔会抖动，产生瞬间的尖峰数据。这种数据的异常比例不高，但数值偏离极大。用Huber回归来做特征校准非常合适，它既不会丢掉真实波动，又不会被偶尔的尖峰带偏。
某电商业务的日活趋势分析，偶尔会有“秒杀活动”导致数据暴涨，这些点不是错误，是真实事件，但对趋势建模来说，它们依然属于“需要弱化”的极端值。我通常会先用Theil-Sen做趋势基线，再看那些点偏离基线的程度，以此识别活动效应。Theil-Sen的优势在这里体现得很明显，参数解释性强，不需要随机种子带来的不确定性。
5.3 几条踩坑心得
写这篇内容的过程中，我回顾了自己用稳健回归踩过的一些坑，整理成几条实在的体会：
一是不要对“RANSAC永远最稳健”抱有不切实际的期望。当异常值本身聚集成了某种“伪模式”时，RANSAC可能把这种伪模式误认为主模式。举个极端例子，如果脏数据也沿另一条直线集中分布，RANSAC会陷入“两条线之争”，结果可能取决于随机采样路径。
二是Theil-Sen在小样本下非常稳，但维度升高后效果快速衰减。它处理多元回归时的核心思路仍然是“两两配对”，但高维空间里点对之间的斜率组合爆炸式增长，而且中位数对多变量关系表达力有限。我一般只在二维或三维特征空间里用它。
三是Huber回归虽然稳健，但对“重尾分布”的容忍度有限。如果你的异常值中混有极其离谱的点（比如偏离几十个标准差），Huber的线性区域虽然不会放大它们，但仍会受一定影响。此时可以考虑对目标变量先做截断或分位数变换，再套Huber回归。
四是所有稳健回归都无法替代“方向正确的特征工程”。如果特征本身和目标的真实关系是非线性的，你再怎么换成稳健回归，也只是在错误的形状上画一条顽固的直线。先做多项式扩展、分箱或领域特征构造，然后再考虑用稳健回归来护航。
我个人在实际项目里反复体会最深的一点是：稳健回归不是万能的银弹，但它确实是数据建模中的“保险丝”。当你收到一批来源复杂、质量未知的数据时，先跑一次RANSAC或Huber，比先用OLS拟合再手动清洗要高效得多。把这三行代码放进你的建模工具箱里，以后遇到脏数据，你就不会只剩下“删点”一个选择了。
