§1. Оценивание плотности

Из §1 гл. 9 известно, что несмещённой и состоятельной оценкой неизвестной функции распределения $F(x)$ случайной величины $X$ является эмпирическая функция распределения

которая возрастает скачками величины $\frac1n$ в точках $X_{(i)}$, где $X_{(1)} \le \dots \le X_{(n)}$ — упорядоченные элементы выборки.

А как оценить плотность $p(x) = F'(x)$ (если она существует)? Так как $\hat F_n(x)$ — кусочно-постоянная функция, то её производная равна 0 всюду, за исключением точек скачков $\hat F_n(x)$, и не годится в качестве оценки для $p(x)$.

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

Основная идея состоит в предварительном сглаживании эмпирического распределения за счёт его свёртки (см. П3) с распределениями, имеющими плотности и сходящимися к сосредоточенному в точке 0 распределению. Точнее: рассмотрим случайную величину $Z_n = X + h_nY$, где случайная величина $Y$ имеет (известную) плотность $q(y)$ и не зависит от $X$, а числа $h_n > 0$ и $h_n \to 0$ при $n \to \infty$.

Согласно следствию из формулы преобразования плотности (П8), плотностью случайной величины $h_nY$ служит $\frac1{h_n}q\bigl(\frac y{h_n}\bigr)$. С учётом формулы свёртки (П3) выразим плотность $r_n(z)$ случайной величины $Z_n$:

При $h_n \to 0$ случайные величины $Z_n$ будут сходиться к $X$ по вероятности (П5), а плотности $r_n(z)$ — к плотности $p(z)$ для почти всех $z$ при выполнении, скажем, условия 1 из приведённой ниже теоремы 1. (Достаточно ограниченности плотности $q(y)$ и выполнения при некоторых $c > 0$ и $\varepsilon > 0$ неравенства $q(y) \le c/|y|^{1+\varepsilon}$.)

Заменив $F(x)$ в формуле (2) на эмпирическую функцию распределения $\hat F_n(x)$, получим для $p(z)$ оценку

которую обычно называют оценкой Розенблатта — Парзена (Rosenblatt, 1956; Parzen, 1962).

По кирпичикам: каждое наблюдение $X_i$ заменяем маленьким «холмиком» $\frac1{h_n}q\bigl(\frac{z - X_i}{h_n}\bigr)$ площади 1 с центром в $X_i$ и шириной порядка $h_n$. Складываем холмики и делим на $n$ — получаем плотность площади 1. Гистограмма — это то же самое, только с прямоугольными «кирпичами», прибитыми к фиксированной сетке.

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

Свойства оценки

Какими статистическими свойствами она обладает?

Теорема 1. Пусть выполнены следующие условия:

1. плотность $q(y)$ непрерывна и ограничена, причём $\alpha = \int q^2(y)\,dy < \infty$;
2. $h_n \to 0$ при $n \to \infty$ так, что $nh_n \to \infty$.

Тогда

$$\hat p_n(z) = r_n(z) + \xi_n(z)\big/\sqrt{nh_n}, \tag{4}$$

где $r_n(z) \to p(z)$ при почти всех $z$, а случайные величины $\xi_n(z)$ асимптотически нормальны: $\xi_n(z) \xrightarrow{d} \xi(z) \sim \mathcal N(0, \alpha p(z))$.

(Доказательство этой теоремы можно найти в учебнике [11, с. 58] у Лагутина.)

Естественно, возникает вопрос об оптимальном выборе $q(y)$ и скорости стремления к нулю последовательности $h_n$ при $n \to \infty$. Ответ на него зависит от свойств гладкости оцениваемой плотности $p(x)$. Ограничимся классом плотностей, для которых верно допущение

Д1. Носителем $p(x)$ (то есть множеством $\{x\colon p(x) > 0\}$) является конечный интервал, на котором $p(x)$ дважды непрерывно дифференцируема с условием $\gamma = \int[p''(x)]^2\,dx < \infty$.

Предположим также, что для $q(y)$ справедливы условия

