Метод максимального правдоподобия

Заметка об оценках судила о них, когда они уже были заданы. Метод максимального правдоподобия (maximum likelihood) их строит: среди всех значений параметра он выбирает то, при котором наблюдённые данные наиболее вероятны. Рецепт работает для любой модели, которая приписывает данным вероятность или плотность, от монеты до языковой модели, и большинство распространённых функций потерь машинного обучения — это он в замаскированном виде. На больших выборках он к тому же почти так же хорош, как вообще может быть оценка. В этой заметке мы определяем правдоподобие, вычисляем оценки для распространённых семейств, выводим из него квадратичную ошибку и кросс-энтропию и объясняем точность оценки через информацию Фишера (Fisher information).

Обозначения. $x_1, \dots, x_n$ — наблюдённые значения; считается, что они н.о.р. с распределением $p_\theta$ (функцией вероятности или плотностью). $L(\theta) = \prod_i p_\theta(x_i)$ — правдоподобие (likelihood), $\ell(\theta) = \log L(\theta)$ — логарифм правдоподобия (log-likelihood), $\hat\theta = \arg\max_\theta L(\theta)$ — оценка максимального правдоподобия (maximum likelihood estimate, MLE).

Правдоподобие#

Правдоподобие — это вероятность наблюдённых данных, рассматриваемая как функция параметра. После 7 орлов в 10 независимых бросках вероятность именно этой последовательности при вероятности орла $p$ равна

$$L(p) = p^7 (1 - p)^3 .$$

Данные фиксированы; меняется $p$. $L$ — не распределение вероятностей на значениях $p$: её интеграл не обязан быть равен 1, и смысл имеют только отношения $L(p_1)/L(p_2)$, которые показывают, насколько одно значение объясняет данные лучше другого. Чтобы превратить её в распределение на значениях $p$, нужно априорное распределение — предмет заметки о байесовском выводе.

Рисунок

Правдоподобие вероятности орла

Максимизация правдоподобия#

Произведения многих малых чисел неудобны, а логарифм превращает их в суммы, не сдвигая максимума. Для монеты $\ell(p) = 7 \log p + 3 \log(1 - p)$, и, приравняв производную к нулю, получаем

$$\ell'(p) = \frac{7}{p} - \frac{3}{1-p} = 0 \quad \Longrightarrow \quad \hat p = \frac{7}{10} .$$

Производная логарифма правдоподобия называется функцией вклада (score). То же вычисление для распространённых семейств даёт оценки, которые задним числом кажутся очевидными:

Модель	MLE	Примечание
$\mathrm{Bernoulli}(p)$	$\hat p = k/n$, доля успехов
$N(\mu, \sigma^2)$	$\hat\mu = \bar x$, $\hat\sigma^2 = \frac1n \sum (x_i - \bar x)^2$	$\hat\sigma^2$ смещена с множителем $(n-1)/n$
$\mathrm{Poisson}(\lambda)$	$\hat\lambda = \bar x$
$\mathrm{Exp}(\lambda)$	$\hat\lambda = 1/\bar x$	смещена вверх: $E[\hat\lambda] = \frac{n}{n-1}\lambda$ при $n > 1$, бесконечно при $n = 1$
$U(0, \theta)$	$\hat\theta = \max_i x_i$	максимум достигается на границе, где $\ell' \ne 0$

Последняя строка — оценка правого конца равномерного распределения из предыдущей заметки: метод максимального правдоподобия не обещает несмещённости. Зато он обещает инвариантность (invariance): оценка максимального правдоподобия любой функции $g(\theta)$ равна $g(\hat\theta)$, поэтому для стандартного отклонения она равна $\sqrt{\hat\sigma^2}$, а для шансов — $\hat p / (1 - \hat p)$.

Функции потерь — это правдоподобия#

