§2. Непараметрическая регрессия

Изучаемая в этом параграфе задача заключается в сглаживании данных $\{(x_i, y_i)\}_{i=1}^n$ при помощи некоторой непрерывной кривой $y = m(x)$. Данные рассматриваются как реализация выборки, описываемой одной из двух моделей.

Смысл по-пацански: в гл. 21 мы заранее говорили: «зависимость — прямая» (или парабола) — и подгоняли коэффициенты. Здесь формулу не заказываем вообще. Пусть точки сами нарисуют кривую, а мы только аккуратно сгладим шум.

Модель со случайным планом эксперимента

Пары $(X_i, Y_i)$ считаются независимыми и одинаково распределёнными. Кривая регрессии определяется как (см. П7)

Это определение корректно, если $\mathsf M|Y| < \infty$. Если у вектора $(X, Y)$ существует плотность $p(x, y)$, то

где $p(x) = \int p(x, y)\,dy$ — маргинальная плотность случайной величины $X$.

Например, пусть $(X, Y) \sim \mathcal N(\boldsymbol\mu, \boldsymbol\Sigma)$ (см. П9), где $\boldsymbol\mu = (\mu_1, \mu_2)$, $\boldsymbol\Sigma = \|\sigma_{ij}\|_{2\times2}$. Несложно подсчитать, что в этом случае

т. е. кривая регрессии линейна.

По кирпичикам: условное распределение $Y$ при $X = x$ в нормальном случае снова нормально, его среднее сдвигается от $\mu_2$ пропорционально отклонению $x - \mu_1$, а коэффициент — ковариация, делённая на дисперсию предиктора $X$, т. е. $\sigma_{12}/\sigma_{11}$. (У Лагутина в знаменателе напечатано $\sigma_{22}$ — это опечатка: дисперсия отклика тут ни при чём, ср. наклон прямой МНК $b = \sigma_{XY}/\sigma_X^2$ из гл. 21.)

Смысл по-пацански: для нормального облака «кривая» регрессии — честная прямая. Всё нелинейное начинается, когда облако не нормальное, — вот для этого и нужна непараметрика.

Модель с фиксированным планом эксперимента

В случае управляемой неслучайной предикторной (предсказывающей) переменной $X$ применяется модель

где $x_i$ — «узлы» (задаваемые значения $X$, в которых производятся измерения отклика $Y$), $\varepsilon_i$ — независимые и одинаково распределённые случайные ошибки, $\mathsf M\varepsilon_i = 0$, $\mathsf D\varepsilon_i = \sigma^2$.

Замечание 1. Если допустимы многократные наблюдения в фиксированной точке $X = x$, то можно оценить $m(x)$ с помощью среднего арифметического соответствующих $Y_i$.

Однако часто этому препятствуют финансовые ограничения: слишком дорого проводить более одного измерения для каждого фиксированного уровня $x$ предиктора $X$. Возможно также, что условия эксперимента нельзя воспроизвести из-за разрушения объекта (например, при изучении безопасности водителя при столкновении автомобилей) и т. п.

Локальное усреднение

Если есть уверенность, что $m$ — гладкая кривая, наблюдения $X_i$ вблизи точки $x$ должны содержать информацию о значении $m$ в точке $x$. Таким образом, представляется возможным использовать нечто вроде локального усреднения близких к $x$ данных для формирования оценки $m(x)$.

— Р. Эйбанк

В основе процедуры оценивания кривой регрессии $m(x)$ лежит идея локального усреднения. Оценка определяется как взвешенное среднее

где веса $w_i(x)$ велики для $X_i$, близких к точке $x$, и малы для остальных $X_i$.

Смысл по-пацански: хочешь узнать, сколько в среднем получают на твоей улице, — спрашивай соседей, а не весь город. Ближнего соседа слушай внимательно, дальнего — вполуха, а с другого конца города — вообще не слушай.

Замечание 2. В соответствии с примером 5 гл. 21 величина $\hat m(x)$ при каждом $x$ является оценкой взвешенных наименьших квадратов, так как она служит решением задачи минимизации функции

По кирпичикам: производная по $\theta$ равна $-2\sum w_i(Y_i - \theta)$; приравняв её нулю, получим $\theta = \sum w_iY_i/\sum w_i$ — это и есть (7). То есть в каждой точке $x$ мы подгоняем константу взвешенным МНК.

Оценка Надарая — Ватсона

Один из наиболее простых способов определения весов $w_i(x)$ опирается на использование некоторого ядра $q(y)$ (см. §1):

Здесь $h = h_n$ обозначает ширину «окна сглаживания». Для многомерной предикторной переменной $\boldsymbol X_i = (X_{i1}, \dots, X_{im})$ можно взять

