Свойства оценок, профильное правдоподобие и примеры анализа данных

Статистический анализ данных

Дмитрий В. Наумов (ОИЯИ)

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

Смещение, дисперсия и средний квадрат ошибки

Смещение: \[ \operatorname{смещение}(\widehat\theta) = \mathbb E[\widehat\theta]-\theta. \]

Средний квадрат ошибки: \[ \mathbb E[(\widehat\theta-\theta)^2]. \]

Он раскладывается как \[ \boxed{ \mathbb E[(\widehat\theta-\theta)^2] = \operatorname{Var}(\widehat\theta) + \operatorname{смещение}^2(\widehat\theta). } \]

Несмещённость не гарантирует оптимальность

Оценка может быть несмещённой, но иметь большую дисперсию: \[ \mathbb E[\widehat\theta]=\theta, \qquad \operatorname{Var}(\widehat\theta)\ \text{велика}. \]

Небольшое смещение иногда уменьшает общий средний квадрат ошибки: \[ \mathbb E[(\widehat\theta-\theta)^2] = \operatorname{Var} + \operatorname{смещение}^2. \]

Свойства процедуры необходимо оценивать вместе, а не выбирать оценку только по одному критерию.

Интерактив: делитель \(N\) или \(N-1\)?

Генерируем гауссовы выборки с истинной дисперсией \(\sigma^2=1\) и сравниваем две оценки: \[ \widehat{\sigma^2}_{\mathrm{ММП}} = \frac1N\sum_i(X_i-\langle X\rangle)^2, \qquad s^2 = \frac1{N-1}\sum_i(X_i-\langle X\rangle)^2. \]

Состоятельность

Оценка называется состоятельной, если при росте выборки она сходится к истинному параметру: \[ \boxed{ \widehat\theta_N \xrightarrow[N\to\infty]{P} \theta. } \]

  • Смещение должно исчезать или становиться пренебрежимо малым.
  • Дисперсия должна стремиться к нулю.
  • Состоятельность — асимптотическое свойство.
  • Для конечной выборки оценка всё ещё может быть смещённой и иметь большую ошибку.

Зачем нужна нижняя граница точности

Для \(X_i\sim\mathrm N(\mu,\sigma^2)\) при известной \(\sigma\) можно придумать несколько оценок \(\mu\): \[ \widehat\mu_1=\langle X\rangle, \qquad \widehat\mu_2=\operatorname{median}(X_1,\ldots,X_N), \qquad \widehat\mu_3=\frac{X_1+X_N}{2}. \]

Все выглядят разумно, но используют информацию по-разному: \[ \operatorname{Var}(\widehat\mu_1)=\frac{\sigma^2}{N}, \qquad \operatorname{Var}(\widehat\mu_2)\approx\frac{\pi}{2}\frac{\sigma^2}{N}, \qquad \operatorname{Var}(\widehat\mu_3)=\frac{\sigma^2}{2}. \]

  • Среднее использует все наблюдения наиболее эффективно для гауссовой модели.
  • Медиана устойчивее к выбросам, но в чистой гауссовой модели платит дисперсией.
  • Оценка по двум точкам почти выбрасывает остальные данные.

Граница Крамера — Рао

Скор-функция показывает, как логарифм правдоподобия реагирует на параметр: \[ U(\theta) = \frac{\partial\ell(\theta)}{\partial\theta}. \]

Информация Фишера: \[ \boxed{ I(\theta) = \mathbb E\!\left[U^2(\theta)\right] = -\mathbb E\!\left[ \frac{\partial^2\ell}{\partial\theta^2} \right]. } \]

При выполнении условий регулярности для несмещённой оценки: \[ \boxed{ \operatorname{Var}(\widehat\theta) \ge \frac{1}{I(\theta)}. } \]

Это нижняя граница: лучше неё несмещённая оценка в данной модели быть не может. Она не утверждает, что ММП всегда несмещён или всегда достигает границы; эти свойства нужно проверять отдельно.