Регрессионная модель $f_w(x)$ с гауссовским шумом утверждает, что $y = f_w(x) + \varepsilon$, $\varepsilon \sim N(0, \sigma^2)$. Отрицательный логарифм правдоподобия обучающей выборки равен

$$-\ell(w) = \sum_{i=1}^n \frac{(y_i - f_w(x_i))^2}{2\sigma^2} + n \log\big(\sigma\sqrt{2\pi}\big) ,$$

и минимизировать его по $w$ — значит минимизировать сумму квадратов ошибок. Классификатор, который выдаёт вероятности классов $q_w(k \mid x)$, приписывает наблюдённой метке $y_i$ вероятность $q_w(y_i \mid x_i)$, поэтому

$$-\ell(w) = -\sum_{i=1}^n \log q_w(y_i \mid x_i) ,$$

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

Рисунок

Функции потерь — это минус логарифм правдоподобия

Этот взгляд объясняет решения, которые иначе выглядят произвольными. Абсолютная ошибка — отрицательный логарифм правдоподобия шума Лапласа — робастна, потому что распределение Лапласа допускает редкие большие остатки; пуассоновская регрессия для счётных данных использует правдоподобие Пуассона; дисперсия шума $\sigma^2$, не влияющая на подобранные $w$, становится важной, как только модель должна сообщать о собственной неопределённости. Наконец, на большой выборке средний отрицательный логарифм правдоподобия приближается к $E_{p}[-\log p_\theta(X)]$ — кросс-энтропии между распределением данных и моделью, поэтому максимизация правдоподобия минимизирует прямую KL-дивергенцию (forward KL) распределения данных относительно модели; этот взгляд принят в курсе генеративных моделей.

Информация Фишера#

Насколько точна $\hat\theta$? Ответ подсказывает рисунок: острый пик логарифма правдоподобия точно задаёт параметр, пологий — нет. Ожидаемая кривизна на одно наблюдение — это информация Фишера

$$I(\theta) = -E_\theta\big[\tfrac{\partial^2}{\partial\theta^2} \log p_\theta(X)\big] = \operatorname{Var}_\theta\big(\tfrac{\partial}{\partial\theta} \log p_\theta(X)\big) ,$$

то есть дисперсия функции вклада. Для одного бернуллиевского наблюдения $I(p) = 1/(p(1-p))$. Кривизна, которую показывает рисунок, — другая величина, зависящая от данных: наблюдаемая информация (observed information) $J_n(\hat\theta) = -\ell_n''(\hat\theta)$ всей выборки; в истинном $\theta$ её математическое ожидание равно $n I(\theta)$, и при условиях ниже при больших $n$ эти величины близки. Для монеты они совпадают в точке оценки: $J_n(\hat p) = n/(\hat p(1 - \hat p)) = n I(\hat p)$. При условиях регулярности (в частности, модель задана верно и идентифицируема, т. е. разные $\theta$ дают разные распределения, носитель не зависит от $\theta$, логарифм правдоподобия гладкий, а истинное $\theta$ не лежит на границе) два результата делают оценку максимального правдоподобия оценкой по умолчанию:

• Асимптотическая нормальность. $\hat\theta \approx N\big(\theta, \frac{1}{n I(\theta)}\big)$ при большом $n$. Для монеты стандартная ошибка равна $\sqrt{p(1-p)/n}$: $0.145$ после 10 бросков и $0.046$ после 100, а симуляция $\hat p$ по 100 броскам даёт разброс $0.0458$.
• Эффективность. Неравенство Крамера–Рао (Cramér–Rao bound) утверждает, что ни у одной несмещённой оценки дисперсия не бывает меньше $1/(n I(\theta))$. Оценка максимального правдоподобия достигает этой границы асимптотически, поэтому на больших выборках ни одна разумная оценка не точнее.