Задавать веса можно и по-другому. Для одномерной предикторной переменной Гассер и Мюллер в 1979 г. предложили

где $X_{(1)} \le X_{(2)} \le \dots \le X_{(n)}$ — упорядоченные по возрастанию $X_i$, а отрезки $[s_{i-1}, s_i]$ разбивают прямую на $n$ кусков, по одному на каждое наблюдение: $s_0 = -\infty$, $s_n = +\infty$, а $s_i = \bigl(X_{(i)} + X_{(i+1)}\bigr)/2$ — середины между соседними узлами. Заметим, что $\sum_{i=1}^n\tilde w_i(x) = 1$ (интегралы по кускам складываются в $\int q_h = 1$), поэтому делить на сумму весов в (7) уже не нужно.

(У Лагутина пределы интегрирования записаны как $X_{(i-1)}$ и $X_{(i)}$ с суммой по $i$ до $n + 1$ — так весов получается на один больше, чем наблюдений. Мы приводим стандартную форму Гассера — Мюллера с серединами; именно она реализована в gasserMuller из нашей библиотеки.)

Смысл по-пацански: Гассер и Мюллер каждому наблюдению выдают свой «участок» на оси и взвешивают его по тому, какая часть ядра на этот участок легла. Где точки густо — участки маленькие, и вес каждой точки меньше; где редко — участок большой, вес больше. Честная делёжка.

В случае весов $w_i(x)$ из формулы (9) правая часть равенства (7) приобретает вид

Выражение (11) определяет так называемую оценку Надарая — Ватсона (Nadaraya, 1964; Watson, 1964). С учётом формулы (3) §1

где $\hat p_n(x)$ — ядерная оценка Розенблатта — Парзена маргинальной плотности $p(x)$ переменной $X$.

По кирпичикам: знаменатель в (11), делённый на $n$, — это ровно оценка плотности из §1. Так что Надарая — Ватсон — это «оценка $\int y\,p(x, y)\,dy$, делённая на оценку $p(x)$», т. е. подстановка ядерных оценок в формулу для $m(x)$ из модели со случайным планом.

Состоятельность и риск

В следующей теореме приведены условия, обеспечивающие состоятельность оценки Надарая — Ватсона $\hat m(x)$.

Теорема 2. Пусть для модели со случайным планом эксперимента и одномерной предикторной переменной

1. $\int|q(y)|\,dy < \infty$,
2. $yq(y) \to 0$ при $|y| \to \infty$,
3. $\mathsf MY^2 < \infty$,
4. $n \to \infty$, $h_n \to 0$, $nh_n \to \infty$.

Тогда $\hat m(x) \xrightarrow{\mathsf P} m(x)$ в любой точке непрерывности функций $m(x)$, $p(x)$ и $\sigma^2(x) = \mathsf D(Y \mid X = x)$ такой, что $p(x) > 0$.

(Доказательство этой теоремы приведено в книге В. Хардле «Прикладная непараметрическая регрессия».)

По кирпичикам: условие $h_n \to 0$ убивает смещение (окно стягивается к точке $x$), а $nh_n \to \infty$ — дисперсию (в окно попадает всё больше точек: их примерно $nh_np(x)$). Отсюда и требование $p(x) > 0$: где точек нет совсем, усреднять нечего.

Сформулированный результат показывает, что оценка $\hat m(x)$ из формулы (11) сходится по вероятности к $m(x)$ — истинной кривой регрессии. Естественно поставить вопрос, какова скорость этой сходимости. Одной из возможных мер точности оценки является квадратичный риск $R(x, h) = \mathsf M[\hat m(x) - m(x)]^2$.

Теорема 3 (Гассер и Мюллер, 1984). Рассмотрим модель с фиксированным планом эксперимента для одномерной предикторной переменной $X$. (Для случайного плана скорость сходимости та же, но с другими константами.) Без ограничения общности можно считать, что «узлы» $x_i \in [0, 1]$. Предположим, что

1. используется ядро $q(y)$ с $\alpha = \int q^2(y)\,dy < \infty$ и $\beta = \int y^2q(y)\,dy < \infty$;
2. $q(y)$ имеет носитель $(-1, 1)$, причём $q(-1) = q(1) = 0$;
3. веса $\tilde w_i(x)$ задаются формулой (10);
4. $m(x)$ дважды непрерывно дифференцируема;
5. $\max_{1\le i\le n+1}(x_{(i)} - x_{(i-1)}) = O(1/n)$ $\;(x_{(0)} = 0,\; x_{(n+1)} = 1)$;
6. $\mathsf D\varepsilon_i = \sigma^2$, $i = 1, \dots, n$;
7. $n \to \infty$, $h_n \to 0$, $nh_n \to \infty$.

