§2. M-оценки

Усечённое среднее $\overline X_\alpha$ служит компромиссом между средним $\overline X$ ($\alpha = 0$) и выборочной медианой $MED$ ($\alpha \to \frac12$). Другой компромисс подсказан тем фактом, что среднее и медиана суть величины, минимизирующие, соответственно, $\sum(X_i - \theta)^2$ и $\sum|X_i - \theta|$ (см. задачу 1 гл. 1 и задачу 4 гл. 7). П. Хьюбер предложил минимизировать меру «разброса»

где функция $\rho_b(x)$ при некотором $b > 0$ задаётся формулой

Она пропорциональна $x^2$ при $|x| \le b$, а вне отрезка $[-b, b]$ заменяет ветви параболы на две прямые линии. Эти куски графика согласованы так, что $\rho_b(x)$ и её производная $\psi_b(x) = \rho_b'(x)$ непрерывны. При увеличении $b$ функция $\rho_b$ будет совпадать с $x^2/2$ на всё большей части её определения, так что оценка Хьюбера будет сближаться с $\overline X$. В свою очередь, при уменьшении $b$ происходит сближение с выборочной медианой. (В качестве компромисса иногда предлагается значение $b = 1{,}5$.)

По кирпичикам: в точке $|x| = b$ парабола $x^2/2$ имеет значение $b^2/2$ и наклон $b$; прямая $b|x| - b^2/2$ имеет там то же значение $b^2 - b^2/2 = b^2/2$ и тот же наклон $b$. Поэтому шов гладкий. Производная $\psi_b(x) = \max\bigl(-b, \min(b, x)\bigr)$ — это $x$, «обрезанный» на уровне $\pm b$.

Смысл по-пацански: среднее штрафует за отклонение квадратом — далёкий чудак получает штраф до небес и перетягивает центр на себя. Медиана штрафует модулем — всем одинаково «по рублю за метр», но и близких она слушает плохо. Хьюбер: близких — по квадрату (точно), далёких — по модулю (не даём перетягивать).

Подобно тому, как усечённые средние $\overline X_\alpha$ ($0 < \alpha < \frac12$) образуют однопараметрическое подмножество класса $L$-оценок, оценки Хьюбера ($0 < b < \infty$) образуют однопараметрическое подмножество класса так называемых $M$-оценок.

Определение и асимптотика

$M$-оценка $M_n$ параметра сдвига $\theta$ определяется для некоторой функции $\rho(x)$ как точка минимума по $\theta$ функции

Если $\rho$ — выпуклая и чётная, то можно показать, что значения $\theta$, минимизирующие сумму (8), образуют некоторый отрезок, а если $\rho$ строго выпукла (П4), то минимизирующее значение единственно.

Если функция $\rho$ имеет производную $\psi = \rho'$, то оценка $M_n$ является решением уравнения

Пусть элементы выборки $X_i$ имеют функцию распределения $F(x - \theta)$, где $F$ обладает чётной плотностью $p(x)$. Тогда при слабых предположениях относительно $\psi$ и $F$

где асимптотическая дисперсия $\sigma^2(F, \psi)$ задаётся формулой

(Доказательство (10) и подробное изложение теории $M$-оценок есть в книге Хьюбера «Робастность в статистике».)

По кирпичикам — откуда (11): разложим $\sum\psi(X_i - M_n) = 0$ около истинного $\theta$: $\sum\psi(X_i - \theta) - (M_n - \theta)\sum\psi'(X_i - \theta) \approx 0$. Отсюда $\sqrt n(M_n - \theta) \approx \frac{n^{-1/2}\sum\psi(X_i - \theta)}{n^{-1}\sum\psi'(X_i - \theta)}$. Числитель по ЦПТ $\approx \mathcal N(0, \mathsf M\psi^2)$, знаменатель по ЗБЧ $\to \mathsf M\psi'$. Получаем дисперсию $\mathsf M\psi^2/(\mathsf M\psi')^2$.

Если $\rho(x) = -\ln p(x)$, то $\psi = \rho' = \gamma$ (см. формулу (4) §1). В этом случае минимизация суммы (8) эквивалентна максимизации произведения $\prod p(X_i - \theta)$, т. е. $M$-оценка параметра сдвига совпадает с оценкой максимального правдоподобия (ОМП) (см. §4 гл. 9). Таким образом, $M$-оценки обобщают ОМП (ей они обязаны своим названием). Известно, что при выполнении условий регулярности (см. §3 гл. 9) ОМП является асимптотически эффективной оценкой с $\sigma^2(F, \gamma) = 1/I(F)$, где $I(F)$ определена в (5).

