§4. Метод максимального правдоподобия: какое θ лучше всего объясняет данные

Метод получил распространение после статьи Р. Фишера 1912 года, где было доказано, что получаемые им оценки — асимптотически наиболее точные при выполнении условий регулярности (теорема 1 ниже).

Функция правдоподобия

Для простоты пусть элементы выборки $X_i$ имеют дискретное распределение $f(x, \boldsymbol\theta) = \mathbf{P}(X_1 = x)$, где $\boldsymbol\theta = (\theta_1, \dots, \theta_m)$ — вектор неизвестных параметров. Тогда совместная вероятность выборки

зависит от $n + m$ аргументов. Рассматриваемая как функция от $\theta_1, \dots, \theta_m$ при фиксированных значениях $x_1, \dots, x_n$, она называется функцией правдоподобия (по-английски likelihood) и обычно обозначается $L(\boldsymbol\theta)$. Её можно считать мерой правдоподобия значения $\boldsymbol\theta$ при заданной реализации $\mathbf x$.

Представляется разумным взять в качестве оценок наиболее правдоподобные значения $\tilde\theta_1, \dots, \tilde\theta_m$ — те, что максимизируют $L(\boldsymbol\theta)$. Такие оценки называют оценками максимального правдоподобия (ОМП).

Часто проще искать максимум $\ln L(\boldsymbol\theta)$ — он совпадает с максимумом $L$, ведь логарифм монотонен. А произведение логарифм превращает в сумму.

Смысл по-пацански: данные — это улики. Перебираешь подозреваемых ($\theta$) и спрашиваешь: при каком из них такие улики выглядели бы наиболее естественно? Того и берёшь.

Примеры

Пример 5. Схема Бернулли. $f(x, \theta) = \theta^x(1 - \theta)^{1 - x}$, $x \in \{0, 1\}$. Функция правдоподобия

— многочлен степени $n$ с горбом внутри $(0, 1)$. Максимизируем логарифм $\ln L = s_n\ln\theta + (n - s_n)\ln(1 - \theta)$. Производная:

ОМП в схеме Бернулли — частота успехов.

Как и в §3, для непрерывных моделей $f(x, \theta)$ — плотность $X_1$, а $L(\theta) = \prod f(x_i, \theta)$.

Пример 6. Сдвиг показательного закона с плотностью $f(x, \theta) = e^{-(x - \theta)}I_{\{x \ge \theta\}}$:

Функция растёт по $\theta$ до точки $x_{(1)}$, а дальше обрывается в ноль — ведь при $\theta > x_{(1)}$ наблюдение $x_{(1)}$ стало бы невозможным. ОМП $\tilde\theta = x_{(1)}$. Она отлична от оценки метода моментов $\hat\theta = \bar x - 1$ из §2 и намного точнее. Заметим, что здесь $L(\theta)$ не гладкая, и ОМП нельзя найти, приравнивая нулю производную.

Уравнение правдоподобия

Если $L(\boldsymbol\theta)$ гладко зависит от $\theta_1, \dots, \theta_m$, ОМП — компоненты решения (вообще говоря, нелинейной) системы

Иногда решение можно найти явно (задача 2), но чаще приходится вычислять его приближённо итерационными методами — например, методом Ньютона из §5.

Пример 7. Закон Вейбулла — Гнеденко с функцией распределения $F_{\alpha, \beta}(x) = 1 - \exp\{-x^\alpha/\beta\}$ при $x > 0$ используется в моделях прочности материалов (глава 4). Система (2) после упрощений приводится к виду

Доказано, что первое уравнение (а значит, и вся система) при любых $x_1, \dots, x_n$, не равных одновременно 1, имеет единственное решение. Решается численно.

Когда правдоподобие неограниченно

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

Пример 8. Смесь нормальных. Рассмотрим смесь плотностей (глава 5)

где $\varepsilon \in (0, 1)$, $\sigma_1, \sigma_2 > 0$, а $p$ — любая положительная плотность, скажем, стандартная нормальная. Плотность выборки — произведение таких сумм, и после раскрытия скобок одно из слагаемых равно

Положим $\mu_1 = x_1$ и устремим $\sigma_1 \to 0$: множитель $\frac{1}{\sigma_1}p(0)$ улетает в бесконечность при любых фиксированных остальных параметрах. Значит, $L$ неограниченна, и глобальной ОМП нет. Хотя условия теоремы 1 ниже для смеси нормальных выполняются, и при достаточно большом $n$ правдоподобие с вероятностью 1 имеет локальный максимум — решение системы (2).

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