Д2. $\int yq(y)\,dy = 0$ (это всегда верно для чётной функции $q$) и $\beta = \int y^2q(y)\,dy < \infty$.

Совершив замену $y = (z - x)/h_n$ во втором интеграле в равенствах (2), применим формулу Тейлора:

Подставляя результат в представление (4), видим, что

поскольку из формул (3) и (2) имеем $\mathbf M\hat p_n(z) = r_n(z)$, то есть $\mathbf M\xi_n(z) = 0$.

Математическое ожидание во второй из формул представляет собой квадратичный риск оценки $\hat p_n(z)$ (см. §3 гл. 6; абсолютный риск $\mathbf M|\hat p_n(z) - p(z)|$ изучается в монографии [21] у Лагутина). Однако риск зависит от $z$ через неизвестные значения $p(z)$ и $p''(z)$. Чтобы избавиться от этого, усредним риск по всем $z$ с помощью интеграла

главная часть которого (после отбрасывания $o(h_n^4)$) равна

При $h_n \to 0$ первое слагаемое в сумме (5), отвечающее за смещение оценки $\hat p_n(z)$, убывает. Второе слагаемое, отражающее дисперсию («разброс») оценки $\hat p_n(z)$, наоборот, возрастает. Минимум функции $J(h_n)$ достигается при $h_n^* = [\alpha/(\beta^2\gamma)]^{1/5}n^{-1/5}$ (убедитесь!).

По кирпичикам: $J'(h) = \beta^2\gamma h^3 - \alpha/(nh^2) = 0 \iff h^5 = \alpha/(\beta^2\gamma n)$.

Подставляя $h_n^*$ в формулу (5), находим минимальное значение

Таким образом, порядок малости оптимального $h_n^*$ есть $n^{-1/5}$. При этом скорость сходимости $\hat p_n(z)$ к $p(z)$, определяемая величиной $\sqrt{J(h_n^*)}$, составляет лишь $n^{-2/5}$ — в отличие от скорости $n^{-1/2}$, которая имеет место для сходимости $\hat F_n(x)$ к $F(x)$ (см. §2 гл. 12). Это можно объяснить тем, что в оценке значения $p(z)$ при фиксированном $z$ принимают участие лишь $X_i$, находящиеся в некоторой уменьшающейся окрестности точки $z$, а не вся выборка.

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

Оптимальное ядро

Теперь выберем оптимальную плотность $q(y)$. Заметим, что $J(h_n^*)$ в формуле (6) зависит от $q(y)$ через $\alpha^2\beta$. Эта величина инвариантна относительно преобразования $\tilde Y = \sigma Y$, $\tilde q(y) = \frac1\sigma q\bigl(\frac y\sigma\bigr)$, где $\sigma > 0$. Поэтому, не ограничивая общности, можно считать, что $\beta = \int y^2q(y)\,dy = 1$. Тогда (ввиду Д2) приходим к задаче минимизации $\alpha = \int q^2(y)\,dy$ при условиях $\int q(y)\,dy = \int y^2q(y)\,dy = 1$, $\int yq(y)\,dy = 0$. Её решением является функция

ранее встречавшаяся в теореме 4 гл. 8. При выборе $\beta = 1/5$ оптимальная плотность $q^*$ имеет более простой вид:

В статистической литературе её обычно именуют ядром Епанечникова (Epanechnikov, 1969), хотя оно было ещё раньше введено Бартлеттом (Bartlett, 1963). (Ядро — это непрерывная ограниченная чётная функция с единичным интегралом $\int q(y)\,dy = 1$, не обязательно неотрицательная.)

Ядра — чистый изумруд…

— А. С. Пушкин, «Сказка о царе Салтане»

Насколько важно использовать именно такое ядро? Конкуренцию ему могут составить ряд ядер, перечисленных в следующей таблице:

