1. 从“最小二乘”到“最小截断平方”一个被低估的稳健回归利器在数据分析、金融风控、传感器校准乃至任何涉及从数据中寻找规律的场景里线性回归都是我们工具箱里最基础、最常用的工具。提到线性回归几乎所有人的第一反应就是“最小二乘法”Ordinary Least Squares, OLS。它优雅、计算高效有着完美的数学解释。但从业时间稍长你就会发现一个残酷的现实现实世界的数据远非教科书里那般纯净。一个异常值哪怕只有一个就足以让OLS辛辛苦苦拟合出的那条“最佳”直线偏离十万八千里。这条直线会为了迎合那个离群点而“背叛”绝大多数正常数据点。这种现象我们称之为OLS对异常值“缺乏稳健性”。那么有没有一种方法既能保持线性回归的简洁框架又能像一位经验丰富的侦探一样自动识别并忽略那些试图“带偏节奏”的异常值呢答案是肯定的这就是我们今天要深入探讨的最小截断平方法。请注意这里有一个关键但常被混淆的概念在稳健统计领域LTS通常指的是Least Trimmed Squares即“最小截断平方”而非Least Squares的简单缩写。它是由统计学家Rousseeuw在1984年提出的经典稳健回归方法。而有时在工程或信号处理语境下LTS也可能指“Least Total Squares”等变体但其核心思想与稳健性关联最强的无疑是Least Trimmed Squares。本文我们将聚焦于这个强大的稳健回归器——LTS。你可以把它想象成一个“有主见”的拟合算法。OLS追求的是所有点的误差平方和最小因此每个点都拥有平等的“投票权”异常点的一票影响力巨大。而LTS的思路则非常直接我怀疑数据里混进了“坏蛋”那我就不看全部数据了。我先假设数据中有一部分是“好的”一部分是“坏的”。LTS的目标是找出一部分“最干净”的数据子集比如50%的数据使得用这个子集拟合出来的模型其误差平方和最小。那部分被“截断”掉的数据就被视为潜在的异常值。这个方法简单、粗暴但异常有效尤其适合那些异常值比例可能较高但通常不超过50%的场景比如工业检测中的瑕疵品数据、金融交易中的欺诈数据清洗、环境监测中的传感器瞬时故障等。2. LTS的核心原理如何优雅地“忽略”异常值要理解LTS我们不能只停留在“它很稳健”的感性认知上必须拆解其数学骨架和算法逻辑。这能帮助我们在应用时真正理解每一个参数的意义而不是当一个“调包侠”。2.1 与OLS的对比目标函数的根本差异我们先回顾一下OLS的目标。对于线性模型y Xβ εOLS求解的参数β_ols满足argmin_β Σ_{i1}^{n} (y_i - x_i^T β)^2这里求和是对所有n个数据点进行的。每个残差r_i y_i - x_i^T β都被平方然后求和。平方项放大了大残差的影响这正是OLS对异常值敏感的病根。LTS则完全不同。它首先对所有数据点的残差平方r_i^2进行排序得到有序序列r_(1)}^2 ≤ r_(2)}^2 ≤ ... ≤ r_(n)}^2。然后LTS的目标是argmin_β Σ_{i1}^{h} r_(i)}^2这里h是一个我们预先设定的参数n/2 ≤ h ≤ n。这个目标函数的意思是我只关心残差最小的那h个点也就是拟合得最好的那部分数据并让这h个点的残差平方和最小。至于残差最大的那n-h个点它们的误差直接被“截断”了根本不进入目标函数的计算。这个h参数是LTS的灵魂。它代表了算法所“信任”的数据比例。h越接近nLTS就越像OLSh越接近n/2通常取floor((n1)/2)或floor((np1)/2)其中p是变量数算法就越“挑剔”只采用最核心的那部分数据稳健性也越强。通常我们会将h设置为能覆盖至少50%数据的值这是一个在效率和稳健性之间较好的平衡点。2.2 LTS的算法实现为什么它计算量更大看到这个目标函数一个很自然的问题是怎么求解OLS有解析解β (X^T X)^{-1} X^T y但LTS的目标函数由于引入了排序和截断变得非光滑、非凸没有封闭形式的解。因此LTS的求解依赖于随机抽样算法。最经典的是PROGRESS 算法或其改进版本。其基本思想是一种“采样-验证”的启发式搜索随机子集采样从n个数据点中随机抽取p个点p是自变量个数对于简单线性回归p2。因为p个点可以唯一确定一个超平面两点确定一条直线。初始拟合用这p个点拟合一个初始回归模型计算参数β_init。计算所有残差用β_init计算所有n个数据点的残差r_i。识别“好”的子集根据残差绝对值|r_i|对所有点进行排序选出残差最小的前h个点构成一个“候选好子集”。重新拟合用这个“候选好子集”重新进行OLS拟合得到新的参数β_new并计算这个子集下的目标函数值即这h个点的残差平方和。迭代改进以β_new作为新的起点重复步骤3-5即计算残差、排序、选取前h个点、重新拟合。这个过程可以进行几步C步C-step迭代通常能快速收敛到一个局部最优解。多次随机启动由于初始采样是随机的单次采样可能陷入不好的局部最优。因此我们需要重复上述过程很多次比如500次或1000次在所有这些随机启动中找到那个使目标函数h个点的残差平方和最小的解作为最终的LTS估计。注意正是这个“多次随机采样”的过程使得LTS的计算量远大于OLS。当数据量很大n很大或变量很多p很大时需要的采样次数会急剧增加以保证找到可靠解这是LTS的主要性能瓶颈。在实际应用中我们常使用其快速近似算法如FAST-LTS。2.3 LTS的输出不止一条回归线LTS算法运行完毕后我们会得到一系列宝贵的输出稳健的回归系数基于最优的h个子集拟合出的系数它不受“被截断”的异常值影响。异常值标识算法自然地将数据分为两部分。残差排在最小的h个点之内的被视为“正常点”剩下的n-h个点则被标记为“潜在异常点”。这是一个非常直观的副产品。稳健的尺度估计我们可以基于最优子集的残差计算一个稳健的尺度估计如标准化残差用于后续的统计推断这个尺度估计同样不受异常值污染。3. 手把手实现从理论到Python代码理解了原理我们来看如何动手实现。虽然scikit-learn没有直接提供LTS但我们可以利用statsmodels库或者基于numpy自己实现一个简化版来加深理解。这里我提供一个结合了statsmodels和自定义抽样逻辑的清晰方案。3.1 环境准备与数据构造首先我们创建一个包含明显异常值的仿真数据集。这样我们可以直观对比OLS和LTS的效果。import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.robust.robust_linear_model import RLM import warnings warnings.filterwarnings(ignore) # 忽略一些不影响运行的警告 # 1. 生成干净数据 np.random.seed(42) # 固定随机种子确保结果可复现 n_points 100 x np.linspace(0, 10, n_points) true_slope 2.5 true_intercept 1.0 y_clean true_slope * x true_intercept np.random.normal(0, 1.5, n_points) # 加入高斯噪声 # 2. 人为加入异常值 outlier_indices [20, 40, 60, 80] # 在第20, 40, 60, 80个点位置加入异常 y_contaminated y_clean.copy() y_contaminated[outlier_indices] y_clean[outlier_indices] np.array([25, -20, 30, -25]) # 大幅扰动 # 准备矩阵形式的输入statsmodels需要添加常数项 X sm.add_constant(x) # 添加常数项列对应截距3.2 传统OLS拟合作为对比基线我们用被污染的数据进行普通最小二乘拟合看看异常值的影响有多大。# 使用OLS拟合被污染的数据 model_ols sm.OLS(y_contaminated, X) results_ols model_ols.fit() ols_slope, ols_intercept results_ols.params[1], results_ols.params[0] print(OLS 拟合结果) print(f斜率: {ols_slope:.4f}, 截距: {ols_intercept:.4f}) print(f真实值 - 斜率: {true_slope}, 截距: {true_intercept}) print(*50)运行后你可能会看到OLS的斜率被严重拉偏比如从真实的2.5变成了2.8或3.0截距也变化很大。这就是异常值的破坏力。3.3 利用Statsmodels实现LTS通过M估计逼近statsmodels的RLMRobust Linear Models模块提供了多种稳健回归方法。虽然它没有直接命名为LTS但其RLM配合HuberT()或更稳健的RamsayE()等权重函数并以M估计另一种稳健估计的结果作为迭代重加权最小二乘IRLS的起点可以达到类似LTS的稳健效果。不过为了更贴近LTS的“截断”思想我们可以使用trimmed mean的思路来自定义一个简化版本。下面我展示一个自己实现的、基于随机抽样的简易LTS核心逻辑。请注意这是一个用于教学理解的简化版生产环境应考虑使用更成熟库如R语言的robustbase包或更复杂的算法。def simple_lts(x, y, h_frac0.5, n_trials500): 简易版LTS实现仅适用于简单线性回归 y ax b 参数 x, y: 数据 h_frac: 保留数据的比例默认0.5 n_trials: 随机抽样次数 返回 best_slope, best_intercept: 最优LTS系数 best_indices: 被选为“好子集”的点的索引 all_costs: 每次试验的成本用于调试 n len(x) h int(n * h_frac) # 确保h至少大于等于参数个数简单线性回归参数为2 h max(h, 2) best_cost float(inf) best_slope, best_intercept None, None best_indices None all_costs [] for _ in range(n_trials): # 1. 随机抽取2个点因为两点确定一条直线 sample_idx np.random.choice(n, size2, replaceFalse) x_sample, y_sample x[sample_idx], y[sample_idx] # 2. 用这两个点拟合一条直线 (直接解方程) # 防止两点x坐标相同导致除零 if np.std(x_sample) 1e-10: continue slope_init (y_sample[1] - y_sample[0]) / (x_sample[1] - x_sample[0]) intercept_init y_sample[0] - slope_init * x_sample[0] # 3. C-step迭代这里只做一步完整算法需迭代至收敛 residuals np.abs(y - (slope_init * x intercept_init)) # 获取残差最小的前h个点的索引 good_idx np.argsort(residuals)[:h] x_good, y_good x[good_idx], y[good_idx] # 4. 对好子集进行OLS拟合使用numpy的lstsq A np.vstack([x_good, np.ones(len(x_good))]).T slope_new, intercept_new np.linalg.lstsq(A, y_good, rcondNone)[0] # 5. 计算当前好子集的损失残差平方和 residuals_new y_good - (slope_new * x_good intercept_new) cost np.sum(residuals_new ** 2) all_costs.append(cost) # 6. 更新最优解 if cost best_cost: best_cost cost best_slope, best_intercept slope_new, intercept_new best_indices good_idx.copy() # 注意使用copy() return best_slope, best_intercept, best_indices, all_costs # 应用我们的简易LTS lts_slope, lts_intercept, lts_good_idx, costs simple_lts(x, y_contaminated, h_frac0.55, n_trials1000) print(简易LTS拟合结果) print(f斜率: {lts_slope:.4f}, 截距: {lts_intercept:.4f}) print(f使用的‘好子集’数据点数量: {len(lts_good_idx)})这个简易实现虽然不如专业算法严谨但它清晰地揭示了LTS的核心流程随机抽样、基于残差筛选子集、在子集上做OLS。运行后你会发现LTS拟合出的斜率和截距远比OLS更接近真实值。3.4 结果可视化与对比让我们把三种情况画在一张图上直观感受差异# 计算拟合线 y_ols_fit ols_slope * x ols_intercept y_lts_fit lts_slope * x lts_intercept y_true_fit true_slope * x true_intercept # 绘图 plt.figure(figsize(12, 8)) # 绘制所有数据点 plt.scatter(x, y_contaminated, alpha0.6, label观测数据 (含异常值), colorgray) # 高亮异常值点 plt.scatter(x[outlier_indices], y_contaminated[outlier_indices], colorred, s100, markerx, linewidths3, label人为添加的异常值) # 高亮LTS选出的好子集 plt.scatter(x[lts_good_idx], y_contaminated[lts_good_idx], facecolorsnone, edgecolorsgreen, s80, linewidths2, labelLTS识别的“好子集”) # 绘制三条拟合线 plt.plot(x, y_true_fit, k--, linewidth3, labelf真实关系 (斜率{true_slope})) plt.plot(x, y_ols_fit, r-, linewidth2, labelfOLS拟合 (斜率{ols_slope:.2f})) plt.plot(x, y_lts_fit, b-, linewidth2, labelfLTS拟合 (斜率{lts_slope:.2f})) plt.xlabel(X) plt.ylabel(Y) plt.title(OLS vs LTS 在含异常值数据上的拟合效果对比) plt.legend(locbest) plt.grid(True, alpha0.3) plt.tight_layout() plt.show()通过这张图你可以清晰地看到红线OLS明显被四个红色的异常值“拉拽”整体偏离了真实的黑色虚线。蓝线LTS几乎与黑色虚线重合它成功地忽略了异常值其拟合所依赖的“好子集”绿色圆圈也完美地避开了红色异常点。绿色圆圈直观展示了LTS算法所“信任”的数据部分。4. 关键参数调优与实战避坑指南在实际项目中应用LTS绝不是调用一个函数那么简单。以下几个参数和细节决定了它是“神器”还是“坑器”。4.1 核心参数h的选择艺术与科学的平衡h参数定义了算法认为的“好数据”的最小比例。这是一个需要先验知识或通过经验设置的参数。理论下限h必须至少大于等于模型参数的数量p对于简单线性回归p2否则无法唯一确定模型。通常为了保证估计的稳定性h需要显著大于p。经验值一个常见的经验法则是设置h floor((n p 1)/2)。这保证了至少使用了一半的数据在效率和稳健性之间取得平衡。例如你有100个点简单线性回归p2那么h floor((10021)/2)51。自适应选择在更严谨的应用中h可以通过“稳健性权重图”或“诊断图”来辅助确定。例如你可以尝试多个h值观察回归系数和识别出的异常值是否稳定。如果在一个合理的h范围内比如从50%到80%结果变化不大说明你的模型是稳健的。陷阱h设置过高如90%LTS会退化失去对异常值的抵抗能力。h设置过低如刚好等于p1虽然非常稳健但估计的方差会很大因为用的数据太少且可能把一些好的杠杆点误判为异常值。实操心得我的习惯是首先使用经验公式h floor(0.75 * n)或h floor((n p 1)/2)作为起点。然后如果有领域知识例如我知道这个生产流程的次品率大约在10%我会将h设置为n * (1 - 预期异常比例)并留出一些余量。最后一定会做敏感性分析看看h在±5%范围内波动时关键结论是否改变。4.2 算法稳定性与随机种子LTS依赖于随机抽样这意味着每次运行的结果可能会有细微差异。虽然在大规模抽样下这种差异会很小但在报告结果时为了可重复性固定随机数种子是必须的。在Python中就是np.random.seed(42)。在正式分析报告中应注明所使用的随机种子。4.3 高维数据与“维数灾难”LTS以及许多稳健统计方法在处理高维数据自变量p很大时会面临严峻挑战这被称为“维数灾难”。抽样组合爆炸随机抽样p个点来构造初始解当p很大时抽到“干净”子集的概率极低。为了保证找到全局最优解需要的抽样次数n_trials呈指数级增长计算变得不可行。数据稀疏性在高维空间所有数据点都可能显得很“远”异常值的概念变得模糊基于距离残差的截断方法效果下降。应对策略变量选择/降维在应用LTS之前先使用主成分分析PCA或变量选择方法降低维度。使用专门的高维稳健方法考虑稀疏LTSSparse LTS它结合了L1正则化LASSO进行变量选择和稳健估计。分步处理先使用更快的异常检测方法如基于马氏距离的方法初步筛选再在“相对干净”的子集上应用LTS。4.4 异常值诊断不要盲目相信LTS的判决LTS给出了一个异常值列表但这只是一个统计判断绝不能不加思考地直接删除。业务复核必须结合业务逻辑检查这些被标记的点。它真的是数据录入错误、传感器故障吗还是它代表了一种罕见但真实的业务模式比如一笔巨额合法交易盲目删除后者会导致模型丢失重要信息。杠杆点 vs. 离群点LTS主要抵抗y方向的离群点。对于x方向的异常值高杠杆点即使它的y值很正常也可能对OLS估计产生巨大影响。LTS对高杠杆点的处理能力相对有限。需要结合库克距离等指标综合判断。建议流程将LTS标记的异常点作为“可疑名单”交给业务专家审核。确认是错误后可以选择删除、修正或者使用稳健回归的结果而不删除数据因为LTS的系数本身已不受其影响。5. 超越基础LTS的变体与在复杂场景中的应用掌握了标准LTS我们可以看看它的几个变体以及如何将它应用到更复杂的现实问题中。5.1 稀疏LTS当变量很多时前面提到高维问题稀疏LTSSparse Least Trimmed Squares就是一个优雅的解决方案。它在LTS的目标函数中加入了L1正则化项argmin_β Σ_{i1}^{h} r_(i)}^2 λ * ||β||_1其中||β||_1是系数向量的L1范数绝对值之和λ是正则化强度参数。这样做的效果是同时进行稳健回归和变量选择。在寻找最优h子集的同时它也会将一些不重要的变量的系数压缩至0。这在基因数据、金融因子模型等变量成千上万的场景中非常有用。在R语言的robustbase包中就有sparseLTS函数。在Python生态中你可以通过sklearn_extra库中的RANSACRegressor配合某些设置来近似实现类似思想或者自己实现坐标下降算法进行求解。5.2 加权LTS与迭代重加权最小二乘有时我们对不同数据点的信任程度不同。加权LTS允许我们为每个数据点赋予一个先验权重w_i。其目标函数变为argmin_β Σ_{i1}^{h} (w_(i) * r_(i)}^2)这里权重w_i可以来自业务知识例如某些传感器的精度更高也可以来自另一轮稳健估计的结果。这引出了迭代重加权最小二乘的思想先进行一次LTS拟合根据残差计算每个点的权重残差大的点权重小然后用加权最小二乘进行下一次拟合如此迭代直至收敛。statsmodels的RLM本质上就是这种思想。5.3 在时间序列与金融数据中的应用案例金融数据中充满了“尖峰厚尾”和结构性断点OLS回归在这里常常失灵。案例股票Beta系数估计。股票的Beta系数通常通过其收益率对市场收益率做回归得到。但市场暴跌或暴涨时异常值OLS估计的Beta会失真。使用LTS回归可以剔除这些极端日子的影响得到一个更稳健、更能代表“正常市场状况”下的Beta值。这对于风险管理和资产定价至关重要。实操步骤获取标的股票和基准指数如沪深300的日收益率序列。以指数收益率为X股票收益率为Y构建回归模型Y α β * X ε。分别用OLS和LTS设置h约为75%进行拟合。对比两个Beta值。通常会发现在波动剧烈的时期OLS的Beta值会显著偏离LTS的Beta值。LTS给出的Beta更稳定受单日极端行情影响小。关键检查查看LTS识别出的异常值日期往往对应着公司特殊公告日、市场黑天鹅事件日等这反过来也验证了模型的业务直觉。5.4 与RANSAC的对比另一种稳健回归思路在计算机视觉和机器学习中另一个常用的稳健回归方法是RANSAC。它和LTS哲学相似但实现不同LTS目标是找到一个固定大小h的“最好”子集。RANSAC随机采样一个最小子集如2个点拟合模型然后计算有多少点符合这个模型即残差小于某个阈值t这些点称为“内点”。重复多次选择“内点”最多的那个模型最后用所有“内点”重新拟合最终模型。对比效率RANSAC在异常值比例很高时可能更高效因为它不需要对全部残差排序只需要计数。参数RANSAC需要设定阈值t来定义“内点”这个阈值有时比LTS的h更难以确定。适用性RANSAC更通用可用于任何可以用参数模型描述的问题如拟合圆、单应性矩阵。LTS是专门为回归问题设计的。结果两者目标类似结果通常接近。在实践中如果数据量不大可以两者都尝试相互印证。我个人在处理纯粹的线性回归稳健拟合问题时更倾向于使用LTS因为它的统计解释更清晰最小化截断平方和。而在处理图像配准、点云匹配等更复杂的模型拟合时RANSAC是首选工具。 
长河君 · 2026/8/27 14:04:51 · 阅读 1.2k 最小截断平方回归：从原理到Python实现，解决异常值干扰的稳健回归方法
1. 从“最小二乘”到“最小截断平方”一个被低估的稳健回归利器在数据分析、金融风控、传感器校准乃至任何涉及从数据中寻找规律的场景里线性回归都是我们工具箱里最基础、最常用的工具。提到线性回归几乎所有人的第一反应就是“最小二乘法”Ordinary Least Squares, OLS。它优雅、计算高效有着完美的数学解释。但从业时间稍长你就会发现一个残酷的现实现实世界的数据远非教科书里那般纯净。一个异常值哪怕只有一个就足以让OLS辛辛苦苦拟合出的那条“最佳”直线偏离十万八千里。这条直线会为了迎合那个离群点而“背叛”绝大多数正常数据点。这种现象我们称之为OLS对异常值“缺乏稳健性”。那么有没有一种方法既能保持线性回归的简洁框架又能像一位经验丰富的侦探一样自动识别并忽略那些试图“带偏节奏”的异常值呢答案是肯定的这就是我们今天要深入探讨的最小截断平方法。请注意这里有一个关键但常被混淆的概念在稳健统计领域LTS通常指的是Least Trimmed Squares即“最小截断平方”而非Least Squares的简单缩写。它是由统计学家Rousseeuw在1984年提出的经典稳健回归方法。而有时在工程或信号处理语境下LTS也可能指“Least Total Squares”等变体但其核心思想与稳健性关联最强的无疑是Least Trimmed Squares。本文我们将聚焦于这个强大的稳健回归器——LTS。你可以把它想象成一个“有主见”的拟合算法。OLS追求的是所有点的误差平方和最小因此每个点都拥有平等的“投票权”异常点的一票影响力巨大。而LTS的思路则非常直接我怀疑数据里混进了“坏蛋”那我就不看全部数据了。我先假设数据中有一部分是“好的”一部分是“坏的”。LTS的目标是找出一部分“最干净”的数据子集比如50%的数据使得用这个子集拟合出来的模型其误差平方和最小。那部分被“截断”掉的数据就被视为潜在的异常值。这个方法简单、粗暴但异常有效尤其适合那些异常值比例可能较高但通常不超过50%的场景比如工业检测中的瑕疵品数据、金融交易中的欺诈数据清洗、环境监测中的传感器瞬时故障等。2. LTS的核心原理如何优雅地“忽略”异常值要理解LTS我们不能只停留在“它很稳健”的感性认知上必须拆解其数学骨架和算法逻辑。这能帮助我们在应用时真正理解每一个参数的意义而不是当一个“调包侠”。2.1 与OLS的对比目标函数的根本差异我们先回顾一下OLS的目标。对于线性模型y Xβ εOLS求解的参数β_ols满足argmin_β Σ_{i1}^{n} (y_i - x_i^T β)^2这里求和是对所有n个数据点进行的。每个残差r_i y_i - x_i^T β都被平方然后求和。平方项放大了大残差的影响这正是OLS对异常值敏感的病根。LTS则完全不同。它首先对所有数据点的残差平方r_i^2进行排序得到有序序列r_(1)}^2 ≤ r_(2)}^2 ≤ ... ≤ r_(n)}^2。然后LTS的目标是argmin_β Σ_{i1}^{h} r_(i)}^2这里h是一个我们预先设定的参数n/2 ≤ h ≤ n。这个目标函数的意思是我只关心残差最小的那h个点也就是拟合得最好的那部分数据并让这h个点的残差平方和最小。至于残差最大的那n-h个点它们的误差直接被“截断”了根本不进入目标函数的计算。这个h参数是LTS的灵魂。它代表了算法所“信任”的数据比例。h越接近nLTS就越像OLSh越接近n/2通常取floor((n1)/2)或floor((np1)/2)其中p是变量数算法就越“挑剔”只采用最核心的那部分数据稳健性也越强。通常我们会将h设置为能覆盖至少50%数据的值这是一个在效率和稳健性之间较好的平衡点。2.2 LTS的算法实现为什么它计算量更大看到这个目标函数一个很自然的问题是怎么求解OLS有解析解β (X^T X)^{-1} X^T y但LTS的目标函数由于引入了排序和截断变得非光滑、非凸没有封闭形式的解。因此LTS的求解依赖于随机抽样算法。最经典的是PROGRESS 算法或其改进版本。其基本思想是一种“采样-验证”的启发式搜索随机子集采样从n个数据点中随机抽取p个点p是自变量个数对于简单线性回归p2。因为p个点可以唯一确定一个超平面两点确定一条直线。初始拟合用这p个点拟合一个初始回归模型计算参数β_init。计算所有残差用β_init计算所有n个数据点的残差r_i。识别“好”的子集根据残差绝对值|r_i|对所有点进行排序选出残差最小的前h个点构成一个“候选好子集”。重新拟合用这个“候选好子集”重新进行OLS拟合得到新的参数β_new并计算这个子集下的目标函数值即这h个点的残差平方和。迭代改进以β_new作为新的起点重复步骤3-5即计算残差、排序、选取前h个点、重新拟合。这个过程可以进行几步C步C-step迭代通常能快速收敛到一个局部最优解。多次随机启动由于初始采样是随机的单次采样可能陷入不好的局部最优。因此我们需要重复上述过程很多次比如500次或1000次在所有这些随机启动中找到那个使目标函数h个点的残差平方和最小的解作为最终的LTS估计。注意正是这个“多次随机采样”的过程使得LTS的计算量远大于OLS。当数据量很大n很大或变量很多p很大时需要的采样次数会急剧增加以保证找到可靠解这是LTS的主要性能瓶颈。在实际应用中我们常使用其快速近似算法如FAST-LTS。2.3 LTS的输出不止一条回归线LTS算法运行完毕后我们会得到一系列宝贵的输出稳健的回归系数基于最优的h个子集拟合出的系数它不受“被截断”的异常值影响。异常值标识算法自然地将数据分为两部分。残差排在最小的h个点之内的被视为“正常点”剩下的n-h个点则被标记为“潜在异常点”。这是一个非常直观的副产品。稳健的尺度估计我们可以基于最优子集的残差计算一个稳健的尺度估计如标准化残差用于后续的统计推断这个尺度估计同样不受异常值污染。3. 手把手实现从理论到Python代码理解了原理我们来看如何动手实现。虽然scikit-learn没有直接提供LTS但我们可以利用statsmodels库或者基于numpy自己实现一个简化版来加深理解。这里我提供一个结合了statsmodels和自定义抽样逻辑的清晰方案。3.1 环境准备与数据构造首先我们创建一个包含明显异常值的仿真数据集。这样我们可以直观对比OLS和LTS的效果。import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.robust.robust_linear_model import RLM import warnings warnings.filterwarnings(ignore) # 忽略一些不影响运行的警告 # 1. 生成干净数据 np.random.seed(42) # 固定随机种子确保结果可复现 n_points 100 x np.linspace(0, 10, n_points) true_slope 2.5 true_intercept 1.0 y_clean true_slope * x true_intercept np.random.normal(0, 1.5, n_points) # 加入高斯噪声 # 2. 人为加入异常值 outlier_indices [20, 40, 60, 80] # 在第20, 40, 60, 80个点位置加入异常 y_contaminated y_clean.copy() y_contaminated[outlier_indices] y_clean[outlier_indices] np.array([25, -20, 30, -25]) # 大幅扰动 # 准备矩阵形式的输入statsmodels需要添加常数项 X sm.add_constant(x) # 添加常数项列对应截距3.2 传统OLS拟合作为对比基线我们用被污染的数据进行普通最小二乘拟合看看异常值的影响有多大。# 使用OLS拟合被污染的数据 model_ols sm.OLS(y_contaminated, X) results_ols model_ols.fit() ols_slope, ols_intercept results_ols.params[1], results_ols.params[0] print(OLS 拟合结果) print(f斜率: {ols_slope:.4f}, 截距: {ols_intercept:.4f}) print(f真实值 - 斜率: {true_slope}, 截距: {true_intercept}) print(*50)运行后你可能会看到OLS的斜率被严重拉偏比如从真实的2.5变成了2.8或3.0截距也变化很大。这就是异常值的破坏力。3.3 利用Statsmodels实现LTS通过M估计逼近statsmodels的RLMRobust Linear Models模块提供了多种稳健回归方法。虽然它没有直接命名为LTS但其RLM配合HuberT()或更稳健的RamsayE()等权重函数并以M估计另一种稳健估计的结果作为迭代重加权最小二乘IRLS的起点可以达到类似LTS的稳健效果。不过为了更贴近LTS的“截断”思想我们可以使用trimmed mean的思路来自定义一个简化版本。下面我展示一个自己实现的、基于随机抽样的简易LTS核心逻辑。请注意这是一个用于教学理解的简化版生产环境应考虑使用更成熟库如R语言的robustbase包或更复杂的算法。def simple_lts(x, y, h_frac0.5, n_trials500): 简易版LTS实现仅适用于简单线性回归 y ax b 参数 x, y: 数据 h_frac: 保留数据的比例默认0.5 n_trials: 随机抽样次数 返回 best_slope, best_intercept: 最优LTS系数 best_indices: 被选为“好子集”的点的索引 all_costs: 每次试验的成本用于调试 n len(x) h int(n * h_frac) # 确保h至少大于等于参数个数简单线性回归参数为2 h max(h, 2) best_cost float(inf) best_slope, best_intercept None, None best_indices None all_costs [] for _ in range(n_trials): # 1. 随机抽取2个点因为两点确定一条直线 sample_idx np.random.choice(n, size2, replaceFalse) x_sample, y_sample x[sample_idx], y[sample_idx] # 2. 用这两个点拟合一条直线 (直接解方程) # 防止两点x坐标相同导致除零 if np.std(x_sample) 1e-10: continue slope_init (y_sample[1] - y_sample[0]) / (x_sample[1] - x_sample[0]) intercept_init y_sample[0] - slope_init * x_sample[0] # 3. C-step迭代这里只做一步完整算法需迭代至收敛 residuals np.abs(y - (slope_init * x intercept_init)) # 获取残差最小的前h个点的索引 good_idx np.argsort(residuals)[:h] x_good, y_good x[good_idx], y[good_idx] # 4. 对好子集进行OLS拟合使用numpy的lstsq A np.vstack([x_good, np.ones(len(x_good))]).T slope_new, intercept_new np.linalg.lstsq(A, y_good, rcondNone)[0] # 5. 计算当前好子集的损失残差平方和 residuals_new y_good - (slope_new * x_good intercept_new) cost np.sum(residuals_new ** 2) all_costs.append(cost) # 6. 更新最优解 if cost best_cost: best_cost cost best_slope, best_intercept slope_new, intercept_new best_indices good_idx.copy() # 注意使用copy() return best_slope, best_intercept, best_indices, all_costs # 应用我们的简易LTS lts_slope, lts_intercept, lts_good_idx, costs simple_lts(x, y_contaminated, h_frac0.55, n_trials1000) print(简易LTS拟合结果) print(f斜率: {lts_slope:.4f}, 截距: {lts_intercept:.4f}) print(f使用的‘好子集’数据点数量: {len(lts_good_idx)})这个简易实现虽然不如专业算法严谨但它清晰地揭示了LTS的核心流程随机抽样、基于残差筛选子集、在子集上做OLS。运行后你会发现LTS拟合出的斜率和截距远比OLS更接近真实值。3.4 结果可视化与对比让我们把三种情况画在一张图上直观感受差异# 计算拟合线 y_ols_fit ols_slope * x ols_intercept y_lts_fit lts_slope * x lts_intercept y_true_fit true_slope * x true_intercept # 绘图 plt.figure(figsize(12, 8)) # 绘制所有数据点 plt.scatter(x, y_contaminated, alpha0.6, label观测数据 (含异常值), colorgray) # 高亮异常值点 plt.scatter(x[outlier_indices], y_contaminated[outlier_indices], colorred, s100, markerx, linewidths3, label人为添加的异常值) # 高亮LTS选出的好子集 plt.scatter(x[lts_good_idx], y_contaminated[lts_good_idx], facecolorsnone, edgecolorsgreen, s80, linewidths2, labelLTS识别的“好子集”) # 绘制三条拟合线 plt.plot(x, y_true_fit, k--, linewidth3, labelf真实关系 (斜率{true_slope})) plt.plot(x, y_ols_fit, r-, linewidth2, labelfOLS拟合 (斜率{ols_slope:.2f})) plt.plot(x, y_lts_fit, b-, linewidth2, labelfLTS拟合 (斜率{lts_slope:.2f})) plt.xlabel(X) plt.ylabel(Y) plt.title(OLS vs LTS 在含异常值数据上的拟合效果对比) plt.legend(locbest) plt.grid(True, alpha0.3) plt.tight_layout() plt.show()通过这张图你可以清晰地看到红线OLS明显被四个红色的异常值“拉拽”整体偏离了真实的黑色虚线。蓝线LTS几乎与黑色虚线重合它成功地忽略了异常值其拟合所依赖的“好子集”绿色圆圈也完美地避开了红色异常点。绿色圆圈直观展示了LTS算法所“信任”的数据部分。4. 关键参数调优与实战避坑指南在实际项目中应用LTS绝不是调用一个函数那么简单。以下几个参数和细节决定了它是“神器”还是“坑器”。4.1 核心参数h的选择艺术与科学的平衡h参数定义了算法认为的“好数据”的最小比例。这是一个需要先验知识或通过经验设置的参数。理论下限h必须至少大于等于模型参数的数量p对于简单线性回归p2否则无法唯一确定模型。通常为了保证估计的稳定性h需要显著大于p。经验值一个常见的经验法则是设置h floor((n p 1)/2)。这保证了至少使用了一半的数据在效率和稳健性之间取得平衡。例如你有100个点简单线性回归p2那么h floor((10021)/2)51。自适应选择在更严谨的应用中h可以通过“稳健性权重图”或“诊断图”来辅助确定。例如你可以尝试多个h值观察回归系数和识别出的异常值是否稳定。如果在一个合理的h范围内比如从50%到80%结果变化不大说明你的模型是稳健的。陷阱h设置过高如90%LTS会退化失去对异常值的抵抗能力。h设置过低如刚好等于p1虽然非常稳健但估计的方差会很大因为用的数据太少且可能把一些好的杠杆点误判为异常值。实操心得我的习惯是首先使用经验公式h floor(0.75 * n)或h floor((n p 1)/2)作为起点。然后如果有领域知识例如我知道这个生产流程的次品率大约在10%我会将h设置为n * (1 - 预期异常比例)并留出一些余量。最后一定会做敏感性分析看看h在±5%范围内波动时关键结论是否改变。4.2 算法稳定性与随机种子LTS依赖于随机抽样这意味着每次运行的结果可能会有细微差异。虽然在大规模抽样下这种差异会很小但在报告结果时为了可重复性固定随机数种子是必须的。在Python中就是np.random.seed(42)。在正式分析报告中应注明所使用的随机种子。4.3 高维数据与“维数灾难”LTS以及许多稳健统计方法在处理高维数据自变量p很大时会面临严峻挑战这被称为“维数灾难”。抽样组合爆炸随机抽样p个点来构造初始解当p很大时抽到“干净”子集的概率极低。为了保证找到全局最优解需要的抽样次数n_trials呈指数级增长计算变得不可行。数据稀疏性在高维空间所有数据点都可能显得很“远”异常值的概念变得模糊基于距离残差的截断方法效果下降。应对策略变量选择/降维在应用LTS之前先使用主成分分析PCA或变量选择方法降低维度。使用专门的高维稳健方法考虑稀疏LTSSparse LTS它结合了L1正则化LASSO进行变量选择和稳健估计。分步处理先使用更快的异常检测方法如基于马氏距离的方法初步筛选再在“相对干净”的子集上应用LTS。4.4 异常值诊断不要盲目相信LTS的判决LTS给出了一个异常值列表但这只是一个统计判断绝不能不加思考地直接删除。业务复核必须结合业务逻辑检查这些被标记的点。它真的是数据录入错误、传感器故障吗还是它代表了一种罕见但真实的业务模式比如一笔巨额合法交易盲目删除后者会导致模型丢失重要信息。杠杆点 vs. 离群点LTS主要抵抗y方向的离群点。对于x方向的异常值高杠杆点即使它的y值很正常也可能对OLS估计产生巨大影响。LTS对高杠杆点的处理能力相对有限。需要结合库克距离等指标综合判断。建议流程将LTS标记的异常点作为“可疑名单”交给业务专家审核。确认是错误后可以选择删除、修正或者使用稳健回归的结果而不删除数据因为LTS的系数本身已不受其影响。5. 超越基础LTS的变体与在复杂场景中的应用掌握了标准LTS我们可以看看它的几个变体以及如何将它应用到更复杂的现实问题中。5.1 稀疏LTS当变量很多时前面提到高维问题稀疏LTSSparse Least Trimmed Squares就是一个优雅的解决方案。它在LTS的目标函数中加入了L1正则化项argmin_β Σ_{i1}^{h} r_(i)}^2 λ * ||β||_1其中||β||_1是系数向量的L1范数绝对值之和λ是正则化强度参数。这样做的效果是同时进行稳健回归和变量选择。在寻找最优h子集的同时它也会将一些不重要的变量的系数压缩至0。这在基因数据、金融因子模型等变量成千上万的场景中非常有用。在R语言的robustbase包中就有sparseLTS函数。在Python生态中你可以通过sklearn_extra库中的RANSACRegressor配合某些设置来近似实现类似思想或者自己实现坐标下降算法进行求解。5.2 加权LTS与迭代重加权最小二乘有时我们对不同数据点的信任程度不同。加权LTS允许我们为每个数据点赋予一个先验权重w_i。其目标函数变为argmin_β Σ_{i1}^{h} (w_(i) * r_(i)}^2)这里权重w_i可以来自业务知识例如某些传感器的精度更高也可以来自另一轮稳健估计的结果。这引出了迭代重加权最小二乘的思想先进行一次LTS拟合根据残差计算每个点的权重残差大的点权重小然后用加权最小二乘进行下一次拟合如此迭代直至收敛。statsmodels的RLM本质上就是这种思想。5.3 在时间序列与金融数据中的应用案例金融数据中充满了“尖峰厚尾”和结构性断点OLS回归在这里常常失灵。案例股票Beta系数估计。股票的Beta系数通常通过其收益率对市场收益率做回归得到。但市场暴跌或暴涨时异常值OLS估计的Beta会失真。使用LTS回归可以剔除这些极端日子的影响得到一个更稳健、更能代表“正常市场状况”下的Beta值。这对于风险管理和资产定价至关重要。实操步骤获取标的股票和基准指数如沪深300的日收益率序列。以指数收益率为X股票收益率为Y构建回归模型Y α β * X ε。分别用OLS和LTS设置h约为75%进行拟合。对比两个Beta值。通常会发现在波动剧烈的时期OLS的Beta值会显著偏离LTS的Beta值。LTS给出的Beta更稳定受单日极端行情影响小。关键检查查看LTS识别出的异常值日期往往对应着公司特殊公告日、市场黑天鹅事件日等这反过来也验证了模型的业务直觉。5.4 与RANSAC的对比另一种稳健回归思路在计算机视觉和机器学习中另一个常用的稳健回归方法是RANSAC。它和LTS哲学相似但实现不同LTS目标是找到一个固定大小h的“最好”子集。RANSAC随机采样一个最小子集如2个点拟合模型然后计算有多少点符合这个模型即残差小于某个阈值t这些点称为“内点”。重复多次选择“内点”最多的那个模型最后用所有“内点”重新拟合最终模型。对比效率RANSAC在异常值比例很高时可能更高效因为它不需要对全部残差排序只需要计数。参数RANSAC需要设定阈值t来定义“内点”这个阈值有时比LTS的h更难以确定。适用性RANSAC更通用可用于任何可以用参数模型描述的问题如拟合圆、单应性矩阵。LTS是专门为回归问题设计的。结果两者目标类似结果通常接近。在实践中如果数据量不大可以两者都尝试相互印证。我个人在处理纯粹的线性回归稳健拟合问题时更倾向于使用LTS因为它的统计解释更清晰最小化截断平方和。而在处理图像配准、点云匹配等更复杂的模型拟合时RANSAC是首选工具。