Раскрываем вторую производную

  • Пусть для одного наблюдения \(s_\theta(x)=\partial_\theta\ln f(x\mid\theta)\).

  • Из нормировки плотности: \[ \int f(x\mid\theta)\,dx=1. \]

  • Дифференцируем один раз: \[ \mathbb E[s_\theta(X)] = \int \partial_\theta f(x\mid\theta)\,dx =0. \]

  • Дифференцируем второй раз: \[ \partial_\theta f = f\,\partial_\theta\ln f = f\,s_\theta, \] \[ \partial_\theta^2 f = \partial_\theta(f s_\theta) = (\partial_\theta f)s_\theta +f\,\partial_\theta s_\theta = f\left[ s_\theta^2+\partial_\theta^2\ln f \right]. \]

Интегрируем результат

  • Интегрируем: \[ 0 = \int \partial_\theta^2 f\,dx = \int f\left[ s_\theta^2+\partial_\theta^2\ln f \right]dx = \mathbb E\!\left[ s_\theta^2(X)+\partial_\theta^2\ln f(X\mid\theta) \right]. \]

  • Поэтому \[ \boxed{ \mathbb E[s_\theta^2(X)] = -\mathbb E[\partial_\theta^2\ln f(X\mid\theta)]. } \]

Пример: гауссово среднее

Для одного наблюдения: \[ \ln f(x\mid\mu) = -\frac{(x-\mu)^2}{2\sigma^2} +\mathrm{const}, \qquad s_\mu(x) = \frac{x-\mu}{\sigma^2}. \]

Тогда \[ I_1(\mu) = \mathbb E[s_\mu^2(X)] = \frac{1}{\sigma^4}\mathbb E[(X-\mu)^2] = \frac{1}{\sigma^2}. \]

Для \(N\) независимых наблюдений: \[ I_N(\mu)=\frac{N}{\sigma^2}, \qquad \operatorname{Var}(\widehat\mu)\ge\frac{\sigma^2}{N}. \]

Среднее \(\langle X\rangle\) имеет ровно такую дисперсию, значит в этой модели достигает границы Крамера — Рао.

Как читать информацию Фишера

  • Информация Фишера измеряет чувствительность распределения к параметру.
  • Чем острее максимум логарифма правдоподобия, тем точнее определяется параметр.
  • Для независимых наблюдений информация складывается: \(I_N=N I_1\).
  • Поэтому типичная ошибка регулярной эффективной оценки уменьшается как \(1/\sqrt N\).
  • Если оценка достигает границы Крамера — Рао, её называют эффективной.

Форма правдоподобия

Кривизна около максимума

Разложим логарифм правдоподобия около максимума: \[ \ell(\theta) \approx \ell(\widehat\theta) - \frac12 (\theta-\widehat\theta)^2 \left[ -\ell''(\widehat\theta) \right]. \]

Отсюда локальная оценка дисперсии: \[ \boxed{ \sigma_{\widehat\theta}^2 \approx \left[ -\frac{\partial^2\ell}{\partial\theta^2} \bigg|_{\widehat\theta} \right]^{-1}. } \]

Это гауссово приближение: оно может нарушаться для малых выборок, границ и асимметричных функций правдоподобия.

Относительное правдоподобие

Абсолютная нормировка \(L\) не важна. Удобно сравнивать с максимумом: \[ \lambda(\theta) = \frac{L(\theta)}{L(\widehat\theta)}. \]

Часто используют \[ \boxed{ q(\theta) = -2\ln\lambda(\theta) = -2\left[\ell(\theta)-\ell(\widehat\theta)\right]. } \]

  • В максимуме \(q(\widehat\theta)=0\).
  • Чем хуже параметр описывает данные, тем больше \(q(\theta)\).
  • Для локально гауссова правдоподобия: \[ q(\theta)\approx\frac{(\theta-\widehat\theta)^2}{\sigma_{\widehat\theta}^2}. \]

