稳健M估计

稳健M估计（robust M-estimation）通过选择对大残差增长较缓的损失或有界得分，降低少量极端观测对拟合的作用。M估计本身是一个更大的类别，普通最小二乘也在其中；稳健性取决于损失、尺度、观测设计以及所估计的总体目标。

以教学读数 $(0, 0, 1, 1, 10)$ 为例，等权平均为 $2.4$ ，中位数为一。若希望小误差仍按平方损失处理，而读数十的拉动不再随其大小无界增加，可以对残差得分作限幅。

固定尺度的位置方程

设尺度 $s > 0$ 事先给定，对位置 $a$ 最小化 $\sum_{i} ρ ((x_{i} - a) / s) .$ 若 $ρ$ 可微，记 $ψ = ρ^{'}$ ，一阶条件为 $\sum_{i} ψ ((x_{i} - a) / s) = 0$ 。胡贝尔得分取 $ψ_{c} (u) = {\begin{matrix} - c, & u < - c, \\ u, & | u | \leq c, \\ c, & u > c, \end{matrix} c > 0.$ 它保持单调，故对应损失凸；得分和随 $a$ 不增，零点就是全局最小点。一般数据可能产生一段零点，不能不检查便宣称严格唯一。位置M估计、得分与回归杠杆的联系见Ruckstuhl讲义。^[1]

阈值一点五的胡贝尔得分在负一点五至正一点五之间是斜率为一的直线，外侧保持正负一点五 — 横轴是标准化残差，纵轴是单条观测的得分。中间段保留小残差的大小差异，两端把得分限制在固定幅值；极端读数仍以非零的限幅贡献参与估计。

主例取 $s = 1, c = 3 / 2$ 。若 $0 \leq a \leq 1$ ，前四个残差均在未截断区，最后一个贡献固定为 $3 / 2$ ，方程成为 $- 2 a + 2 (1 - a) + 3 / 2 = 0.$ 解得 $\hat{a} = 7 / 8 = 0.875$ ，确实在所假设区间内。附近得分和的斜率为负四，故该零点唯一；凸性保证全局最优。十若改成一百，尾部贡献仍为 $3 / 2$ ，结果保持 $7 / 8$ 。

五个读数的胡贝尔得分和在0.875处过零，极端点十改为一百后在这一区域不改变方程 — 横轴为候选位置，纵轴为限幅残差和。零点附近四个普通观测决定斜率，极端观测只贡献固定上限。

估计的总体位置与影响函数

对应总体目标应明确定义为 $𝔼 ψ_{c} ((X - a) / s) = 0$ 的适当根。若误差分布关于 $a_{0}$ 对称、 $ψ_{c}$ 为奇函数，则在 $a_{0}$ 处期望为零；若该总体方程还有唯一性，才可称目标与对称中心一致。对偏斜分布，这个目标通常不同于总体均值。

假设分布在两截断点无原子、期望可微，且 $A = 𝔼 {ψ_{c}}^{'} ((X - a_{0}) / s) > 0$ 。沿定点污染路径，令根为 $a_{ε}$ ，对方程求导得 $- \frac{A}{s} {\frac{d a_{ε}}{d ε} |}_{0} + ψ_{c} ((z - a_{0}) / s) = 0.$ 因此影响函数为 $s ψ_{c} ((z - a_{0}) / s) / A$ ，绝对值至多 $s c / A$ 。截断点处的导数、分母是否为零以及尺度是否固定，都是这项证明的条件。尺度若也从数据估计，需要把尺度方程一起线性化。

单调求根的可复算程序

以下代码返回固定尺度位置最优解的一项数值近似，原始样本不被修改。最小读数处得分和非负，最大读数处非正；中点若得分为正，根应在右半段，否则在左半段。每步保留含根区间，精确实数运算下区间宽度减半。

Python 3

import math

def huber_location(values, cutoff=1.5, scale=1.0, steps=80):
    x = tuple(float(v) for v in values)
    cutoff, scale = float(cutoff), float(scale)
    if not x or any(not math.isfinite(v) for v in x):
        raise ValueError("nonempty finite sample required")
    if not (math.isfinite(cutoff) and cutoff > 0
            and math.isfinite(scale) and scale > 0):
        raise ValueError("positive finite cutoff and scale required")
    if type(steps) is not int or steps < 1:
        raise ValueError("positive integer steps required")
    from fractions import Fraction
    exact_x = tuple(Fraction(v) for v in x)
    exact_cutoff, exact_scale = Fraction(cutoff), Fraction(scale)
    lo, hi = min(x), max(x)
    for _ in range(steps):
        mid = max(lo, min(hi, lo / 2 + hi / 2))
        exact_mid = Fraction(mid)
        score = sum((max(-exact_cutoff, min(exact_cutoff,
                     (v-exact_mid)/exact_scale)) for v in exact_x), Fraction(0))
        if score == 0:
            return mid
        if score > 0:
            lo = mid
        else:
            hi = mid
    return max(lo, min(hi, lo / 2 + hi / 2))

print(huber_location([0, 0, 1, 1, 10]))  # 0.875
print(huber_location([0, 0, 1, 1, 100])) # 0.875

执行 $K$ 轮时，需要 $O (n K)$ 次有理数算术，样本副本占 $O (n)$ 个数的存储；整数位数对时间和空间的影响还须另计。有限浮点输入本身都是精确的二进制有理数，Fraction 按这些已转换的值计算得分符号，避免极大差值溢出、极小尺度比下溢而误判为零。二分点和返回值仍用浮点数，最终可表示精度及迭代轮数限制返回位置的准确度。这一教学实现优先保证符号判断，速度不等同于专用稳健统计库。若零点构成区间，程序可以逼近其中一项，而不是声称目标严格凸。

变式。 将十改成一，数据变为 $(0, 0, 1, 1, 1)$ ，均值为 $0.6$ ，全部残差均小于 $1.5$ ，所以M估计也为 $0.6$ 。限幅只在残差达到阈值后起作用，不是对所有样本都强行取中位数。

在回归中，方程常含 $\sum_{i} x_{i} ψ (r_{i} / s)$ 。即使 $ψ$ 有界，任意大的解释变量 $x_{i}$ 仍可能制造高杠杆作用。位置模型的有界影响结论不能原样搬成“对任意回归异常点都稳健”。

参考文献

↑ Andreas Ruckstuhl，Robust Fitting: Introduction to Robust Estimation Techniques，第15—18张位置M估计；后续回归M估计中得分含设计变量的讨论。

[1] Andreas Ruckstuhl，Robust Fitting: Introduction to Robust Estimation Techniques，第15—18张位置M估计；后续回归M估计中得分含设计变量的讨论。

[1]