最小截平方和法(LTS)：高崩溃点稳健回归原理与Python实现

1. 项目概述从“拟合”到“稳健”的回归进化在数据分析和数学建模的世界里回归分析是当之无愧的基石。无论是预测房价、分析广告点击率还是研究药物剂量与疗效的关系我们都在试图用一个或多个变量自变量去解释另一个变量因变量的变化规律。最经典、最广为人知的方法莫过于最小二乘法。它通过最小化预测值与真实值之间误差的平方和找到一条“最佳”的拟合直线或曲线。然而但凡在实际项目中用过经典最小二乘法的朋友大概率都踩过同一个坑它对异常值过于敏感。几个偏离主流的“捣蛋鬼”数据点就能把整条回归线拽得偏离正轨导致模型预测失准。这就像用一把普通的尺子去量一块边缘有毛刺的木板几个毛刺点就能让测量结果谬以千里。为了解决这个问题统计学家们提出了各种稳健回归方法。今天我们要深入探讨的最小截平方和法就是其中一种极具代表性的“抗异常值”利器。LTS不像经典最小二乘法那样试图讨好所有数据点而是“聪明”地寻找数据中的“大多数”共识。它只利用一部分最“干净”、最“一致”的数据子集来拟合模型从而有效抵御异常值的干扰。理解并掌握LTS意味着你在处理真实世界往往充满噪声和异常值的数据时手中多了一件可靠的工具。无论你是数学建模竞赛的选手还是从事数据分析、金融风控、工业质量控制的工程师这项技术都能让你构建的模型更加稳健、可信。2. 核心原理为什么LTS比OLS更“抗造”要理解LTS我们必须先回顾一下它的“前辈”——普通最小二乘法的软肋。2.1 经典最小二乘法的阿喀琉斯之踵普通最小二乘法的目标函数是最小化所有残差观测值减去预测值的平方和。用公式表示就是Minimize Σ(y_i - ŷ_i)²其中y_i是真实值ŷ_i是模型预测值。这个方法的优点是数学性质优美有解析解计算高效。但其致命弱点就藏在“平方”二字里。平方操作会放大较大误差的影响。假设有一个异常点其残差为10那么它对目标函数的“贡献”是10²100。而一个正常点残差为1贡献仅为1。这意味着为了最小化总平方和模型会不惜扭曲整体趋势去“迎合”那个贡献了100的异常点哪怕牺牲99个正常点的拟合精度也在所不惜。这种现象在统计上被称为缺乏崩溃点即极少量的严重异常值就足以让估计结果完全失效。2.2 LTS的稳健哲学寻找“干净的多数”LTS采用了完全不同的策略。它的核心思想是从全部n个数据点中选出所有可能的、包含h个点的子集h通常大于n/2例如取n的75%。对每一个这样的子集用普通最小二乘法拟合一个模型并计算这个模型在整个数据集上的残差平方和但只取其中最小的h个残差平方进行求和。最终选择那个使得“最小h个残差平方和”最小的子集所对应的模型参数作为LTS的估计结果。简单来说LTS不是在找一条让“所有点”的误差平方和最小的线而是在找一条让“大多数好点”的误差平方和最小的线。它主动忽略掉那些误差最大的“坏点”异常值。参数h控制了模型的“稳健性”和“效率”之间的权衡h越大越接近n模型效率越高估计更精确但稳健性越差h越小模型越稳健但效率会损失。通常h的取值在floor((np1)/2)到n之间其中p是自变量个数。一个常用的经验值是h floor(0.75 * n)即使用75%的数据。注意LTS的“最小截平方和”中的“截”指的就是截取一部分h个残差最小的点。这使它拥有了高达50%的崩溃点即即使有接近一半的数据是异常值LTS依然能给出相对合理的估计。这是经典最小二乘法无法企及的优势。2.3 与相关热词的辨析LTS不是“长期支持版”在展开具体实现前有必要澄清一个由热词搜索带来的常见误解。网络热词中出现了大量如“Ubuntu 22.04 LTS”、“Photoshop 2025 LTS”等。这里的LTS是“Long-Term Support”的缩写意为“长期支持”是软件版本的一种发布策略与我们的统计方法最小截平方和法风马牛不相及。同样热词中的“随机森林回归算法”是一种基于决策树集成的强大机器学习方法虽然也用于回归且对异常值有一定鲁棒性但其原理Bagging随机特征与LTS这种基于子集搜索的稳健统计方法截然不同。我们在学习和交流时务必注意区分上下文。3. 算法实现与核心步骤拆解理解了原理我们来看如何将LTS从数学公式变为可运行的代码。由于其核心是组合优化问题从n个点中选h个暴力枚举所有子集在数据量稍大时就是计算灾难C(n, h)增长极快。因此实际中普遍采用随机抽样迭代算法来逼近最优解。3.1 算法流程详解一个典型的LTS算法实现包含以下步骤初始化与参数设置输入数据矩阵X包含常数项如果有的话响应向量y以及参数h子集大小。设置最大迭代次数max_iter如500次和随机子集采样数量n_subsets如1000个。这些参数是为了在计算时间和求解精度间取得平衡。随机子集采样与初步拟合循环进行n_subsets次 a. 从全部n个观测中完全随机地抽取一个大小为p1的子集p是自变量个数p1是能唯一确定一个线性模型的最小样本数。 b. 用这个小子集的数据通过普通最小二乘法计算出一组初始回归系数β_init。 c. 用这组β_init去预测所有n个点的值计算所有n个残差r_i y_i - X_iβ_init。 d. 对这n个残差取绝对值并排序找出残差绝对值最小的前h个观测点构成一个新的“候选干净子集”。基于候选子集的精确拟合与迭代改进对上一步得到的“候选干净子集”大小为h使用普通最小二乘法进行拟合得到一组新的系数β_candidate。计算使用β_candidate时所有n个点的残差并同样取绝对值最小的前h个残差平方和记为SSR_h。记录下到目前为止得到的最小的SSR_h及其对应的系数β_best。C-Step迭代收敛对于表现较好的候选解例如SSR_h较小的前10%可以进行称为“C-Step”的迭代改进 a. 给定当前系数β计算全样本残差选出残差绝对值最小的h个点构成新子集。 b. 用新子集重新拟合OLS得到新的β。 c. 重复a和b直到子集成员不再变化或SSR_h不再减小。这个过程能快速将随机得到的初始解“抛光”到局部最优。输出最终结果在所有随机采样和C-Step改进完成后选择那个使得SSR_h最小的系数β_best作为LTS的最终估计。同时算法可以输出被识别为异常值的点即最终未进入最优h子集的那些观测。3.2 关键参数选择与调优经验h子集大小这是LTS最重要的参数。我的经验是在没有任何先验信息时从h floor(0.75 * n)开始尝试是一个稳健的起点。如果你对数据的污染程度有一个粗略估计例如认为最多有20%的异常值那么可以设置h floor((1 - 污染比例) * n)。在实践中可以尝试几个不同的h值如0.7, 0.75, 0.8观察模型系数和异常值识别结果的稳定性。n_subsets随机子集数理论上采样越多找到全局最优解的概率越大。但计算成本也线性增加。对于中小规模数据n 10001000-5000次采样通常足够。对于大规模数据可能需要根据时间预算进行调整。一个技巧是可以先用一个较小的n_subsets如500快速运行如果结果不稳定多次运行结果差异大再增加采样次数。随机种子由于算法包含随机采样为了结果可复现务必固定随机数生成器的种子。3.3 一个清晰的Python实现示例下面我们抛开复杂的统计包用NumPy从头实现一个简化版的LTS以便彻底理解其每一步。我们将使用一个包含明显异常值的合成数据集来演示。import numpy as np import matplotlib.pyplot as plt def least_trimmed_squares(X, y, h, n_subsets1000, max_csteps10, random_state42): 最小截平方和法LTS的简单实现。 参数: X : numpy array, 形状 (n_samples, n_features) 包含常数项如有。 y : numpy array, 形状 (n_samples,) h : int, 用于拟合的子集大小 n/2 h n。 n_subsets : int, 随机初始子集的数量。 max_csteps : int, C-Step迭代的最大次数。 random_state : int, 随机种子。 返回: beta_best : numpy array, 最优的回归系数。 inlier_mask : boolean array, 形状 (n_samples,) True表示内点在最优h子集中。 best_ssr : float, 最优的截断残差平方和。 np.random.seed(random_state) n_samples, n_features X.shape # 参数校验 if not (n_samples//2 h n_samples): raise ValueError(fh must be between n/2 and n. Got h{h}, n{n_samples}) beta_best None best_ssr np.inf best_inlier_mask None # 1. 随机采样多个初始子集 for _ in range(n_subsets): # 随机选择 p1 个点作为初始子集 (确保非奇异) random_indices np.random.choice(n_samples, sizen_features, replaceFalse) X_sub X[random_indices] y_sub y[random_indices] # 2. 初始OLS拟合 (使用伪逆避免奇异矩阵问题) try: beta_init np.linalg.lstsq(X_sub, y_sub, rcondNone)[0] except np.linalg.LinAlgError: continue # 如果初始子阵奇异跳过此次迭代 beta_current beta_init.copy() # 3. C-Step 迭代改进 for _ in range(max_csteps): # 计算所有残差 residuals y - X.dot(beta_current) abs_residuals np.abs(residuals) # 找到残差绝对值最小的 h 个点的索引 inlier_indices np.argpartition(abs_residuals, h-1)[:h] # 用这 h 个点重新拟合 OLS X_h X[inlier_indices] y_h y[inlier_indices] beta_new np.linalg.lstsq(X_h, y_h, rcondNone)[0] # 检查收敛系数是否基本不变 if np.linalg.norm(beta_new - beta_current) 1e-8: beta_current beta_new break beta_current beta_new # 计算当前解对应的截断残差平方和 residuals_final y - X.dot(beta_current) squared_residuals residuals_final ** 2 # 取最小的 h 个残差平方和 ssr_h np.sum(np.sort(squared_residuals)[:h]) # 4. 更新最优解 if ssr_h best_ssr: best_ssr ssr_h beta_best beta_current.copy() # 确定最终的内点掩码 abs_residuals_final np.abs(residuals_final) best_inlier_mask np.argsort(abs_residuals_final) h return beta_best, best_inlier_mask, best_ssr # --- 生成示例数据包含异常值--- np.random.seed(0) n 100 X np.linspace(0, 10, n) # 真实关系y 2*X 5 噪声 y_true 2 * X 5 noise np.random.randn(n) * 1.5 y_clean y_true noise # 故意添加几个异常值 outlier_indices [20, 40, 60, 80] y y_clean.copy() y[outlier_indices] np.array([25, -20, 30, -25]) # 大幅扰动 # 准备设计矩阵添加常数项 X_design np.column_stack([np.ones_like(X), X]) # --- 应用经典OLS和我们的LTS --- # 经典OLS (对所有数据) beta_ols np.linalg.lstsq(X_design, y, rcondNone)[0] y_pred_ols X_design.dot(beta_ols) # LTS (假设我们认为有约20%的异常值h取80) h int(0.8 * n) beta_lts, inlier_mask, ssr_lts least_trimmed_squares(X_design, y, h, n_subsets500) y_pred_lts X_design.dot(beta_lts) # --- 可视化对比 --- plt.figure(figsize(12, 6)) plt.scatter(X, y, alpha0.6, label数据点 (含异常值), cgray) plt.scatter(X[inlier_mask], y[inlier_mask], alpha0.8, labelLTS识别的内点, cgreen) plt.scatter(X[outlier_indices], y[outlier_indices], markerx, s100, label真实异常值位置, cred, linewidths2) x_plot np.linspace(0, 10, 100) X_plot_design np.column_stack([np.ones_like(x_plot), x_plot]) plt.plot(x_plot, X_plot_design.dot(beta_ols), r--, linewidth2, labelf经典OLS: y{beta_ols[1]:.2f}x{beta_ols[0]:.2f}) plt.plot(x_plot, X_plot_design.dot(beta_lts), b-, linewidth2, labelfLTS (h{h}): y{beta_lts[1]:.2f}x{beta_lts[0]:.2f}) plt.plot(x_plot, 2*x_plot5, k:, linewidth1.5, label真实关系: y2x5) plt.xlabel(X) plt.ylabel(y) plt.title(经典OLS vs. 最小截平方和法(LTS)在含异常值数据上的表现) plt.legend() plt.grid(True, alpha0.3) plt.show() # 打印结果对比 print( 模型系数对比 ) print(f真实模型: 截距 5.00, 斜率 2.00) print(f经典OLS估计: 截距 {beta_ols[0]:.2f}, 斜率 {beta_ols[1]:.2f}) print(fLTS估计 (h{h}): 截距 {beta_lts[0]:.2f}, 斜率 {beta_lts[1]:.2f}) print(f\nLTS识别出 {np.sum(inlier_mask)} 个内点{n - np.sum(inlier_mask)} 个异常值。) print(f真实异常值中被LTS正确识别的比例: {np.mean(np.isin(outlier_indices, np.where(~inlier_mask)[0]))*100:.1f}%)运行这段代码你会直观地看到红色的经典OLS回归线如何被四个异常值“拉偏”而蓝色的LTS回归线则几乎完全不受影响紧密贴合了真实的绿色内点群和黑色真实关系线。这就是稳健回归的力量。4. 实战应用场景与选型指南LTS不是万能的它在特定场景下光芒四射在其他场景下可能不如其他方法。理解其适用边界是成为高手的必经之路。4.1 LTS的典型应用场景金融数据清洗与建模金融时间序列中常包含因市场剧烈波动、数据录入错误产生的异常值。在构建风险模型如VaR或量化因子模型前使用LTS进行初步拟合可以识别并处理这些异常值防止它们扭曲风险参数估计。工业质量控制与传感器数据分析生产线上传感器数据可能因设备间歇性故障、电磁干扰产生跳变。用LTS拟合传感器读数与工艺参数的关系可以更稳健地监控过程是否处于统计受控状态避免误报警。地理空间数据与遥感分析在地理加权回归中局部区域的异常观测如因云层覆盖导致的错误遥感反射率会影响局部参数估计。采用LTS框架的稳健地理加权回归能提升反演精度。生物医学与化学计量学在光谱分析如近红外光谱预测成分含量或剂量反应分析中个别样本的制备误差或测量失误会产生离群点。LTS能帮助建立更可靠的校准模型。任何探索性数据分析的初始步骤在对一个新数据集建立复杂模型如神经网络、梯度提升树之前先用简单的LTS线性模型跑一遍。它不仅能提供一个对异常值不敏感的基准其输出的异常值标识本身就是极佳的数据质量诊断报告。4.2 与其他稳健回归方法的对比选型LTS是稳健回归家族的重要成员但非唯一选择。下表对比了几种常见方法方法核心思想优点缺点适用场景最小截平方和法寻找一个子集使其残差平方和最小。崩溃点高可达50%概念直观对y方向异常值极佳。计算量较大需随机采样迭代对高杠杆点X异常的稳健性不如某些方法。数据中存在响应变量(y)异常且异常值比例可能较高时。M-估计用增长慢于平方函数的ρ函数代替平方损失。计算相对高效有现成迭代重加权最小二乘算法。崩溃点较低依赖初始值通常30%对高杠杆点敏感。异常值比例不高且希望计算效率高时。常作为其他方法的初始值。S-估计最小化残差尺度的M-估计该尺度本身具有高崩溃点。同时估计回归系数和残差尺度具有高崩溃点。计算复杂效率有时低于LTS。需要高崩溃点且对效率要求不是极端苛刻时。MM-估计先用高崩溃点方法如LTS获得初始估计和残差尺度再用高效的M-估计进行“抛光”。兼具高崩溃点和高统计效率是当前推荐的综合选择。实现稍复杂需要两个阶段。通用推荐。在需要同时保证稳健性和估计精度时首选。随机森林回归基于决策树集成通过平均多棵树预测来降低方差。对非线性关系建模能力强对异常值有一定包容性无需假设数据分布。模型可解释性差计算和存储成本高是“黑箱”模型。关系复杂、非线性且预测精度优先于模型解释时。选型建议新手入门或快速诊断从LTS开始。它的结果内点/异常点划分非常直观能让你迅速了解数据质量。生产环境或严肃分析优先考虑MM-估计。它在R的robustbase包和Python的statsmodels中都有成熟实现平衡了稳健性和效率。已知异常值主要在y方向LTS表现优异。担心高杠杆点X异常考虑广义M-估计或S-估计它们对设计空间中的异常更稳健。追求极致预测精度且关系复杂可以尝试随机森林等集成方法但务必结合交叉验证并意识到其解释成本。5. 高级话题与性能优化当数据量变大或维度变高时基础的随机采样LTS算法会面临挑战。以下是几个进阶方向。5.1 处理大规模数据算法加速技巧当样本量n很大时计算所有点的残差并排序可能成为瓶颈。可以采用以下策略分块计算将数据随机分成多个块在每个块上独立运行LTS或计算残差然后合并结果。这类似于“分而治之”的思想。改进的随机采样策略不是完全随机采样p1个点而是先通过一些快速异常检测方法如基于马氏距离剔除明显异常的点在“干净”候选集中进行采样提高初始子集的质量减少无效迭代。利用稀疏性如果数据或设计矩阵是稀疏的可以使用专门的稀疏矩阵运算库来加速矩阵乘法和求解。近似算法采用Fast-LTS算法它使用一系列巧妙的初始子集选择方法如六边形等来代替完全随机采样能大幅减少达到相同精度所需的迭代次数。5.2 从线性到非线性LTS思想的扩展经典的LTS针对线性模型。但其“寻找一个干净子集”的核心思想可以推广广义线性模型对于逻辑回归、泊松回归等可以将残差平方和替换为相应的偏差Deviance寻求最小化部分观测的偏差之和。非线性回归对于形如y f(X, β) ε的非线性模型算法框架类似。关键在于对于每个候选子集需要使用非线性优化方法如Levenberg-Marquardt来拟合参数β计算成本会显著增加。分位数回归LTS关注的是中心趋势条件均值的稳健估计。而分位数回归如中位数回归本身对异常值就稳健。可以将LTS思想与分位数回归结合进一步提升其在尾部的稳健性。5.3 统计推断如何为LTS结果计算置信区间经典OLS的统计推断假设检验、置信区间建立在正态误差和同方差等假设上。LTS抛弃了这些假设其抽样分布更加复杂。常用的推断方法有自助法这是最实用、最通用的方法。从原始数据中有放回地重复抽样对每个自助样本计算LTS估计然后用这些自助估计的分布来近似原估计的抽样分布从而计算标准误和置信区间。基于权重的近似推断LTS的最终解可以看作是一个加权最小二乘解内点权重为1异常点权重为0。基于这个加权方案可以推导近似的协方差矩阵但这种方法通常比较粗糙。基于稳健尺度的推断先通过LTS获得残差然后用高崩溃点方法如Qn或MAD估计残差尺度再利用这个尺度进行类似t检验的推断。这种方法相对简单但前提是残差分布大致对称。实操心得在实际项目中如果需要进行严格的统计推断我强烈推荐使用自助法。虽然计算量大但它对模型假设要求最低结果也最可靠。对于大多数探索性分析或预测任务直接报告LTS的点估计和识别出的异常值列表往往已经足够支持决策。6. 常见陷阱、问题排查与实战心得即使理解了原理和算法在实际编码和应用中依然会遇到各种坑。下面是我总结的一些典型问题及解决方案。6.1 算法不收敛或结果不稳定现象多次运行LTS得到的系数差异很大。可能原因与解决随机采样次数n_subsets太少这是最常见的原因。增加n_subsets比如从1000增加到5000或10000给算法更多探索机会。数据中存在大量异常值超过了算法的崩溃点检查你的h值设置。如果真实异常值比例超过(n-h)/nLTS可能失效。尝试减小h值例如从0.75n降到0.6n但要注意这会损失效率。更好的做法是结合业务知识预先审查并处理最明显的异常值。初始子集(p1)点共线或近似共线在多元回归中随机抽到的p1个点可能几乎落在同一个超平面上导致初始OLS拟合失败或不稳定。可以在代码中增加判断如果初始子阵的条件数过大则跳过此次采样。C-Step陷入循环极少数情况下C-Step可能在两个相近的解之间震荡。在代码中设置最大迭代次数max_csteps如10或20和收敛容差如系数变化小于1e-8可以避免无限循环。6.2 误伤与漏报异常值识别不准现象LTS将一些看起来正常的点判为异常或者漏掉了一些明显的异常点。可能原因与解决参数h设置不当h是控制敏感度的阀门。h设得太大模型过于“宽容”会漏报异常值h设得太小模型过于“苛刻”会误伤正常点。没有银弹需要通过可视化如残差图、杠杆值-残差图结合业务理解来调整。可以尝试运行多个h值观察异常点列表的变化选择一个使结果在业务上最合理的h。存在高杠杆点LTS对y方向的异常稳健但对X空间的异常点高杠杆点识别能力有限。一个在X空间远离主体但y值恰好落在回归线上的点可能不会被LTS判为异常但它会极大地影响回归线的斜率。解决方法是结合诊断图。在拟合LTS后计算每个点的杠杆值hat value和LTS残差绘制残差-杠杆图。落在图右上或右下方区域的点需要高度警惕。数据存在集群或分组结构如果数据本身来自多个不同的群体混合分布LTS可能会把其中一个群体全部判为异常。这时线性模型可能已不适用需要考虑混合回归模型或聚类分析。6.3 计算效率低下现象数据量稍大如n10000时程序运行非常慢。优化策略向量化操作确保核心计算如矩阵乘法X.dot(beta)、残差计算使用NumPy的向量化操作避免Python层级的循环。使用更快的排序算法np.partition比np.sort更快因为我们只需要最小的h个值而不是全排序。我们的示例代码已经使用了np.argpartition。降维如果自变量很多可以考虑先使用主成分分析进行降维在低维空间进行稳健拟合但这会损失可解释性。调用优化库对于生产环境直接使用高度优化的库如Pythonstatsmodels中的RLM提供了M、S、MM估计或R语言中的robustbase和MASS包。它们底层由C/Fortran实现效率远高于自编的Python循环。6.4 与现有工作流的整合LTS不应是一个孤立的步骤。一个完整的稳健建模流程应该是数据可视化绘制散点图、箱线图对数据分布有一个直观认识。运行经典OLS作为一个基准并计算其残差、杠杆值等诊断统计量。运行LTS或MM估计获得稳健的系数估计和异常值候选列表。对比分析比较OLS和LTS的系数。如果差异巨大说明数据受异常值影响严重LTS结果更可信。诊断调查仔细审查被LTS标记为异常的点。结合业务逻辑判断它们是数据错误需修正或删除、特殊事件需单独建模还是模型缺陷如缺失重要变量、非线性。决策与报告根据诊断结果决定是清洗数据后使用OLS还是直接报告和使用LTS模型。在报告中必须明确说明使用了稳健方法以及处理异常值的策略。最后记住一点任何模型都是对现实的简化。LTS帮助我们抵御异常值的干扰但它不能替代对业务逻辑的深入理解。最强大的模型永远是“统计方法”与“领域知识”的结合。当你看到一个异常值时第一反应不应该是简单地删除它而是问一句“这个点为什么会在这里” 答案或许会引领你发现一个全新的业务洞察或数据质量问题。