Квадратичный профиль среднего гаусса

Для \(N\) наблюдений с известной \(\sigma\): \[ \widehat\mu=\langle x\rangle, \qquad \sigma_{\widehat\mu}=\frac{\sigma}{\sqrt N}. \]

Профиль относительного правдоподобия точно квадратичен: \[ \boxed{ q(\mu) = \frac{(\mu-\langle x\rangle)^2} {\sigma^2/N}. } \]

Точки \(q=1\) расположены на расстоянии одной стандартной ошибки от максимума. Интерпретация таких порогов вне гауссова случая требует отдельного обоснования.

Пример: время жизни частицы

Пусть времена распада описываются экспоненциальным распределением: \[ T_i\sim\operatorname{Exp}(\tau), \qquad f(t\mid\tau)=\frac1\tau e^{-t/\tau}, \qquad t\ge0. \]

Функция правдоподобия: \[ L(\tau) = \prod_{i=1}^{N} \frac1\tau e^{-t_i/\tau} = \tau^{-N} \exp\!\left( -\frac1\tau\sum_{i=1}^{N}t_i \right). \]

Логарифм правдоподобия: \[ \ell(\tau) = -N\ln\tau -\frac1\tau\sum_{i=1}^{N}t_i +\mathrm{const}. \]

Условие максимума: \[ \frac{\partial\ell}{\partial\tau} = -\frac{N}{\tau} +\frac{1}{\tau^2}\sum_{i=1}^{N}t_i =0 \quad\Longrightarrow\quad \boxed{\widehat\tau_{\mathrm{ММП}}=\langle t\rangle.} \]

Интерактив: правдоподобие времени жизни

Границы и малые выборки

Для пуассоновского наблюдения \(n\): \[ q(\mu) = 2\left[ \mu-n+n\ln\frac n\mu \right], \] где при \(n=0\) последний член принимается равным нулю.

  • Форма не симметрична относительно \(\widehat\mu=n\).
  • При \(n=0\) максимум находится на границе \(\widehat\mu=0\).
  • Кривизна около максимума может плохо описывать дальние области.
  • Симметричная запись \(\widehat\mu\pm\sigma\) может быть физически бессмысленной.

Несколько параметров

Параметры интереса и мешающие параметры

Пусть модель зависит от \[ L(\theta,\boldsymbol\eta), \] где:

  • \(\theta\) — параметр интереса, который хотим измерить;
  • \(\boldsymbol\eta\) — мешающие параметры: фон, эффективность, калибровка, разрешение;
  • все параметры влияют на форму правдоподобия;
  • игнорирование неопределённости мешающих параметров обычно занижает ошибку \(\theta\).

Совместный максимум

Полная подгонка находит \[ (\widehat\theta,\widehat{\boldsymbol\eta}) = \underset{\theta,\boldsymbol\eta} {\operatorname{arg\,max}}\, L(\theta,\boldsymbol\eta). \]

  • Корреляции параметров образуют вытянутые области правдоподобия.
  • Изменение \(\theta\) может частично компенсироваться изменением \(\boldsymbol\eta\).
  • Поэтому неопределённость \(\theta\) определяется не только сечением правдоподобия при фиксированных мешающих параметрах.

Профильное правдоподобие

Для каждого фиксированного \(\theta\) повторно оптимизируем мешающие параметры: \[ \widehat{\widehat{\boldsymbol\eta}}(\theta) = \underset{\boldsymbol\eta}{\operatorname{arg\,max}}\, L(\theta,\boldsymbol\eta). \]

Профильное правдоподобие: \[ \boxed{ L_{\mathrm{p}}(\theta) = L\!\left( \theta, \widehat{\widehat{\boldsymbol\eta}}(\theta) \right). } \]

Двойная шляпка подчёркивает: мешающие параметры оптимизированы при фиксированном \(\theta\).