Смысл по-пацански: «M» — от «maximum likelihood». Взял штраф $\rho = -\ln p$ — получил ОМП. Взял любой другой разумный штраф — получил родственника ОМП, может, чуть менее точного, зато не такого доверчивого.

Пример 2. Толерантность M-оценок

Для закона $\mathcal N(0, 1)$ функция $\gamma(x) = x$, поэтому оптимальная оценка $M_n$ совпадает с $\overline X$. Нулевая толерантность этой оценки (см. §4 гл. 8) связана с неограниченностью функции $\gamma(x)$.

Для логистического распределения с $F(x) = 1/(1 + e^{-x})$ имеем $\gamma(x) = \operatorname{th}(x/2)$. В силу строгой монотонности и ограниченности $\gamma(x)$ толерантность решения уравнения (9) равна $\frac12$.

В случае закона Лапласа с плотностью $p(x) = \frac12e^{-|x|}$ весовая функция $\rho(x) = |x|$ не имеет производной в нуле. Полагая $\gamma(0) = 0$, получаем $\gamma(x) = \operatorname{sign}x$. Соответствующая ОМП — выборочная медиана $MED$, имеющая $\sigma^2(F, \gamma) = 1$ и обладающая асимптотической толерантностью $\frac12$.

Для распределения Коши с плотностью $p(x) = 1/[\pi(1 + x^2)]$ функция $\gamma(x) = 2x/(1 + x^2)$ оказывается немонотонной: растёт до $x = 1$, потом убывает к нулю. Поэтому уравнение (9) может иметь несколько решений (см. вопрос 3 в гл. 9). Коллинз (1976) и Кларк (1984) показали, что если а) брать ближайшее к выборочной медиане решение или б) пользоваться методом Ньютона, выбирая в качестве начального приближения $MED$ (см. теорему 4 гл. 9), то оценка будет асимптотически эффективной, причём от медианы будет унаследована толерантность $\frac12$.

Наконец, рассмотрим случай «наименее благоприятного распределения» Хьюбера с плотностью $p(x)$, пропорциональной функции $\exp\{-\rho_b(x)\}$, где $\rho_b(x)$ задаётся формулой (7). Для этого закона $\gamma(x)$ совпадает с $\psi_b(x) = \rho_b'(x)$. Толерантность оценки Хьюбера также равна $\frac12$.

По кирпичикам — общее правило: ограниченная $\psi$ даёт положительную толерантность (одна точка не может сдвинуть сумму $\sum\psi$ больше чем на $\sup|\psi|$); неограниченная — нулевую.

Пример 3. Минимаксный подход Хьюбера

Вместо того, чтобы предполагать функцию распределения $F$ известной, более реалистично допустить, что она известна лишь приближённо. Одной из таких моделей, предложенных Хьюбером, является представление $F$ в виде смеси

где функция распределения $G(x)$ и $0 \le \varepsilon \le 1$ заданы, а $H(x)$ — это произвольная неизвестная функция распределения (предполагается только, что $G$ и $H$ симметричны относительно нуля, см. гл. 8).