При $d$ параметрах $I(\theta)$ — матрица $d \times d$, а обратная к ней матрица, делённая на объём выборки, — асимптотическая ковариационная матрица оценки. Роль информации Фишера в машинном обучении не ограничивается интервалами погрешности: методы естественного градиента (natural gradient) используют её для предобуславливания, а методы непрерывного обучения (continual learning), например elastic weight consolidation, — чтобы измерить, какие веса важны.

Когда метод максимального правдоподобия не работает#

• Малые выборки. После 3 орлов в 3 бросках $\hat p = 1$: решка объявляется невозможной. Слово, ни разу не встретившееся при обучении, получает в n-граммной модели нулевую вероятность, поэтому такие модели сглаживают, как в заметке об n-граммных языковых моделях. Априорные распределения решают эту проблему систематически.
• Неограниченное правдоподобие. В смеси двух гауссиан компонента с центром в одной точке данных и дисперсией, стремящейся к нулю, неограниченно увеличивает правдоподобие: для пяти точек логарифм правдоподобия равен $-7.64$, $-3.06$ и $3.84$ при ширине компоненты $10^{-1}$, $10^{-3}$, $10^{-6}$. Максимума не существует, и процедуры подгонки должны ограничивать дисперсии снизу.
• Слишком много параметров. Когда число параметров растёт вместе с объёмом данных, оценка максимального правдоподобия может быть несостоятельной. Если пар измерений много, у каждой пары своё среднее, а дисперсия общая, то оценка максимального правдоподобия дисперсии сходится к половине её истинного значения (задача Неймана–Скотта, Neyman–Scott problem); перепараметризованные модели без регуляризации точно так же подгоняются под шум.
• Неверная модель. Если истинное распределение не совпадает ни с одним $p_\theta$, оценка максимального правдоподобия сходится к члену семейства, ближайшему по KL-дивергенции, а он может плохо описывать то, что важно, например хвосты.

Что изменилось#

	До	С методом максимального правдоподобия
Поиск оценки	новая идея для каждой задачи	один рецепт: максимизировать $\ell(\theta) = \sum \log p_\theta(x_i)$
Функции потерь	выбираются по привычке	отрицательные логарифмы правдоподобия модели шума: MSE, MAE, кросс-энтропия
Насколько точна оценка	выясняется в каждом случае отдельно	$\hat\theta \approx N(\theta, 1/(nI(\theta)))$; неравенство Крамера–Рао
Когда не работает		малые выборки, неограниченное правдоподобие, слишком много параметров, неверная модель

Ограничения#

• Модель нужно задать. Метод максимального правдоподобия находит лучший параметр внутри семейства; выбор семейства и его проверка лежат вне метода.
• Асимптотическим гарантиям нужны регулярность и большое $n$. При малом числе наблюдений, оценках на границе или параметрах, которые меняют носитель (правый конец равномерного распределения), нормальное приближение и неравенство Крамера–Рао неприменимы.
• Точечная оценка скрывает неопределённость. Сама по себе оценка максимального правдоподобия не говорит, насколько ей можно доверять; для этого нужны информация Фишера, кривая правдоподобия, доверительный интервал или апостериорное распределение.

Что нужно запомнить#

• Правдоподобие $L(\theta) = \prod_i p_\theta(x_i)$ — вероятность данных как функция параметра; смысл имеют только отношения.
• Оценка максимального правдоподобия максимизирует $\ell(\theta) = \sum_i \log p_\theta(x_i)$; для распространённых семейств это доля, среднее, дисперсия с $1/n$, $1/\bar x$, максимум.
• Квадратичная ошибка, абсолютная ошибка и кросс-энтропия — отрицательные логарифмы правдоподобия гауссовской, лапласовской и категориальной моделей; метод максимального правдоподобия минимизирует KL-дивергенцию распределения данных относительно модели.
• Информация Фишера $I(\theta) = \operatorname{Var}(\text{вклад})$; $\hat\theta \approx N(\theta, 1/(nI(\theta)))$; ни одна несмещённая оценка не имеет дисперсии меньше $1/(nI(\theta))$.
• Малые выборки, неограниченное правдоподобие, растущее число параметров и неверно специфицированные модели ломают метод.