Модель для аплета

В аплете используется двумерная гауссова форма правдоподобия около максимума \((0,0)\): \[ q(\theta,\eta) = -2\ln\frac{L(\theta,\eta)}{L_{\max}} = \frac{\theta^2-2\rho\theta\eta+\eta^2}{1-\rho^2}, \qquad \frac{L(\theta,\eta)}{L_{\max}} = \exp\!\left[-\frac12 q(\theta,\eta)\right]. \]

Здесь \(\theta\) — интересующий параметр, \(\eta\) — мешающий параметр, а \(\rho\) задаёт их корреляцию: \[ \Sigma= \begin{pmatrix} 1 & \rho\\ \rho & 1 \end{pmatrix}, \qquad q= \begin{pmatrix}\theta & \eta\end{pmatrix} \Sigma^{-1} \begin{pmatrix}\theta\\ \eta\end{pmatrix}. \]

Профилирование выбирает лучшее значение \(\eta\) при каждом фиксированном \(\theta\): \[ \frac{\partial q}{\partial\eta}=0 \quad\Rightarrow\quad \widehat{\widehat\eta}(\theta)=\rho\theta, \qquad q_{\mathrm{prof}}(\theta)=q(\theta,\widehat{\widehat\eta}(\theta))=\theta^2. \]

Если вместо профилирования зафиксировать \(\eta=0\), получаем более крутую кривую: \[ q_{\eta=0}(\theta)=q(\theta,0)=\frac{\theta^2}{1-\rho^2}. \]

Интерактив: что делает профилирование

Цвет показывает совместное относительное правдоподобие \(L(\theta,\eta)/L_{\max}\). Белая линия проходит через лучшее значение \(\eta\) для каждого фиксированного \(\theta\).

Что показывает аплет

  1. При фиксированном \(\theta\) синяя вертикальная линия задаёт допустимые значения мешающего параметра \(\eta\).
  2. Оранжевая точка выбирает на этой линии максимум правдоподобия.
  3. Повторение для всех \(\theta\) образует белую профильную траекторию.
  4. Профильная кривая шире сечения при фиксированном \(\eta=0\).
  5. Чем сильнее корреляция, тем сильнее фиксирование мешающего параметра занижает неопределённость.

Пример: реакторные нейтрино

Игрушечная модель: два флейвора, вакуум, реакторные \(\bar\nu_e\): \[ P_{ee}(E; A, \Delta m^2) = 1 - A\, \sin^2\!\left( 1.267\,\Delta m^2\,\frac{L}{E_{\mathrm{ГэВ}}} \right), \qquad A=\sin^2 2\theta. \]

Ожидание в энергетическом бине: \[ \lambda_i(A,\Delta m^2) = C\,\phi(E_i)\sigma(E_i)P_{ee}(E_i;A,\Delta m^2)\Delta E_i, \qquad n_i\sim\operatorname{Pois}(\lambda_i). \]

Правдоподобие: \[ \ell(A,\Delta m^2) = \sum_i\left[n_i\ln\lambda_i(A,\Delta m^2)-\lambda_i(A,\Delta m^2)\right] \quad \text{с точностью до константы.} \]

Параметризация Huber–Mueller аппроксимирует логарифм спектра полиномом по \(E\): \[ \phi(E) = \sum_k f_k \exp\!\left[ \sum_{p=0}^{5} a_{kp} E^p \right], \qquad k\in\{^{235}\mathrm U,\,^{238}\mathrm U,\,^{239}\mathrm{Pu},\,^{241}\mathrm{Pu}\}. \]

\(E\) в МэВ; \(a_{kp}\) размерны, показатель экспоненты безразмерен.