• Квартическое ядро, в отличие от ядра Епанечникова, дифференцируемо в точках $-1$ и $1$.
• В случае треугольного ядра можно быстро производить пересчёт оценки плотности $\hat p_n(z)$ при увеличении переменной $z$ с заданным шагом $\Delta z$. Это важно при интерактивном построении графиков $\hat p_n(z)$ для «окон сглаживания» $h_n$ разной ширины (см. пример 1 ниже).
• Ядро Гаусса (или нормальное) бесконечно дифференцируемо на всей оси. Однако оценка $\hat p_n(z)$ вычисляется медленно по причине многократного подсчёта значений экспонент.
• Достоинством прямоугольного ядра является его простой вид. Формально эту функцию нельзя называть ядром из-за разрывов в точках $-1$ и $1$, но её можно рассматривать как предел последовательности ядер, имеющих форму трапеции.

В последнем столбце таблицы указаны эффективности

ядер по отношению к оптимальному ядру Епанечникова $q^*$ (см. формулу (6)). На основе этой таблицы приходим к заключению, что выбор ядра мало влияет на величину среднеквадратичной ошибки оценки $\hat p_n(z)$.

Смысл по-пацански: спорить о форме ядра — всё равно что спорить, лопатой или совком насыпать кучку песка. Разница — проценты. А вот сколько песка в кучке (ширина окна) — решает всё.

Ядра высших порядков

Отметим, что если неизвестная плотность $p(x)$ имеет непрерывные производные более высокого чётного порядка $m > 2$, то можно добиться и большей, чем $n^{-2/5}$, скорости сходимости за счёт использования ядер, принимающих отрицательные значения. При выполнении условий $\int q(y)\,dy = 1$, $\int y^kq(y)\,dy = 0$ при $1 \le k < m$, $\int y^mq(y)\,dy = c \ne 0$ прежние рассуждения приводят к задаче минимизации функционала

решением которой на отрезке $[-1, 1]$ является некоторый полином степени $m$. В частности, для $m = 4$ (при $c = -1/21$) получается ядро вида $(15/32)(7y^4 - 10y^2 + 3)I_{\{|y|\le1\}}$. Скорость сходимости соответствующей ядерной оценки равна $n^{-4/9}$.

В общем случае погрешность оптимальной оценки имеет порядок малости $n^{-m/(2m+1)}$, который приближается к $n^{-1/2}$ при увеличении $m$. Это объясняется тем, что для более гладких плотностей к оцениванию значения $p(x)$ привлекаются элементы выборки, лежащие во всё более широких окрестностях точки $x$.

Недостатком таких оценок является возможность появления у них отрицательных значений, что противоречит тому, что $p(x) \ge 0$.

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

Пример 1. Выбор окна

При практическом оценивании плотности $p(x)$ по заданной реализации выборки важен не столько выбор вида ядра, сколько правильное определение ширины «окна сглаживания» $h_n$.

С помощью таблицы случайных чисел T1 у Лагутина была смоделирована выборка размера $n = 100$ из распределения с плотностью $p(z) = 2(I_{[0{,}1;\,0{,}4]} + I_{[0{,}6;\,0{,}8]})$. На рис. 3 приведены графики оценок $\hat p_n(z)$, вычисленных на основе треугольного ядра для следующих значений $h_n$: 0,1; 0,02; 0,5.

Выбор слишком малого $h_n$ приводит к быстро меняющейся, неустойчивой оценке, так как $\hat p_n(z)$ опирается только на небольшое количество наблюдений из узкой окрестности точки $z$ (велико второе слагаемое в формуле (5), отвечающее за дисперсию оценки).

С другой стороны, слишком большое значение $h_n$ влечёт чрезмерное сглаживание оцениваемой плотности, что не позволяет выявить её характерные особенности (такие, например, как наличие нескольких максимумов). Оценка $\hat p_n(z)$ в этом случае оказывается сильно смещённой.

В виджете — своя модельная выборка из той же плотности (таблицы T1 под рукой нет); жми кнопки $h = 0{,}02$, $0{,}1$, $0{,}5$:

Смысл по-пацански: при $h = 0{,}02$ — «расчёска» из пиков на каждой точке, при $h = 0{,}5$ — один пологий бугор, и два плато слились. Истина — посередине.

Проверь себя