Обозначим через $\Omega(G, \varepsilon)$ семейство распределений, удовлетворяющих (12) для фиксированных $G$ и $\varepsilon$. (Легко видеть, что любое распределение вида $(1 - \varepsilon')G(x) + \varepsilon'H(x)$, где $0 \le \varepsilon' \le \varepsilon$, также входит в $\Omega(G, \varepsilon)$.) В соответствии с минимаксным подходом представляет интерес нахождение $M$-оценки (функции $\psi$), на которой достигается

В частном случае, когда $G$ совпадает с функцией распределения стандартного нормального закона $\Phi$ с плотностью $\varphi(x) = (2\pi)^{-1/2}e^{-x^2/2}$, решением этой задачи служит оценка Хьюбера с $\psi_b(x)$, у которой $b = b(\varepsilon)$ определяется из соотношения

Некоторые значения $b(\varepsilon)$ приведены в таблице (из книги Хьюбера):

В частности, видим, что 5%-му «загрязнению» отвечает $b \approx 1{,}4$. (Решая (14) численно, мы получаем 2,633 и 1,398 для $\varepsilon = 0{,}001$ и $0{,}05$; остальные значения совпадают — расхождение в третьем знаке, видимо, от округлений в исходной таблице.)

По ходу доказательства устанавливается, что при $\psi = \psi_b$ наибольшее значение $\sigma^2(F, \psi_b)$ достигается на распределении $F$ с плотностью $p(x)$, пропорциональной $\exp\{-\rho_b(x)\}$, ввиду чего его и называют «наименее благоприятным распределением» Хьюбера.

По кирпичикам: правая часть (14) — это полная масса плотности $\varphi(x)$ на $[-b, b]$ плюс масса двух экспоненциальных «хвостов» $\varphi(b)e^{-b(|x|-b)}$ снаружи. Левая часть $1/(1-\varepsilon)$ — во сколько раз должна быть больше общая масса, чтобы после умножения на $(1 - \varepsilon)$ получилась вероятностная плотность. Чем больше $\varepsilon$, тем тяжелее хвосты и тем меньше $b$.

Смысл по-пацански: минимакс — это стратегия параноика. «Пусть 5 % данных подсунул злой гений, причём самым вредным образом. Какую оценку взять, чтобы даже в худшем случае потерять поменьше?» Ответ Хьюбера: обрезай вклад на уровне $b \approx 1{,}4$.

W-оценки Тьюки

Замечание 2. Дж. Тьюки в 1970 г. предложил для параметра сдвига $\theta$ так называемую $W$-оценку $\tilde\theta$, представляющую собой решение уравнения

где $w(u)$ — некоторая весовая функция. Правая часть (15) является взвешенным средним наблюдений $X_i$ с весами $w_i = w(X_i - \tilde\theta)$, $i = 1, \dots, n$.

Для подсчёта $\tilde\theta$ обычно прибегают к использованию итерационного алгоритма. Начиная с $\tilde\theta_0 = MED$, вычисляются величины

до тех пор, пока не будет достигнута стабилизация $\tilde\theta_j$. (Процедура итерационного перевзвешивания часто применяется в регрессионном анализе — см., например, алгоритм LOWESS в §2 гл. 22.)

Отметим, что $W$-оценки являются на самом деле $M$-оценками специального вида: в силу (15) $\sum_{i=1}^n(X_i - \tilde\theta)w(X_i - \tilde\theta) = 0$, и значит, $\sum_{i=1}^n\psi_w(X_i - \tilde\theta) = 0$ с $\psi_w(u) = uw(u)$.

Из-за простоты вычисления популярны одношаговые $W$-оценки $\tilde\theta_1$ в (16), имеющие очень хорошие характеристики: они наследуют толерантность $\frac12$ выборочной медианы и обладают асимптотической дисперсией $M$-оценки с $\psi = \psi_w$, задаваемой формулой (11).

По кирпичикам: оценка Хьюбера — тоже $W$-оценка, с весом $w(u) = \psi_b(u)/u = \min(1, b/|u|)$: близким точкам вес 1, далёким — обратно пропорционально расстоянию. Популярный «бивес» Тьюки $w(u) = (1 - (u/c)^2)^2$ при $|u| < c$ и 0 снаружи вообще обнуляет далёкие точки; при $c \approx 4{,}685$ (в единицах масштаба) его эффективность для нормального закона — 95 %.

Смысл по-пацански: $W$-оценка — это «посчитай среднее, но сначала понизь голос тем, кто далеко от предыдущего ответа». Повтори пару раз. Начинаешь с медианы — значит, выбросы с самого начала не задают тон.

Масштаб

Замечание 3. Если для описания данных используется не модель сдвига, а модель сдвига-масштаба, то прежде чем оценивать сдвиг, надо оценить параметр масштаба с помощью какой-нибудь робастной оценки, например, (нормированной) медианы абсолютных отклонений

где $\Phi^{-1}(x)$ — функция, обратная к функции распределения $\mathcal N(0, 1)$, $1/\Phi^{-1}(3/4) \approx 1{,}483$ (она ранее встречалась в §1 гл. 19). Толерантность $MAD$ равна $\frac12$.

Затем для некоторой $\psi$ оценка параметра сдвига $\hat\theta$ находится как решение уравнения

С другой стороны, можно также совместно оценить параметры сдвига и масштаба при помощи обобщения $M$-оценок на случай векторного параметра.

По кирпичикам: зачем делить на масштаб? Порог $b = 1{,}5$ — это «полтора чего?». Если данные в миллиметрах, а не в метрах, без масштаба оценка Хьюбера превратилась бы в среднее (все отклонения $\ll b$) или в медиану (все $\gg b$). Нормировка на $MAD$ делает $b$ безразмерным: «полтора типичных отклонения». А $1{,}483$ подобран так, чтобы для нормальных данных $MAD$ оценивал именно $\sigma$.

В виджете — выборка с подвижным выбросом и итерации (16) для Хьюбера и Тьюки, масштаб по $MAD$:

Проверь себя