изотоп \(f_k\) \(a_0\) \(a_1\) \(a_2\) \(a_3\) \(a_4\) \(a_5\)
\(^{235}\)U 0.58 4.367 -4.577 2.100 -0.5294 0.06186 -0.002777
\(^{238}\)U 0.07 0.4833 0.1927 -0.1283 -0.006762 0.002233 -0.0001536
\(^{239}\)Pu 0.30 4.757 -5.392 2.563 -0.6596 0.07820 -0.003536
\(^{241}\)Pu 0.05 2.990 -2.882 1.278 -0.3343 0.03905 -0.001754

В аплете фиксированы:

величина значение
расстояние \(L\) \(1.6\,\text{км}\)
диапазон энергий \(1.8\ldots 8.0\,\text{МэВ}\)
сетка по \(A=\sin^2 2\theta\) \(0\ldots 0.20\), 61 точка
сетка по \(\Delta m^2\) \((1.6\ldots 3.4)\cdot 10^{-3}\,\text{эВ}^2\), 61 точка

Пользователь меняет \(N_0\), \(A_0\), \(\Delta m^2_0\) и число энергетических бинов \(N_{\mathrm{bins}}\).

Для сечения обратного бета-распада используется простая форма \[ \sigma(E)\propto E_e p_e, \qquad E_e=E-1.293\,\text{МэВ}. \]

Нормировка \(C\) выбирается так, чтобы при \(A=0\) ожидалось \(N_0\) событий во всём диапазоне энергий.

Карта показывает \[ q(A,\Delta m^2) = -2\left[ \ell(A,\Delta m^2)-\ell_{\max} \right]. \]

  1. По заданным \((A_0,\Delta m^2_0)\) считаем \(\lambda_i^0\) во всех энергетических бинах.
  2. Генерируем псевдоданные: \(n_i\sim\operatorname{Pois}(\lambda_i^0)\).
  3. На сетке \((A,\Delta m^2)\) считаем \(\ell(A,\Delta m^2)\) и находим максимум \(\ell_{\max}\).
  4. Строим 2D-карту \(q(A,\Delta m^2)\). Белая область \(q\le 2.30\) — ориентир для совместной области \(1\sigma\) для двух параметров.
  5. Строим профили: \[ q_{\mathrm p}(A)=\min_{\Delta m^2}q(A,\Delta m^2), \qquad q_{\mathrm p}(\Delta m^2)=\min_A q(A,\Delta m^2). \]
  6. Одномерный ориентир \(1\sigma\): пересечение профильной кривой с \(q_{\mathrm p}=1\).

Интерактив: осцилляционные параметры

Что показывает пример с осцилляциями

  1. Даже в простой модели параметры не измеряются независимо: похожие спектры могут получаться при разных \((A,\Delta m^2)\).
  2. Совместная область показывает эту связь напрямую.
  3. Одномерная ошибка для \(A\) получается из профильной кривой \(q_{\mathrm p}(A)\), где для каждого \(A\) заново выбирается лучшее \(\Delta m^2\).
  4. Аналогично, ошибка для \(\Delta m^2\) получается из \(q_{\mathrm p}(\Delta m^2)\).
  5. Это тот же принцип, что в абстрактном аплете, но теперь профиль появляется из физической модели спектра.

Правдоподобие для примеров анализа данных

Ограничивающая информация

Мешающий параметр часто имеет внешнее измерение: \[ \eta_{\mathrm{всп}}\sim\mathrm N(\eta,\sigma_\eta^2). \]

Тогда полное правдоподобие содержит два множителя: \[ \boxed{ L_{\mathrm{полн}}(\theta,\eta) = L_{\mathrm{данные}}(\theta,\eta) \, L_{\mathrm{всп}}(\eta). } \]

  • Основные данные измеряют \(\theta\) и одновременно ограничивают \(\eta\).
  • Вспомогательное измерение добавляет независимую информацию об \(\eta\).
  • Слово «ограничение» не означает произвольный штраф: это часть вероятностной модели данных.

Расширенное правдоподобие