Код#

Правдоподобие для монеты, метод наименьших квадратов как метод максимального правдоподобия, смещение на малых выборках и неограниченное правдоподобие:

import numpy as np, math
from scipy import stats, optimize

rng = np.random.default_rng(0)

# 7 heads in 10 and 70 in 100: the MLE, the likelihood ratio of p = 0.5, the standard error from Fisher information
for k, n in ((7, 10), (70, 100)):
    loglik = lambda p: k * math.log(p) + (n - k) * math.log(1 - p)
    p_hat = optimize.minimize_scalar(lambda p: -loglik(p), bounds=(1e-6, 1 - 1e-6), method='bounded').x
    print(n, round(p_hat, 3), f'{math.exp(loglik(0.5) - loglik(0.7)):.2g}', round(math.sqrt(0.7 * 0.3 / n), 3))
# 10 0.7 0.44 0.145 | 100 0.7 0.00027 0.046
print(round(float((rng.binomial(100, 0.7, 200_000) / 100).std()), 4))   # 0.0458: the simulated spread of p-hat

# least squares is maximum likelihood under Gaussian noise: both fits of a line give the same coefficients
x = rng.uniform(0, 1, 50)
y = 2 + 3 * x + 0.5 * rng.normal(size=50)
ls = np.polyfit(x, y, 1)
nll = lambda w: -stats.norm.logpdf(y, w[1] + w[0] * x, 0.5).sum()
ml = optimize.minimize(nll, [0.0, 0.0]).x
print(ls.round(4), ml.round(4))                                   # [2.7233 2.1203] [2.7233 2.1203]

# small-sample bias of MLEs, n = 5: the normal variance (true 1) and the exponential rate (true 1)
print(round(float(rng.normal(size=(200_000, 5)).var(axis=1).mean()), 3), round(float((1 / rng.exponential(1, (200_000, 5)).mean(axis=1)).mean()), 3))   # 0.798 1.249

# an unbounded likelihood: a two-component Gaussian mixture with one component shrinking onto a data point
data = np.array([-1.2, 0.3, 0.5, 1.1, 2.4])
for s in (1e-1, 1e-3, 1e-6):
    print(s, round(float(np.log(0.5 * stats.norm.pdf(data, data[0], s) + 0.5 * stats.norm.pdf(data, data.mean(), data.std())).sum()), 2))
# 0.1 -7.64 | 0.001 -3.06 | 1e-06 3.84

Дальше#

Метод максимального правдоподобия возвращает одно значение и считает параметр фиксированным, но неизвестным. После 3 орлов в 3 бросках он объявляет решку невозможной, потому что не может выразить, что монета, которая выглядит честной, априори, скорее всего, примерно честна. Следующая заметка считает сам параметр случайным, с априорным распределением, которое формула Байеса обновляет до апостериорного, и показывает, что распространённые регуляризаторы — это замаскированные априорные распределения.

Литература#

• R. A. Fisher. On the Mathematical Foundations of Theoretical Statistics. Philosophical Transactions of the Royal Society A, 222, 1922.
• L. Wasserman. All of Statistics. Springer, 2004. Глава 9 (метод максимального правдоподобия, информация Фишера, асимптотическая нормальность, эффективность).
• G. Casella, R. L. Berger. Statistical Inference, 2nd ed. Duxbury, 2002. Раздел 7.2 (метод максимального правдоподобия), раздел 7.3 (неравенство Крамера–Рао), глава 10 (асимптотика).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Разделы 1.2.4–1.2.5 (метод максимального правдоподобия для гауссовского распределения, подгонка кривой как метод максимального правдоподобия), раздел 9.2.1 (сингулярности правдоподобия смеси).
• J. Neyman, E. L. Scott. Consistent Estimates Based on Partially Consistent Observations. Econometrica, 16(1), 1948.