M估计双平方法（M-Estimation Bisquare Method）是一种鲁棒回归技术，旨在通过动态调整数据点的权重来减少异常值对模型参数估计的影响。其核心思想是赋予正常数据点更高的权重，而对远离模型的异常点赋予低或零权重。以下从核心原理、数学公式、步骤、优缺点及应用场景等方面详细解释：

一、核心原理

1. M估计（M-Estimator）框架

• 目标：替代传统最小二乘法（OLS）的平方损失函数，使用更鲁棒的损失函数，降低异常值的敏感性。

• 损失函数选择：双平方（Bisquare/Tukey）函数，其特点是当残差超过阈值时，损失不再增长，权重降为零。

1. 双平方权重函数

• 通过非线性权重分配，使大残差对应的权重趋近于零，从而在迭代过程中逐步剔除异常值的影响。

二、数学公式

1. 双平方损失函数（Bisquare Loss Function）  
   

• ( $ r = \frac{y_i - \theta^T x_i}{\sigma} $ )：标准化残差（残差除以尺度估计，如中位数绝对偏差MAD）。

• ( c )：调节常数（通常取4.685），控制对异常值的敏感度。

1. 权重函数（Weight Function）  
   由损失函数导出的权重函数为：
   $$
   w(r) = \begin{cases}
   \left(1 - \left(\frac{r}{c}\right)^2\right)^2 & \text{if } |r| \leq c, \
   0 & \text{otherwise}.
   \end{cases}
   $$

• 特点：当残差绝对值超过阈值 ( c )，权重直接置零，完全排除异常值。

三、算法步骤

1. 初始化

• 使用普通最小二乘法（OLS）或中位数回归得到初始参数估计 ( $\theta^{(0)} $)。

• 计算残差的尺度估计 ($ \sigma $)，通常采用中位数绝对偏差（MAD）：
  $$
  \sigma = 1.4826 \cdot \text{median}(|r_i - \text{median}(r_i)|).
  $$

1. 迭代重加权最小二乘法（IRLS）  
   重复以下步骤直至收敛：

• 步骤1：计算当前残差 ($ r_i^{(k)} = y_i - \theta^{(k)T} x_i $)，并标准化 ($ r_i^{(k)} = r_i^{(k)} / \sigma$ ).

• 步骤2：根据权重函数计算每个数据点的权重 ($ w_i^{(k)} = w(r_i^{(k)})$ ).

• 步骤3：通过加权最小二乘更新参数：
  $$
  \theta^{(k+1)} = \arg\min_{\theta} \sum_{i=1}^n w_i^{(k)} (y_i - \theta^T x_i)^2.
  $$

• 步骤4：检查参数变化量是否小于阈值（如 ($ |\theta^{(k+1)} - \theta^{(k)}| < 10^{-6}$ )），若满足则停止迭代。

四、关键参数与调节

1. 调节常数 ( c )

• 作用：控制对异常值的剔除程度。较小的 ( c ) 会剔除更多数据点（更鲁棒但可能损失信息），较大的 ( c ) 则保留更多点（接近OLS）。

• 经验值：通常取 ( c = 4.685 )，在正态分布假设下可保持95%的统计效率。

1. 尺度估计 ( $\sigma$ )

• 必须使用鲁棒的尺度估计（如MAD），避免异常值对方差估计的污染。

五、优缺点

优点

1. 强鲁棒性：对极端异常值完全免疫（权重为零）。

2. 数据驱动：无需预先假设异常值比例，自动调整权重。

3. 灵活性：适用于多种噪声分布（如厚尾分布、污染数据）。

缺点

1. 计算复杂度：需迭代求解，计算量高于OLS。

2. 参数敏感：结果依赖 ( c ) 的选择，需通过交叉验证或经验调整。

3. 非凸优化：双平方损失函数非凸，可能陷入局部最优（需良好初始值）。

六、与其他鲁棒方法的对比

七、应用场景

1. 金融数据分析

• 处理股票价格中的极端波动或黑天鹅事件。

1. 传感器信号处理

• 消除传感器噪声中的偶发异常测量值。

1. 计算机视觉

• 鲁棒拟合图像中的几何模型（如去除误匹配点）。

1. 生物医学统计

• 分析存在离群值的实验数据（如基因表达异常样本）。

八、代码示例（Python）

使用 statsmodels 库实现双平方M估计：

import statsmodels.api as sm
import numpy as np

# 生成含异常值的数据
np.random.seed(0)
X = np.random.randn(100, 2)
y = 1.5 * X[:, 0] + 0.8 * X[:, 1] + np.random.randn(100)
y[::10] += 20  # 添加异常值（每10个点加一个极端值）

# 双平方M估计（TukeyBiweight）
model = sm.RLM(y, sm.add_constant(X), M=sm.robust.norms.TukeyBiweight())
results = model.fit()

# 输出结果
print("估计参数（含截距项）：", results.params)
print("收敛信息：", results.converged)

九、总结

M估计双平方法通过动态权重调整，在保留正常数据信息的同时彻底排除极端异常值，是处理污染数据的有效工具。其核心优势在于平衡了鲁棒性和计算效率，适用于实际场景中常见的非高斯噪声和离群点问题。然而，使用时需注意参数调节和初始值选择，以确保收敛到全局最优解。