Если модель предсказывает и число событий \(\nu(\theta)\), и распределение их признаков \(f(x\mid\theta)\), используют \[ \boxed{ L_{\mathrm{расш}}(\theta) = \operatorname{Pois}(N\mid\nu(\theta)) \prod_{i=1}^{N} f(x_i\mid\theta). } \]

  • Пуассоновский множитель использует информацию о полном числе событий.
  • Произведение плотностей использует форму распределения.
  • Для смеси сигнала и фона в \(f(x\mid\theta)\) входят оба компонента.

Бинированное пуассоновское правдоподобие

Если данные представлены числами событий \(n_j\) в независимых бинах, а модель предсказывает \(\nu_j(\boldsymbol\theta)\): \[ \boxed{ L(\boldsymbol\theta) = \prod_j \operatorname{Pois} \!\left( n_j\mid\nu_j(\boldsymbol\theta) \right). } \]

Соответствующий логарифм правдоподобия: \[ \ell(\boldsymbol\theta) = \sum_j \left[ n_j\ln\nu_j(\boldsymbol\theta) - \nu_j(\boldsymbol\theta) - \ln n_j! \right]. \]

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

Если \[ \widehat\theta_{\mathrm{ММП}} = \underset{\theta}{\operatorname{arg\,max}}\,L(\theta), \] а \(\phi=g(\theta)\) — взаимно однозначное преобразование, то \[ \boxed{ \widehat\phi_{\mathrm{ММП}} = g(\widehat\theta_{\mathrm{ММП}}). } \]

  • Можно подгонять удобную параметризацию.
  • Результат затем переводится в физически интересующую величину.
  • Но ошибка при нелинейном преобразовании обычно становится асимметричной.

Практический алгоритм подгонки

  1. Явно записать вероятностную модель всех наблюдений.
  2. Определить параметры интереса, мешающие параметры и физические границы.
  3. Реализовать устойчивый логарифм правдоподобия.
  4. Найти совместный максимум несколькими начальными приближениями.
  5. Проверить градиенты, сходимость и положение относительно границ.
  6. Построить одномерные и двумерные сканы правдоподобия.
  7. Проверить процедуру на псевдоэкспериментах.
  8. Сравнить подогнанную модель с данными.

Что проверять на псевдоэкспериментах

Для известных истинных параметров \(\boldsymbol\theta_0\):

  • распределение \(\widehat{\boldsymbol\theta}\);
  • смещение и дисперсию оценки;
  • частоту попадания подгонки в границы;
  • устойчивость оптимизатора;
  • покрытие интервалов;
  • распределение статистик отношения правдоподобия;
  • влияние ошибочной модели.

Формула правдоподобия задаёт статистическую модель. Псевдоэксперименты проверяют, как ведёт себя вся процедура анализа.

Типичные проблемы

Проблема Симптом
локальный максимум результат зависит от начальной точки
параметр на границе несимметричное правдоподобие, неверная ошибка по кривизне
слабая идентифицируемость длинная плоская долина правдоподобия
сильная корреляция нестабильная подгонка, большая профильная ошибка
неверная модель смещение и плохое согласие с данными
численная нестабильность NaN, бесконечности, ложная сходимость

Минимизатор сообщил успех. Достаточно?

Нет. После численной оптимизации необходимо проверить:

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

Что важно запомнить

Резюме

  • Оценка — случайная функция данных; оценённое значение — её реализация.
  • Правдоподобие рассматривает вероятностную модель как функцию неизвестных параметров при фиксированных данных.
  • Оценка максимального правдоподобия выбирает параметры, максимизирующие правдоподобие.
  • Хорошую оценку характеризуют смещение, дисперсия, средний квадрат ошибки, состоятельность и эффективность.
  • Кривизна логарифма правдоподобия даёт локальную оценку точности только в регулярном почти гауссовом случае.
  • Профилирование учитывает компенсацию между параметром интереса и мешающими параметрами.
  • Подгонку в анализе данных необходимо проверять сканами, диагностикой оптимизатора и псевдоэкспериментами.