Пример 9. Пусть в модели сдвига-масштаба $\frac1\sigma p\!\left(\frac{x - \mu}{\sigma}\right)$ плотность такова, что $x^{1+\varepsilon}p(x) \to \infty$ при $x \to \infty$ для любого $\varepsilon > 0$ — хвосты тяжелее любой степени, например $p(x) = \frac12(1 + |x|)^{-1}[1 + \ln(1 + |x|)]^{-2}$. Тогда при $n > 1$ глобальная ОМП не существует. Идея: поставим $\mu = x_{(1)} - c\sigma$ с $p(c) > 0$ и устремим $\sigma \to 0$. Одна точка «сидит» у центра, а остальные улетают в хвост, который убывает так медленно, что не компенсирует рост множителя $\sigma^{-n}$. (Обратное тоже доказано: если $p$ ограничена, непрерывна и $|x|^{1+\varepsilon}p(x)$ ограничена при некотором $\varepsilon > 0$, то с вероятностью 1 глобальная ОМП существует при большом $n$ и сходится к истине.)

Свойства ОМП

Координаты точки максимума зависят от реализации: $\tilde\theta_j = \tilde\theta_j(x_1, \dots, x_n)$. Подставив случайные $X_i$, получим случайные величины. Каковы их свойства? Ограничимся скалярным параметром. Пусть помимо R1–R4 из §3 выполнены дополнительные условия регулярности:

• R5) распределения $F_\theta$ различны при разных $\theta \in \Theta$;
• R6) плотность $f(x, \theta)$ при каждом $x \in A$ трижды непрерывно дифференцируема по $\theta$;
• R7) $\int f(x, \theta)\,dx$ можно дважды дифференцировать по $\theta$ под знаком интеграла;
• R8) существует функция $h(x)$ такая, что при всех $x \in A$ и $\theta \in \Theta$ выполнено $\left|\frac{\partial^3}{\partial\theta^3}\ln f(x, \theta)\right| \le h(x)$, причём $\mathbf{M}h(X_1) < \infty$.

Теорема 1. При выполнении условий R1–R8 для достаточно большого $n$ с вероятностью 1 существует решение $\tilde\theta = \tilde\theta_n$ уравнения правдоподобия

$$\sum_{i=1}^n \frac{\partial}{\partial\theta}\ln f(X_i, \theta) = 0, \tag{3}$$

дающее сильно состоятельную оценку: $\tilde\theta_n \to \theta$ почти наверное, причём

$$\sqrt n\,(\tilde\theta_n - \theta) \xrightarrow{d} \xi \sim N\bigl(0,\ 1/I_1(\theta)\bigr),$$

где $I_1(\theta)$ — информация Фишера одного наблюдения (§3).

Теорема 2. Пусть выполнены условия R1–R8 и $\hat\theta$ — некоторая асимптотически нормальная оценка: $\sqrt n\,(\hat\theta_n - \theta) \xrightarrow{d} \zeta \sim N(0, \sigma^2(\theta))$, где асимптотическая дисперсия $\sigma^2(\theta)$ непрерывно зависит от $\theta$. Тогда

$$\sigma^2(\theta) \ge 1/I_1(\theta) \quad \text{при всех } \theta \in \Theta. \tag{4}$$

Из неравенства (4) заключаем, что ОМП (при выполнении R1–R8) имеет наименьшую возможную асимптотическую дисперсию $1/I_1(\theta)$. Это свойство называется асимптотической эффективностью.

Тонкость: Ле Кам в 1953 году доказал, что для разрывных $\sigma^2(\theta)$ неравенство (4) может нарушаться, но не более чем на множестве меры нуль. Контрпример — в задаче 7.

Смысл по-пацански: при большой выборке ОМП — чемпион среди всех «приличных» оценок. Её разброс ровно на потолке, который разрешает информация Фишера.

Пощупай руками

Пощёлкай модели:

• Бернулли — горб, максимум ровно в частоте; подвигай $n$ и смотри, как горб сужается;
• нормальный — парабола в логарифме, максимум в $\bar X$;
• сдвиг показательного — функция растёт и обрывается на $X_{(1)}$: максимум на изломе;
• Коши при $n = 5$ — часто несколько горбов, уравнение (3) имеет несколько корней (вопрос 3). Подвигай $n$ — лишние горбы исчезают.

Проверим теорему 1 на Коши: ОМП против медианы. По теореме 1 асимптотическая дисперсия ОМП $1/I_1 = 2$, у медианы $\pi^2/4 \approx 2{,}47$. Значит, ОМП примерно на 23 % эффективнее ($e = \pi^2/8 \approx 1{,}234$).

Проверь себя