Тогда

Видим, что так же, как и в формуле (5) §1, квадратичный риск состоит из двух частей: квадрата смещения и дисперсии, где смещение (как функция от $h_n$) возрастает, а дисперсия — убывает. В. Хардле пишет: «Это качественное соображение раскрывает сущность задачи сглаживания: баланс между дисперсией и квадратом смещения».

Смысл по-пацански: та же песня, что и с плотностью. Смещение растёт как $h^4$ и сильнее там, где кривая круче изгибается ($m''$ большая): широкое окно «срезает» горбы. Дисперсия падает как $1/(nh)$. Оптимум снова $h \sim n^{-1/5}$, ошибка $\sim n^{-2/5}$.

Пример 2. Обнаружение локального максимума кривой регрессии

На рис. 4 у Лагутина представлены моделированные данные. Использовалась модель с фиксированным планом эксперимента: $n = 50$, $x_i = i/n$, $i = 1, \dots, n$. Точки порождены «зашумлением» с помощью взятой из таблицы T1 реализации $\varepsilon_i \sim \mathcal N(0, 1)$ кривой регрессии

(т. е. на линейную функцию $y = 1 - x$ накладывается нормальный «горб» с центром в $1/2$ и шириной примерно равной $2\sigma = 1/10$).

Из-за большой интенсивности «шума» локальный максимум $m(x)$ практически не различим. Однако он хорошо заметен на графике оценки $\hat m(x)$ из (11), для которой использовалось гладкое квартическое ядро $q(y) = (15/16)(1 - y^2)^2I_{\{|y|\le1\}}$ и был определён методом пропуска (см. ниже) оптимальный размер «окна сглаживания» $h = 0{,}09$.

По кирпичикам: $\exp\{-200(x - 1/2)^2\} = \exp\{-(x - 1/2)^2/(2\sigma^2)\}$ при $2\sigma^2 = 1/200$, т. е. $\sigma = 0{,}05$, $2\sigma = 0{,}1$. Высота горба — 1, а стандартное отклонение шума — тоже 1. Горб «утоплен» в шуме по самую макушку.

Смысл по-пацански: на глаз в облаке точек горба не видно — шум того же роста, что и горб. А сглаживание вытаскивает его наружу: соседи усредняются, шум гасится, горб остаётся.

Таблицы T1 у нас под рукой нет, поэтому в виджете — своя реализация шума (зерно 185 подобрано так, чтобы метод пропуска, как и у Лагутина, дал $h = 0{,}09$). Жми «Новая выборка» — и увидишь, что на другой реализации горб может проступать хуже, а оптимальное окно — уехать.

Оценка производной

Методы ядерного сглаживания могут быть также применены и для оценивания производных кривой регрессии. Так, первую производную $m'(x)$ можно оценить, используя веса

получаемые дифференцированием правой части равенства (11) по $x$.

По кирпичикам: (11) — дробь $A/B$ с $A = \sum q(z_i)Y_i$, $B = \sum q(z_j)$; производная по $x$ равна $(A'B - AB')/B^2$, и $dz_i/dx = 1/h$. Собираем: $\hat m'(x) = \sum_iw_i^*(x)Y_i \big/ \bigl(h\,B^2\bigr)$. Ядро для этого нужно дифференцируемое — квартическое подходит, у Епанечникова излом на краях.

Пример 3. Зависимость скорости роста от возраста

