最小截平方和法(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采用了完全不同的策略。它的核心思想是：

1. 从全部n个数据点中，选出所有可能的、包含h个点的子集（h通常大于n/2，例如取n的75%）。
2. 对每一个这样的子集，用普通最小二乘法拟合一个模型，并计算这个模型在整个数据集上的残差平方和，但只取其中最小的h个残差平方进行求和。
3. 最终，选择那个使得“最小h个残差平方和”最小的子集所对应的模型参数，作为LTS的估计结果。

简单来说，LTS不是在找一条让“所有点”的误差平方和最小的线，而是在找一条让“大多数好点”的误差平方和最小的线。它主动忽略掉那些误差最大的“坏点”（异常值）。参数h控制了模型的“稳健性”和“效率”之间的权衡：h越大（越接近n），模型效率越高（估计更精确），但稳健性越差；h越小，模型越稳健，但效率会损失。通常，h的取值在floor((n+p+1)/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算法实现包含以下步骤：

1. 初始化与参数设置：

   • 输入：数据矩阵X（包含常数项，如果有的话），响应向量y，以及参数h（子集大小）。
   • 设置最大迭代次数max_iter（如500次）和随机子集采样数量n_subsets（如1000个）。这些参数是为了在计算时间和求解精度间取得平衡。

2. 随机子集采样与初步拟合：

   • 循环进行n_subsets次： a. 从全部n个观测中，完全随机地抽取一个大小为p+1的子集（p是自变量个数，p+1是能唯一确定一个线性模型的最小样本数）。 b. 用这个小子集的数据，通过普通最小二乘法计算出一组初始回归系数β_init。 c. 用这组β_init去预测所有n个点的值，计算所有n个残差r_i = y_i - X_iβ_init。 d. 对这n个残差取绝对值并排序，找出残差绝对值最小的前h个观测点，构成一个新的“候选干净子集”。

3. 基于候选子集的精确拟合与迭代改进：

   • 对上一步得到的“候选干净子集”（大小为h），使用普通最小二乘法进行拟合，得到一组新的系数β_candidate。
   • 计算使用β_candidate时，所有n个点的残差，并同样取绝对值最小的前h个残差平方和，记为SSR_h。
   • 记录下到目前为止，得到的最小的SSR_h及其对应的系数β_best。

4. C-Step迭代收敛：

   • 对于表现较好的候选解（例如SSR_h较小的前10%），可以进行称为“C-Step”的迭代改进： a. 给定当前系数β，计算全样本残差，选出残差绝对值最小的h个点构成新子集。 b. 用新子集重新拟合OLS，得到新的β。 c. 重复a和b，直到子集成员不再变化，或SSR_h不再减小。
   • 这个过程能快速将随机得到的初始解“抛光”到局部最优。

5. 输出最终结果：

   • 在所有随机采样和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 < 1000），1000-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_subsets=1000, max_csteps=10, random_state=42): """ 最小截平方和法（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(f"h 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): # 随机选择 p+1 个点作为初始子集 (确保非奇异) random_indices = np.random.choice(n_samples, size=n_features, replace=False) X_sub = X[random_indices] y_sub = y[random_indices] # 2. 初始OLS拟合 (使用伪逆避免奇异矩阵问题) try: beta_init = np.linalg.lstsq(X_sub, y_sub, rcond=None)[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, rcond=None)[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, rcond=None)[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_subsets=500) y_pred_lts = X_design.dot(beta_lts) # --- 可视化对比 --- plt.figure(figsize=(12, 6)) plt.scatter(X, y, alpha=0.6, label='数据点 (含异常值)', c='gray') plt.scatter(X[inlier_mask], y[inlier_mask], alpha=0.8, label='LTS识别的内点', c='green') plt.scatter(X[outlier_indices], y[outlier_indices], marker='x', s=100, label='真实异常值位置', c='red', linewidths=2) 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--', linewidth=2, label=f'经典OLS: y={beta_ols[1]:.2f}x+{beta_ols[0]:.2f}') plt.plot(x_plot, X_plot_design.dot(beta_lts), 'b-', linewidth=2, label=f'LTS (h={h}): y={beta_lts[1]:.2f}x+{beta_lts[0]:.2f}') plt.plot(x_plot, 2*x_plot+5, 'k:', linewidth=1.5, label='真实关系: y=2x+5') plt.xlabel('X') plt.ylabel('y') plt.title('经典OLS vs. 最小截平方和法(LTS)在含异常值数据上的表现') plt.legend() plt.grid(True, alpha=0.3) plt.show() # 打印结果对比 print("=== 模型系数对比 ===") print(f"真实模型: 截距 = 5.00, 斜率 = 2.00") print(f"经典OLS估计: 截距 = {beta_ols[0]:.2f}, 斜率 = {beta_ols[1]:.2f}") print(f"LTS估计 (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的典型应用场景

1. 金融数据清洗与建模：金融时间序列中常包含因市场剧烈波动、数据录入错误产生的异常值。在构建风险模型（如VaR）或量化因子模型前，使用LTS进行初步拟合，可以识别并处理这些异常值，防止它们扭曲风险参数估计。
2. 工业质量控制与传感器数据分析：生产线上传感器数据可能因设备间歇性故障、电磁干扰产生跳变。用LTS拟合传感器读数与工艺参数的关系，可以更稳健地监控过程是否处于统计受控状态，避免误报警。
3. 地理空间数据与遥感分析：在地理加权回归中，局部区域的异常观测（如因云层覆盖导致的错误遥感反射率）会影响局部参数估计。采用LTS框架的稳健地理加权回归能提升反演精度。
4. 生物医学与化学计量学：在光谱分析（如近红外光谱预测成分含量）或剂量反应分析中，个别样本的制备误差或测量失误会产生离群点。LTS能帮助建立更可靠的校准模型。
5. 任何探索性数据分析的初始步骤：在对一个新数据集建立复杂模型（如神经网络、梯度提升树）之前，先用简单的LTS线性模型跑一遍。它不仅能提供一个对异常值不敏感的基准，其输出的异常值标识本身就是极佳的数据质量诊断报告。

4.2 与其他稳健回归方法的对比选型

LTS是稳健回归家族的重要成员，但非唯一选择。下表对比了几种常见方法：

选型建议：

• 新手入门或快速诊断：从LTS开始。它的结果（内点/异常点划分）非常直观，能让你迅速了解数据质量。
• 生产环境或严肃分析：优先考虑MM-估计。它在R的robustbase包和Python的statsmodels中都有成熟实现，平衡了稳健性和效率。
• 已知异常值主要在y方向：LTS表现优异。
• 担心高杠杆点（X异常）：考虑广义M-估计或S-估计，它们对设计空间中的异常更稳健。
• 追求极致预测精度且关系复杂：可以尝试随机森林等集成方法，但务必结合交叉验证，并意识到其解释成本。

5. 高级话题与性能优化

当数据量变大或维度变高时，基础的随机采样LTS算法会面临挑战。以下是几个进阶方向。

5.1 处理大规模数据：算法加速技巧

当样本量n很大时，计算所有点的残差并排序可能成为瓶颈。可以采用以下策略：

• 分块计算：将数据随机分成多个块，在每个块上独立运行LTS（或计算残差），然后合并结果。这类似于“分而治之”的思想。
• 改进的随机采样策略：不是完全随机采样p+1个点，而是先通过一些快速异常检测方法（如基于马氏距离）剔除明显异常的点，在“干净”候选集中进行采样，提高初始子集的质量，减少无效迭代。
• 利用稀疏性：如果数据或设计矩阵是稀疏的，可以使用专门的稀疏矩阵运算库来加速矩阵乘法和求解。
• 近似算法：采用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，得到的系数差异很大。
• 可能原因与解决：
  1. 随机采样次数n_subsets太少：这是最常见的原因。增加n_subsets（比如从1000增加到5000或10000），给算法更多探索机会。
  2. 数据中存在大量异常值，超过了算法的崩溃点：检查你的h值设置。如果真实异常值比例超过(n-h)/n，LTS可能失效。尝试减小h值（例如从0.75n降到0.6n），但要注意这会损失效率。更好的做法是结合业务知识，预先审查并处理最明显的异常值。
  3. 初始子集(p+1)点共线或近似共线：在多元回归中，随机抽到的p+1个点可能几乎落在同一个超平面上，导致初始OLS拟合失败或不稳定。可以在代码中增加判断，如果初始子阵的条件数过大，则跳过此次采样。
  4. C-Step陷入循环：极少数情况下，C-Step可能在两个相近的解之间震荡。在代码中设置最大迭代次数max_csteps（如10或20）和收敛容差（如系数变化小于1e-8）可以避免无限循环。

6.2 误伤与漏报：异常值识别不准

• 现象：LTS将一些看起来正常的点判为异常，或者漏掉了一些明显的异常点。
• 可能原因与解决：
  1. 参数h设置不当：h是控制敏感度的阀门。h设得太大，模型过于“宽容”，会漏报异常值；h设得太小，模型过于“苛刻”，会误伤正常点。没有银弹，需要通过可视化（如残差图、杠杆值-残差图）结合业务理解来调整。可以尝试运行多个h值，观察异常点列表的变化，选择一个使结果在业务上最合理的h。
  2. 存在高杠杆点：LTS对y方向的异常稳健，但对X空间的异常点（高杠杆点）识别能力有限。一个在X空间远离主体，但y值恰好落在回归线上的点，可能不会被LTS判为异常，但它会极大地影响回归线的斜率。解决方法是结合诊断图。在拟合LTS后，计算每个点的杠杆值（hat value）和LTS残差，绘制残差-杠杆图。落在图右上或右下方区域的点，需要高度警惕。
  3. 数据存在集群或分组结构：如果数据本身来自多个不同的群体（混合分布），LTS可能会把其中一个群体全部判为异常。这时，线性模型可能已不适用，需要考虑混合回归模型或聚类分析。

6.3 计算效率低下

• 现象：数据量稍大（如n>10000）时，程序运行非常慢。
• 优化策略：
  1. 向量化操作：确保核心计算（如矩阵乘法X.dot(beta)、残差计算）使用NumPy的向量化操作，避免Python层级的循环。
  2. 使用更快的排序算法：np.partition比np.sort更快，因为我们只需要最小的h个值，而不是全排序。我们的示例代码已经使用了np.argpartition。
  3. 降维：如果自变量很多，可以考虑先使用主成分分析进行降维，在低维空间进行稳健拟合，但这会损失可解释性。
  4. 调用优化库：对于生产环境，直接使用高度优化的库，如Pythonstatsmodels中的RLM（提供了M、S、MM估计），或R语言中的robustbase和MASS包。它们底层由C/Fortran实现，效率远高于自编的Python循环。

6.4 与现有工作流的整合

LTS不应是一个孤立的步骤。一个完整的稳健建模流程应该是：

1. 数据可视化：绘制散点图、箱线图，对数据分布有一个直观认识。
2. 运行经典OLS：作为一个基准，并计算其残差、杠杆值等诊断统计量。
3. 运行LTS（或MM估计）：获得稳健的系数估计和异常值候选列表。
4. 对比分析：比较OLS和LTS的系数。如果差异巨大，说明数据受异常值影响严重，LTS结果更可信。
5. 诊断调查：仔细审查被LTS标记为异常的点。结合业务逻辑，判断它们是数据错误（需修正或删除）、特殊事件（需单独建模）还是模型缺陷（如缺失重要变量、非线性）。
6. 决策与报告：根据诊断结果，决定是清洗数据后使用OLS，还是直接报告和使用LTS模型。在报告中，必须明确说明使用了稳健方法以及处理异常值的策略。

最后，记住一点：任何模型都是对现实的简化。LTS帮助我们抵御异常值的干扰，但它不能替代对业务逻辑的深入理解。最强大的模型，永远是“统计方法”与“领域知识”的结合。当你看到一个异常值时，第一反应不应该是简单地删除它，而是问一句：“这个点为什么会在这里？” 答案或许会引领你发现一个全新的业务洞察或数据质量问题。