M-估计量(M-Estimator)
在经典参数估计中,我们通常先指定一个概率模型,再通过极大似然估计(Maximum Likelihood Estimation,MLE)得到参数估计量。这一方法在模型设定正确时往往具有很高的统计效率,但当真实数据存在厚尾、测量错误或少量异常观测时,估计结果可能对模型偏离非常敏感。稳健统计(Robust Statistics)研究的核心问题,正是如何在保留理想模型下统计效率的同时,控制模型偏离对估计结果的影响。
M-估计量(M-Estimator)为这一问题提供了统一的分析框架。它以损失函数构造估计量,再通过得分函数、估计方程和渐近理论分析估计量的行为。Peter J. Huber 在 1964 年发表的 Robust Estimation of a Location Parameter 是这一领域的重要奠基性论文:文章研究了凸损失 M-估计量的渐近性质,并证明了在特定正态污染邻域和极小极大渐近方差准则下,截断线性的得分函数具有最优意义。本文以该论文为理论主线,结合后来形成的影响函数表述与现代数值计算方法,建立一套可用于继续学习稳健回归和广义 M-估计的知识体系。
从样本均值出发:估计为何会受到异常值影响
正态模型下的最小二乘估计
设 $X_1,\ldots,X_n$ 是独立同分布样本,来自正态位置模型 $X_i\sim N(\theta,\sigma^2)$,其中位置参数 $\theta$ 未知,尺度 $\sigma>0$ 已知。其对数似然函数为:
\[\ell(\theta) =-\frac{n}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2}\sum_{i=1}^{n}(X_i-\theta)^2.\]最大化对数似然等价于最小化平方残差和,因此正态位置模型下的极大似然估计量满足:
\[\hat\theta_{\mathrm{MLE}} =\underset{\theta\in\mathbb R}{\arg\min}\, \sum_{i=1}^{n}(X_i-\theta)^2 =\bar X.\]在正态模型正确时,样本均值的方差为 $\sigma^2/n$,并且达到相应的 Cramér–Rao 下界。问题在于,平方损失对残差的惩罚随残差幅度二次增长,导致极端观测具有很强的参数牵引力。假设把第 $j$ 个观测 $X_j$ 替换为一个极端值 $M$,那么新的样本均值变为:
\[\bar X_{\mathrm{new}} =\bar X+\frac{M-X_j}{n}.\]当 $M\to+\infty$ 时,均值的变化也趋于无穷。这并不是优化算法不稳定,而是样本均值的统计定义决定了它会对单个极端值保持无界的响应。即使其他样本完全不变,一个异常观测也能把估计结果推得任意远。
中位数的稳健性与效率代价
一个自然替代方案是最小化绝对偏差损失。样本中位数属于下列最优化问题的解集合:
\[\hat\theta_{\mathrm{med}} \in\underset{\theta\in\mathbb R}{\arg\min}\, \sum_{i=1}^{n}\lvert X_i-\theta\rvert.\]绝对值损失只随残差线性增长,因此一个观测离中心越来越远,不会像平方损失那样产生持续增大的边际惩罚。在正态分布下,样本中位数的渐近方差约为 $\pi\sigma^2/(2n)$,而样本均值的方差为 $\sigma^2/n$;因此,以样本均值为基准,中位数的渐近相对效率为 $2/\pi\approx0.637$。中位数在极端污染下具有更强的抗干扰能力,却在纯正态模型下损失了一部分估计精度。
由此产生一个关键的统计设计问题:能否使小残差按照平方损失处理,以保留中心观测中的信息,同时使大残差只受到线性惩罚,从而降低极端观测的影响?Huber M-估计正是这一思路的经典实现。
M-估计量的定义:从损失函数到估计方程
位置参数的 M-估计
字母 M 表示 Maximum Likelihood Type,强调这种估计方式在结构上类似极大似然估计,但损失函数不必来自某个完全指定的概率密度。给定损失函数 $\rho:\mathbb R\to\mathbb R$,一维位置参数的 M-估计量定义为:
\[\boxed{ \hat\theta_n \in\underset{\theta\in\mathbb R}{\arg\min}\, \sum_{i=1}^{n}\rho(X_i-\theta) }\]实际分析中常对残差进行尺度标准化,即使用 $u_i=(X_i-\theta)/s$,其中 $s>0$ 是已知尺度或预先给定的尺度估计值。此时的目标函数为:
\[Q_n(\theta) =\sum_{i=1}^{n} \rho\left(\frac{X_i-\theta}{s}\right).\]若 $\rho$ 可微,定义得分函数 $\psi(u)=\rho’(u)$。对 $Q_n$ 关于 $\theta$ 求导,得到:
\[\frac{\partial Q_n(\theta)}{\partial\theta} =-\frac{1}{s}\sum_{i=1}^{n} \psi\left(\frac{X_i-\theta}{s}\right).\]因此,内部极小点必须满足 M-估计方程(Estimating Equation):
\[\boxed{ \sum_{i=1}^{n} \psi\left(\frac{X_i-\hat\theta_n}{s}\right)=0 }\]这一推导说明了 $\rho$ 与 $\psi$ 各自的角色:损失函数 $\rho$ 定义优化目标,得分函数 $\psi$ 决定每个观测如何推动参数估计。在凸且可微的情形,满足一阶条件的解是全局最小化解;若目标函数并非严格凸,解可能不唯一;若使用非凸损失,则估计方程的零点甚至可能只是局部驻点。
均值和中位数都属于 M-估计框架
对于平方损失 $\rho(u)=u^2/2$,有 $\psi(u)=u$,估计方程化为 $\sum_i(X_i-\theta)=0$,其解就是样本均值。对于绝对值损失 $\rho(u)=\lvert u\rvert$,除零点外的导数为 $\psi(u)=\operatorname{sign}(u)$,利用次梯度可以得到样本中位数的最优化条件。前者使用残差幅度,后者主要使用残差的符号,因此两者对极端观测的响应方式不同。
必须强调,M-估计量并不自动意味着稳健估计量。样本均值也是 M-估计量,却对极端观测具有无界敏感性。真正决定稳健行为的是得分函数在尾部如何增长,以及总体估计方程在目标参数附近是否可识别。
Huber 损失:在中心区域保留效率,在尾部限制影响
分段损失函数与截断得分
Huber 损失使用阈值 $c>0$ 区分中心残差与尾部残差,其标准形式为:
\[\boxed{ \rho_c(u)= \begin{cases} \frac{1}{2}u^2, & \lvert u\rvert\le c,\\ c\lvert u\rvert-\frac{1}{2}c^2, & \lvert u\rvert>c. \end{cases} }\]这里的分段函数只在中心区使用平方损失;当残差超过阈值后,损失函数改为线性增长。常数项 $-c^2/2$ 使两段在连接点连续,并保证一阶导数连续。Huber 损失的得分函数为:
\[\boxed{ \psi_c(u)= \begin{cases} -c, & u<-c,\\ u, & -c\le u\le c,\\ c, & u>c. \end{cases} }\]等价地,可以写作 $\psi_c(u)=\operatorname{clip}(u,-c,c)$。与平方损失的 $\psi(u)=u$ 不同,Huber 得分函数满足 $\lvert\psi_c(u)\rvert\le c$。因此,无论观测有多极端,其在估计方程中产生的得分贡献都被限制在固定范围内。这种机制称为得分截断(Score Clipping),并不是将异常观测从样本中删除。
截断阈值控制效率与稳健性的权衡
当 $c$ 较大时,大多数残差处于二次损失区域,Huber 估计与样本均值接近,理想正态模型下的效率较高,但对异常值的抑制相对较弱;当 $c$ 较小时,较多观测进入线性损失区域,估计行为更接近中位数型估计,异常值的得分贡献更早受到限制。
在极限意义下,随着 $c\to\infty$,对于每个固定的残差 $u$,有 $\psi_c(u)\to u$,所以 Huber 估计趋向最小二乘估计。另一方面,当 $c\to0$ 时,经过对整个估计方程乘以 $1/c$ 的重新缩放,得分函数在非零残差处趋向符号函数,这解释了 Huber 估计与中位数之间的联系。这里的极限解释需要注意样本取值、尺度与解的唯一性等条件,不能将它误认为所有数据下的无条件逐点等价。
渐近理论:M-估计量的方差如何确定
总体估计方程与一致性
设观测变量 $X$ 来自分布 $F$,固定尺度为 $s$。与样本估计方程相对应的总体 M-估计泛函 $T(F)$ 满足:
\[\mathbb E_F\left[ \psi\left(\frac{X-T(F)}{s}\right) \right]=0.\]记总体方程的唯一解为 $\theta_0=T(F)$。如果经验目标函数或经验估计方程能够在适当意义下一致逼近总体对应对象,同时满足可识别性和正则条件,则可以证明 $\hat\theta_n\xrightarrow{P}\theta_0$。对于关于 $\theta_0$ 对称的分布和奇函数 $\psi$,总体得分在 $\theta_0$ 处通常为零,因此 $\theta_0$ 可以解释为分布的对称中心。
需要特别区分总体 M-估计泛函与通常意义下的均值。当分布不对称或受到单侧污染时,总体估计方程的根不一定等于 $\mathbb E_F[X]$,甚至不一定等于污染前模型的位置参数。因此,关于一致性的论断必须明确其收敛目标,不能只说估计量“收敛到真实位置”而不说明位置参数如何定义。
渐近正态性的推导
设 $Z=(X-\theta_0)/s$,定义 $A=\mathbb E_F[\psi’(Z)]$ 和 $B=\mathbb E_F[\psi(Z)^2]$。假设 $\mathbb E_F[\psi(Z)]=0$、$0<A<\infty$、$B<\infty$,并且估计量一致、估计方程允许适当的一阶展开。写出样本估计方程:
\[0=\frac{1}{n}\sum_{i=1}^{n} \psi\left(\frac{X_i-\hat\theta_n}{s}\right).\]在 $\theta_0$ 附近进行一阶展开,并令 $Z_i=(X_i-\theta_0)/s$,可以得到:
\[0= \frac{1}{n}\sum_{i=1}^{n}\psi(Z_i) -\frac{\hat\theta_n-\theta_0}{s} \left(\frac{1}{n}\sum_{i=1}^{n}\psi'(Z_i)\right) +o_p(n^{-1/2}).\]重新整理后,估计误差满足渐近线性近似:
\[\sqrt n(\hat\theta_n-\theta_0) =s\frac{ n^{-1/2}\sum_{i=1}^{n}\psi(Z_i) }{ n^{-1}\sum_{i=1}^{n}\psi'(Z_i) }+o_p(1).\]分子可以应用中心极限定理,分母可以应用大数定律,再结合 Slutsky 定理,得到:
\[\boxed{ \sqrt n(\hat\theta_n-\theta_0) \xrightarrow{d} N\left(0,\,s^2\frac{B}{A^2}\right) }\]因此,M-估计量的大样本方差近似为:
\[\boxed{ \operatorname{Var}(\hat\theta_n) \approx \frac{s^2}{n} \frac{\mathbb E_F[\psi(Z)^2]} {\left(\mathbb E_F[\psi'(Z)]\right)^2} }\]该公式是一维夹心方差(Sandwich Variance)的典型形式。分子衡量单个观测得分的随机波动,分母衡量总体估计方程对参数变化的平均响应强度;如果平均导数过小,即使得分函数本身有界,估计量的方差也可能变大。严格地说,Huber 得分在 $\pm c$ 处不可微,但当分布在这些位置没有点质量时,可以使用几乎处处导数,并在适当正则条件下得到相同的渐近结论。
正态模型下的渐近相对效率
对于 $Z\sim N(0,1)$,若使用平方损失的得分 $\psi(Z)=Z$,则 $A=1$、$B=1$,渐近方差就是 $\sigma^2/n$。对于 Huber 得分函数,令尺度 $s=\sigma$,在标准正态分布下有:
\[A(c)=\mathbb P(\lvert Z\rvert\le c) =2\Phi(c)-1,\]其中 $\Phi$ 表示标准正态分布函数。相应的得分二阶矩为:
\[B(c) =\mathbb E\left[ Z^2\mathbf 1_{\{\lvert Z\rvert\le c\}} \right] +c^2\mathbb P(\lvert Z\rvert>c).\]利用标准正态密度 $\phi$ 的积分公式,可进一步写为:
\[B(c) =2\Phi(c)-1-2c\phi(c) +2c^2[1-\Phi(c)].\]由于样本均值在正态位置模型下的渐近方差系数为一,Huber 估计相对于样本均值的渐近相对效率为:
\[\boxed{ \operatorname{ARE}(c) =\frac{A(c)^2}{B(c)} }\]当 $c\approx1.345$ 时,正态模型下的渐近相对效率约为 $95\%$。这意味着该阈值在理想正态模型下只付出相对较小的效率代价,却使异常观测的得分贡献受到截断。需要强调的是,$1.345$ 是由给定效率目标产生的常用选择,而不是对所有污染机制都最优的固定常数。
影响函数:从局部污染理解稳健性
影响函数的定义
现代稳健统计通过影响函数(Influence Function)衡量估计泛函对无限小比例污染的局部响应。设 $\Delta_z$ 为集中于点 $z$ 的概率分布,向原始分布 $F$ 中加入比例为 $\varepsilon$ 的点污染,得到:
\[F_{\varepsilon,z} =(1-\varepsilon)F+\varepsilon\Delta_z.\]估计泛函 $T$ 在分布 $F$ 处的影响函数定义为:
\[\boxed{ \operatorname{IF}(z;T,F) =\lim_{\varepsilon\downarrow0} \frac{T(F_{\varepsilon,z})-T(F)}{\varepsilon} }\]它回答的问题是:若有一个比例趋近于零的污染成分出现在位置 $z$,估计目标将发生怎样的一阶变化?这是局部敏感性概念,不等同于有限比例污染下的整体稳健性,也不能单独替代崩溃点(Breakdown Point)分析。
M-估计的影响函数推导
对于固定尺度 M-估计,令 $\theta_0=T(F)$,总体方程可表示为:
\[\int \psi\left(\frac{x-\theta_0}{s}\right) \,dF(x)=0.\]把 $F$ 换成污染后的 $F_{\varepsilon,z}$,再对 $\varepsilon$ 进行一阶微分。利用总体方程在 $\varepsilon=0$ 时等于零,并假设允许交换微分与积分,得到:
\[\boxed{ \operatorname{IF}(z;T,F) =s\frac{ \psi\left(\frac{z-\theta_0}{s}\right) }{ \mathbb E_F\left[ \psi'\left(\frac{X-\theta_0}{s}\right) \right] } }\]这一表达式将估计量的局部稳健性直接连接到得分函数。对于均值估计,$\psi(u)=u$ 且分母为一,因此 $\operatorname{IF}(z;T,F)=z-\theta_0$,当污染点远离中心时,影响函数无界。对于 Huber 估计,只要 $A=\mathbb E_F[\psi_c’(Z)]>0$,就有以下上界:
\[\sup_{z\in\mathbb R} \lvert\operatorname{IF}(z;T,F)\rvert \le\frac{sc}{A}<\infty.\]因而 Huber 估计具有有界的局部影响函数。这里最值得记住的并非公式本身,而是其中的结构性结论:在常规正则条件下,有界得分函数使 M-估计量的单点污染影响受到控制。对于回归问题,这一结论还必须考虑解释变量的杠杆效应,不能把一维位置模型的结论直接照搬到所有回归场景。
影响函数与渐近方差是同一理论的两面
若 M-估计量具有渐近线性表示,则可以写为:
\[\sqrt n(\hat\theta_n-\theta_0) =\frac{1}{\sqrt n} \sum_{i=1}^{n} \operatorname{IF}(X_i;T,F)+o_p(1).\]因此,渐近方差可以由影响函数的二阶矩给出:
\[\operatorname{Var}(\hat\theta_n) \approx\frac{1}{n} \mathbb E_F\left[ \operatorname{IF}(X;T,F)^2 \right].\]代入前面得到的影响函数表达式,正好恢复 M-估计的夹心方差公式。换句话说,影响函数描述的是单个微小污染源的作用,而渐近方差描述的是全部随机观测的影响在大样本下如何累积;二者建立在同一个一阶线性化结构上。
Huber(1964)的关键贡献:污染模型与极小极大估计
为什么要引入污染分布类
经典正态模型假定全部观测都来自同一个指定分布。Huber 的思路不是去猜测异常值究竟服从哪一种分布,而是用一个包含多种可能分布的邻域来刻画模型不确定性。对于标准正态分布 $\Phi$,考虑污染模型:
\[\boxed{ F=(1-\varepsilon)\Phi+\varepsilon H }\]其中 $0\le\varepsilon<1$ 为污染比例,$H$ 是未知的污染分布。记 $\mathcal P(\mathbb R)$ 为实数轴上的所有概率分布所构成的集合,对应的分布类定义为:
\[\mathcal F_{\varepsilon} =\left\{ (1-\varepsilon)\Phi+\varepsilon H: H\in\mathcal P(\mathbb R) \right\}.\]当 $\varepsilon=0$ 时,模型回到理想正态分布;当 $\varepsilon>0$ 时,部分观测允许来自其他分布。这样一来,估计方法就不能只依据正态模型下的效率评价,还需要考察其在所允许的分布偏离范围内最差能表现到什么程度。
极小极大渐近方差问题
Huber 在位置估计及相应的正则性、对称性与等变性条件下,研究了污染分布类上的极小极大渐近方差。以得分函数 $\psi$ 表示估计方案,问题的基本结构可以概括为:
\[\boxed{ \inf_{\psi\in\Psi} \sup_{F\in\mathcal F_{\varepsilon}} V(\psi,F) }\]这里 $\Psi$ 是允许的得分函数类,$V(\psi,F)$ 是估计量在分布 $F$ 下的渐近方差系数。固定单位尺度时,该系数具有如下形式:
\[V(\psi,F) =\frac{ \int\psi(x)^2\,dF(x) }{ \left(\int\psi'(x)\,dF(x)\right)^2 }.\]该表达式默认总体得分已中心化并满足相应的正则条件。尤其要注意,任意非对称污染可能使估计泛函偏离原始正态位置参数,因此固定污染比例下的偏差问题与围绕估计泛函计算渐近方差的问题不能混为一谈。Huber 原论文分别研究了污染影响、最坏情况渐近方差及相关的偏差现象。
最不利分布为什么导出 Huber 得分
在正态污染邻域的极小极大构造中,存在一个关键的最不利分布(Least Favorable Distribution)。在单位尺度、对称情形下,其密度可以写为:
\[f_0(x)= \begin{cases} (1-\varepsilon)\phi(x), & \lvert x\rvert\le c,\\ (1-\varepsilon)\phi(c) \exp\bigl[-c(\lvert x\rvert-c)\bigr], & \lvert x\rvert>c. \end{cases}\]中心区域保留高斯密度形状,尾部则转为指数衰减。为了使 $f_0$ 的总积分等于一,污染比例 $\varepsilon$ 与截断点 $c$ 必须满足:
\[\frac{1}{1-\varepsilon} =2\Phi(c)-1+\frac{2\phi(c)}{c}.\]这个归一化条件说明,在极小极大问题中,最优截断阈值并不是独立于污染比例随意给定的。若允许的污染比例趋于零,对应的截断点趋于无穷,从而逐渐恢复正态模型下的线性得分;允许更大的污染比例,则需要更早截断极端残差的贡献。
对最不利密度取负对数导数,有:
\[-\frac{d}{dx}\log f_0(x) =\begin{cases} -c, & x<-c,\\ x, & -c\le x\le c,\\ c, & x>c. \end{cases}\]右侧恰好是 Huber 得分函数 $\psi_c(x)$。对它积分,得到的损失函数与 Huber 损失只差一个与位置参数无关的常数。这是 Huber(1964)理论中的核心联系:截断线性的得分函数并不只是为了方便处理异常值而提出的启发式设计,它在指定污染邻域和极小极大渐近方差准则下具有严格的统计决策意义。
还应区分两类常见的阈值选择方法:一类根据污染比例求得极小极大截断点,另一类根据理想正态模型下的渐近相对效率选择截断点。常用的 $c\approx1.345$ 属于后者,即约 $95\%$ 正态渐近效率的选择,不能直接解释为适用于任何问题的最优污染比例。
数值计算:M-估计如何求解
从估计方程导出迭代重加权最小二乘
对于一维 Huber 位置估计,得分函数连续且单调,因此可以直接对估计方程使用二分法或 Brent 方法求根。在更一般的问题中,迭代重加权最小二乘(Iteratively Reweighted Least Squares,IRLS)提供了一种直观的计算方案。
当 $u\ne0$ 时定义权重函数 $w(u)=\psi(u)/u$,并在零点按连续极限赋值。Huber 得分对应的权重为:
\[w_c(u)= \begin{cases} 1, & \lvert u\rvert\le c,\\ \frac{c}{\lvert u\rvert}, & \lvert u\rvert>c. \end{cases}\]由恒等式 $\psi(u)=w(u)u$,原来的估计方程可以改写为加权残差方程:
\[\sum_{i=1}^{n} w_c\left(\frac{X_i-\theta}{s}\right) (X_i-\theta)=0.\]若暂时固定由当前估计 $\theta^{(t)}$ 计算的权重 $w_i^{(t)}=w_c((X_i-\theta^{(t)})/s)$,就得到如下更新公式:
\[\boxed{ \theta^{(t+1)} =\frac{ \sum_{i=1}^{n}w_i^{(t)}X_i }{ \sum_{i=1}^{n}w_i^{(t)} } }\]这一形式说明,Huber 估计可以解释为随残差大小自适应调整权重的估计:中心观测保留完整权重,大残差观测的权重按残差幅度的倒数衰减。对于非凸稳健损失,类似的迭代方案可能出现多个驻点,因此初始化与收敛诊断更重要;而 Huber 损失为凸函数,一维位置问题的求解相对稳定。
一个能够手工验证的异常值案例
考虑样本 $x=(0,0,1,1,2,20)$。其中前五个数集中在零到二之间,最后一个观测明显远离主体数据。样本均值为 $4$,中位数为 $1$。现在令固定尺度 $s=1$,Huber 阈值 $c=1.5$,求对应的 Huber M-估计量。
假设解附近的前五个残差都位于二次区域,而异常值 $20$ 对应的得分已经达到上界 $1.5$。估计方程便成为:
\[(0-\theta)+(0-\theta) +(1-\theta)+(1-\theta) +(2-\theta)+1.5=0.\]整理得到 $5.5-5\theta=0$,从而 $\hat\theta_{\mathrm{Huber}}=1.1$。将它代回得分函数,可以验证前五个标准化残差的绝对值都不超过 $1.5$,而最后一个残差超过该阈值,所以分段假设是自洽的。三种估计量分别为:均值 $4$、中位数 $1$、Huber 估计 $1.1$。这一比较清楚地表明,Huber 估计既不会完全忽略异常观测,也不会让异常观测的巨大数值直接主导估计。
下面的 Python 代码使用 SciPy 中的 Brent 求根算法计算位置 M-估计量。对于一维凸 Huber 损失,这种做法不需要手工推导 IRLS 的收敛控制,也能直接检查估计方程的零点:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
from scipy.optimize import brentq
def huber_location(x, scale=1.0, c=1.5):
x = np.asarray(x, dtype=float)
if x.size == 0:
raise ValueError("The sample must not be empty")
if scale <= 0 or c <= 0:
raise ValueError("scale and c must be positive")
def score(theta):
residuals = (x - theta) / scale
return np.clip(residuals, -c, c).sum()
lower = float(x.min() - c * scale)
upper = float(x.max() + c * scale)
return brentq(score, lower, upper)
x = np.array([0, 0, 1, 1, 2, 20], dtype=float)
print("Mean:", np.mean(x))
print("Median:", np.median(x))
print("Huber:", huber_location(x))
这段代码的预期输出是 Mean: 4.0、Median: 1.0 和 Huber: 1.1。如果将最后一个观测进一步增大,样本均值还会继续变化,而在其余观测仍属于相同得分区域时,Huber 估计的解可以保持不变。这正是得分截断在有限样本估计方程中的直接体现。
实际使用中的尺度估计与适用边界
尺度参数为何会改变截断行为
在标准化残差 $u=(X-\theta)/s$ 中,阈值 $c$ 是无量纲的;换回原始观测单位后,截断发生在残差幅度约为 $cs$ 的地方。因此,如果 $s$ 的选择不合理,Huber 损失在实际数据上的截断强度也会失去预期意义。经典的稳健尺度估计方法之一是中位数绝对偏差(Median Absolute Deviation,MAD):
\[\operatorname{MAD} =\operatorname{median}_i \left( \lvert X_i-\operatorname{median}_jX_j\rvert \right).\]在正态模型下,可以用 $\hat s\approx1.4826\operatorname{MAD}$ 获得经过一致性校正的尺度估计。但尺度估计并非完全独立于位置估计:若两者由同一样本联合估计,需要相应的联合估计方程与渐近协方差分析。在某些对称性条件下,尺度估计的误差对位置估计的一阶影响可以消失,但这属于需要证明的模型性质,不能无条件假定。
当大量观测完全相同时,MAD 可能为零,使标准化残差无法正常计算。实际实现应检查这种退化情形,选择合适的替代尺度估计或预先说明对退化数据的处理规则,避免直接套用正常尺度下的公式。
有界影响函数并不等于对一切异常机制都稳健
Huber 得分函数在尾部被截断为常数,而不是随残差增大重新趋近于零,因此极端观测仍然保留有限但非零的得分贡献。某些重降型得分函数(Redescending Score)可以进一步削弱极端残差的贡献,但相应损失函数往往是非凸的,可能带来多个局部极小点和更困难的理论分析。
此外,一维位置估计中的有界影响函数结论,不能直接推出稳健回归可以自动抵御高杠杆解释变量。稳健回归除了考虑响应变量方向上的大残差,还需要分析设计矩阵的异常结构。影响函数、崩溃点、污染偏差和最坏情况方差描述的也是不同层次的稳健性,选择估计方法时必须首先明确数据可能以何种方式偏离理想模型。
总结:如何形成对 M-估计的完整理解
M-估计理论可以沿着一条清晰的逻辑链条来掌握:先指定损失函数 $\rho$,由其导数构造得分函数 $\psi$;再利用样本估计方程得到参数估计量,利用总体估计方程确定估计的统计目标;随后通过一阶展开获得渐近正态性与夹心方差,通过影响函数分析微小污染对估计结果的作用,最后在污染分布类上研究估计方案的最坏情况表现。
这条链条中的几个核心关系可概括为:
\[\rho \;\xrightarrow{\text{differentiate}}\; \psi \;\xrightarrow{\text{construct}}\; \text{estimating equation} \;\xrightarrow{\text{linearize}}\; \text{asymptotic distribution}.\]而得分函数同时决定影响函数的形式,因此也连接了统计效率与局部稳健性。Huber(1964)的深刻之处在于,他将得分截断放入一个可证明的极小极大统计决策问题中:在特定正态污染邻域中,中心线性、尾部截断的得分函数不仅具有直观的异常值防御作用,还与最不利分布及最坏情况渐近方差相联系。理解了这一结构,再去研究稳健线性回归、联合位置尺度估计、Z-估计和广义 M-估计,就有了统一而可复用的理论基础。
参考文献
Huber, P. J. (1964). Robust Estimation of a Location Parameter. The Annals of Mathematical Statistics, 35(1), 73–101. 原文第 76–78 页讨论凸损失 M-估计的渐近正态性与渐近方差,第 80–92 页展开污染模型、极小极大理论与最不利分布,第 93–100 页进一步讨论尺度估计和计算问题。
Huber, P. J., & Ronchetti, E. M. (2009). Robust Statistics (2nd ed.). Wiley. 适合在理解原始论文后系统学习稳健统计理论,包括位置估计、回归估计和渐近推断。
Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J., & Stahel, W. A. (1986). Robust Statistics: The Approach Based on Influence Functions. Wiley. 适合深入理解影响函数、局部稳健性与不同稳健性准则之间的关系。
van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press. 可用于补充一般 M-估计和 Z-估计的一致性、渐近正态性与经验过程理论。
本文关于凸损失估计、渐近方差及污染模型下极小极大估计的主要理论脉络来自 Huber(1964);关于影响函数的统一表述、夹心方差术语与 Python 实现,则是为教学目的加入的后续理论及计算方法。