На рис. 5 у Лагутина представлены ядерные оценки кривых средней скорости $V$ (в см/год) роста $H$ (т. е. $V = H'$) в зависимости от возраста $T$ (в годах) для мальчиков и для девочек.

У девочек наблюдается локальный максимум скорости роста около 12 лет. У мальчиков аналогичный пик ещё более выражен, но происходит примерно на 2 года позже.

Применение метода непараметрической регрессии позволило обнаружить и у мальчиков, и у девочек дополнительный локальный максимум $V$, так называемый средний скачок роста, в возрасте около 8 лет. Другие подходы, основанные на априорной фиксации параметрических моделей, приводят к значительным трудностям в обнаружении этого второго пика.

Смысл по-пацански: заранее загнал рост в красивую формулу с одним пиком — второй пик формула и не покажет, она его просто не умеет. Ядро ничего заранее не знает — поэтому и видит то, что есть.

Метод пропуска (кросс-проверка)

Так же, как и при оценивании плотности, решающую роль играет не вид функции $q(y)$, а правильный выбор ширины «окна сглаживания» $h$. Часто в этом помогает метод пропуска, или кросс-проверка.

Минимизируется по $h$ функция

где $\hat m_i(X_i) = \sum_{j\ne i}w_j(X_i)Y_j \big/ \sum_{j\ne i}w_j(X_i)$ представляет собой оценку для $m(X_i)$ по выборке, в которой $i$-е наблюдение пропущено. Здесь $S(h)$ — сумма квадратов ошибок предсказания $Y_i$ (невязок) на основе выборки без $(X_i, Y_i)$.

По кирпичикам: зачем выкидывать саму точку? Если её оставить, то при $h \to 0$ оценка в $X_i$ равна самому $Y_i$, невязки обнуляются, и «оптимальным» окажется нулевое окно — чистая подгонка под шум. Пропуск заставляет предсказывать точку по соседям, как будто её мы ещё не видели.

Смысл по-пацански: экзамен честный, только если билет не видел заранее. Каждую точку прячем, угадываем по остальным, считаем промахи. Какое окно меньше всех мажет — то и берём.

Для данных из примера 2 у Лагутина был подсчитан ряд значений $S(h)$:

Заметим, что $S(h)$ имеет два локальных минимума: при $h \approx 0{,}1$ и $h \approx 0{,}3$. Уточним оптимальное значение $h$, проводя вычисления по сетке с меньшим шагом:

Таким образом, с точностью до 0,01 функция $S(h)$ имеет глобальный минимум при $h = 0{,}09$.

По кирпичикам: два минимума — не случайность. Малое окно ($\approx 0{,}1$) — это ширина горба, при нём горб «виден». Большое ($\approx 0{,}3$) горб уже размазывает, зато хорошо ловит линейный тренд $1 - x$. Сумма $S(h)$ у обоих близка, поэтому и нужна мелкая сетка. Заметим также, что $S(h)/n \approx 46{,}4/50 \approx 0{,}93$ — близко к $\sigma^2 = 1$, как и должно быть: предсказать сам шум не может никто.

Посчитай сам — вот кросс-проверка Надарая — Ватсона на модельных данных примера 2 (свой генератор шума, поэтому числа не совпадут с таблицей Лагутина один в один; зерно 67 даёт оптимум $h = 0{,}10$ и $S(h^*)/n \approx 1{,}01$ — попробуй другие зёрна, и окно поплывёт):

Окно по ближайшим соседям

Альтернативным методом, применяемым в случае, когда $X_1, \dots, X_n$ располагаются очень неравномерно («где пусто, а где густо»), является модификация ядерной оценки, у которой $h$ зависит от $x$. В качестве такого $h(x)$ часто берут расстояние от точки $x$ до её $k$-го (в смысле близости к $x$) соседа $X_i$. Значение параметра $k$ ($1 \le k \le n$) задаётся исследователем.

По кирпичикам: при постоянном окне в густой области в окно влезает много точек (дисперсия мала, можно было бы и сузить), а в пустой — мало или ни одной (оценка скачет или вообще не определена: $0/0$ в (11)). Окно «до $k$-го соседа» в любой точке содержит ровно $k$ наблюдений — дисперсия везде примерно одинакова.

Смысл по-пацански: в толпе спрашиваешь пятерых ближайших — они в шаге от тебя. В степи пятеро ближайших — за километр. Не фиксируй расстояние, фиксируй число людей, которых опрашиваешь.

Устойчивое сглаживание и LOWESS

В заключение главы обсудим проблему уменьшения влияния выделяющихся наблюдений («выбросов») на ядерные оценки кривой регрессии. К сожалению, локальное усреднение, осуществляемое этими оценками, имеет неограниченную возможность реагировать на «выбросы».

Например, график оценки $\hat m(x)$ в правой половине рис. 4 у Лагутина особенно извилист из-за необходимости компенсировать квадраты остатков в точках с большими значениями ошибок $\varepsilon_i$ (таких, как $\varepsilon_{35} = 1{,}73$ и $\varepsilon_{39} = -2{,}36$).

Смысл по-пацански: взвешенное среднее — это всё ещё среднее. Один горлопан в окне перекрикивает всех соседей, и кривая к нему тянется. Как и в гл. 21 с МНК: квадраты невязок наказывают выбросы так сильно, что оценка ради них прогибается.

Повышение устойчивости сглаживания может быть достигнуто за счёт уменьшения весов наблюдений с большими невязками. Кливленд (Cleveland W. S., 1979) предложил следующий простой алгоритм, реализующий эту идею. Процедура начинается с вычисления пробной оценки, итеративно пересчитываются веса и несколько раз производится повторное сглаживание.

Алгоритм LOWESS (модификация для $m$-мерных предикторных переменных $\boldsymbol X_i = (X_{i1}, \dots, X_{im})$)

(В оригинале LOcally WEighted Scatter plot Smoothing — локально взвешенное сглаживание диаграммы рассеяния.)

1. Задаётся $k$ ($1 \le k \le n$). Для каждого $\boldsymbol X_i$ вычисляется расстояние $h(\boldsymbol X_i)$ до $k$-го соседа в $\mathbb R^m$ и по формуле (9) определяется соответствующий вес $w_i(\boldsymbol X_i)$. Веса всех наблюдений запоминаются в массиве.
2. В окрестности каждой точки $\boldsymbol X_i$ строится локальная линейная аппроксимация поверхности регрессии $m(\boldsymbol X)$ на основе взвешенного метода наименьших квадратов (см. §5 гл. 21) путём минимизации по $\boldsymbol\theta = (\theta_1, \dots, \theta_m)$ функции

3. Используя оценки невязок $\hat\varepsilon_i$, вычисляется выборочная медиана $\hat\sigma = \mathrm{MED}\{|\hat\varepsilon_i|,\ i = 1, \dots, n\}$ (оценка масштаба) и определяются коэффициенты устойчивости $c_i = q\bigl(\hat\varepsilon_i/(6\hat\sigma)\bigr)$, где $q(y) = (15/16)(1 - y^2)^2I_{\{|y|\le1\}}$ — квартическое ядро.
4. Строится локальная линейная аппроксимация, как в пункте 2, но с весами $c_iw_i(\boldsymbol X_i)$.
5. Пункты 3 и 4 повторяются до стабилизации величин $c_i$.

По кирпичикам:

• Шаг 2 вместо константы подгоняет прямую. Локальная прямая лучше локальной константы у краёв данных: у края все соседи лежат с одной стороны, и среднее «съезжает» внутрь, а прямая продолжает наклон.
• Шаг 3 — медиана, а не среднее. Масштаб невязок оценивается устойчиво: выбросы не раздувают $\hat\sigma$.
• Порог $6\hat\sigma$. Медиана $|\hat\varepsilon|$ для нормального шума $\approx 0{,}67\sigma$, так что $6\hat\sigma \approx 4\sigma$. Точка с невязкой больше $\approx 4\sigma$ получает $c_i = 0$ — выкинута совсем; с невязкой около $2\sigma$ — вес примерно $0{,}6$ от нормального (в нашей реализации $c_i$ ещё делится на $q(0) = 15/16$, чтобы «нормальная» точка имела вес 1).
• Шаг 4 — тот же взвешенный МНК, но выбросам веса срезаны.

Смысл по-пацански: сперва сгладил как получится. Потом посмотрел, кто от кривой дальше всех, — и понизил тому голос. Сгладил заново. Ещё раз проверил. Через пару кругов горлопаны молчат, а кривая идёт по нормальным пацанам.

Пример 4. LOWESS на зашумлённой прямой

У Лагутина (по Хардле) показано применение алгоритма LOWESS к моделированным данным. Точки генерировались в соответствии с моделью $Y_i = (1/50)x_i + \varepsilon_i$ ($n = 50$), $x_i = i$, $\varepsilon_i \sim \mathcal N(0, 1)$. Кривой регрессии является прямая с малым коэффициентом наклона, причём дисперсия ошибок $\mathsf D\varepsilon_i$ настолько велика, что систематический прирост $Y$ при увеличении $X$ (так называемый тренд) почти не заметен из-за интенсивного «шума».

Для $k = 25$, $m = 2$ и $\boldsymbol X_i = (1, x_i)$ (т. е. локальная прямая со свободным членом) на рис. 6 у Лагутина представлен также результат сглаживания после двух итераций алгоритма. Построенная оценка уверенно игнорирует «выбросы», расположенные вблизи границ рисунка, и довольно точно воспроизводит теоретический линейный тренд.

По кирпичикам: за всю выборку $Y$ вырастает в среднем лишь на $50/50 = 1$ — ровно одно стандартное отклонение шума. Тренд в одну сигму на всём диапазоне глазом почти не ловится.

В виджете — своя модельная выборка (с четырьмя подсаженными выбросами у краёв): сравни простую локальную прямую и LOWESS после 1–3 итераций.

Сам алгоритм в 40 строк — без numpy, для $m = 2$ (локальная прямая):

Смысл по-пацански: итерация 0 — обычная локальная прямая, у краёв её утаскивают выбросы. После робастификации выбросы получают веса, близкие к нулю, и кривая ложится на тренд.

Проверь себя