Качество согласия и статистическая значимость

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

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

После подгонки

Что мы уже умеем

  • К этой лекции мы умеем формулировать вероятностную модель данных и правило подгонки.

  • Для гауссовых ошибок мы минимизируем \[ \chi^2(\boldsymbol\theta) = \sum_{i=1}^{N} \frac{ \left[y_i-f(x_i;\boldsymbol\theta)\right]^2 }{\sigma_i^2}. \]

  • Минимум даёт оценку параметров: \[ \widehat{\boldsymbol\theta} = \operatorname{arg} \underset{\boldsymbol\theta}{\operatorname{min}}\, \chi^2(\boldsymbol\theta). \]

  • В окрестности минимума форма \(\chi^2\) даёт ковариационную матрицу параметров.

Что дальше?

  • Минимум существует почти всегда.

  • Даже плохая модель обычно имеет наилучшие параметры: \[ \chi^2_{\min} = \chi^2(\widehat{\boldsymbol\theta}). \]

  • Следующий важный вопрос:

    • насколько хорошо модель описывает наблюдения?
    • Это вопрос о качестве согласия модели с данными.

Пример: калибровочная прямая

  • Пусть измеряется отклик детектора в нескольких точках энергии: \[ y_i=a+b x_i+\varepsilon_i, \qquad \varepsilon_i\sim \mathrm N(0,\sigma_i^2). \]

  • Подгонка прямой всегда найдёт числа \(\widehat a\) и \(\widehat b\).

  • Но это ещё не означает, что данные действительно похожи на прямую.

  • Остатки после подгонки: \[ r_i = y_i-\widehat a-\widehat b x_i. \]

  • Именно остатки отвечают на вопрос о согласии.

Нормированные остатки

  • Если модель верна и ошибки оценены правильно, то величины \[ z_i = \frac{ y_i-f(x_i;\widehat{\boldsymbol\theta}) }{\sigma_i} \] должны быть типичными гауссовыми отклонениями.

  • Сумма квадратов нормированных остатков равна \[ \chi^2_{\min} = \sum_i z_i^2. \]

  • Большой вклад одной точки означает большой остаток в единицах её ошибки.

  • Большой суммарный вклад означает, что модель плохо описывает данные в целом.

Статистика

  • Теперь нас интересует не только \[ \widehat{\boldsymbol\theta}, \] но и \[ \chi^2_{\min}. \]

  • \(\chi^2_{\min}\) - это функция случайных переменных, значит сама - случайная переменная.

  • \(\chi^2_{\min}\) - это пример статистики

  • Статистика - измеримая числовая функция от выборки, не зависящая от неизвестных параметров распределения элементов выборки.

  • Вопрос лекции:

    • Если модель верна, какие значения \(\chi^2_{\min}\) являются обычными, а какие уже редкими?

  • Для ответа нужно знать распределение \(\chi^2_{\min}\).

Фольклор физиков

Из физического фольклора

  • Часто физики говорят: \[ \frac{\chi^2_{\min}}{\nu}\approx 1 \quad\Longrightarrow\quad \text{согласие хорошее}, \] где \[ \nu=N-m \] это число степеней свободы, \(N\) - число точек, \(m\) - число подогнанных параметров.

  • Это неправильно. Совершенно неправильно.

  • Вероятность одного и того же отношения зависит от числа степеней свободы \(\nu\) (об этом дальше).

Вопрос дня

  • Рассмотрим два случая:
    1. \(\chi^2_{\min}=12\) и \(\nu=10\) и
    2. \(\chi^2_{\min}=1200\) и \(\nu=1000\).
  • В обоих случаях \[ \chi^2_{\min}/\nu = 1.2 \]
  • В первом случае есть хорошее согласие модели с данными.
  • Во втором случае согласие модели с данными маловероятно.
  • Почему? Давайте разбираться!

Распределение \(\chi^2\)

Один нормированный остаток

  • Пусть \[ Z\sim \mathrm N(0,1). \]

  • Плотность \(Z\) равна \[ f_Z(z) = \frac{1}{\sqrt{2\pi}}e^{-z^2/2}. \]

  • Квадрат отклонения \[ T=Z^2\ge 0 \]

  • Большое значение \(T\) означает, что наблюдение оказалось далеко от ожидания в единицах своей ошибки.

Вывод: ветви и якобиан

  • Для \(t>0\) уравнение \[ z^2=t \] имеет два корня: \[ z=\sqrt t, \qquad z=-\sqrt t. \]

  • Для каждого корня якобиан равен \[ \left| \frac{dz}{dt} \right| = \frac{1}{2\sqrt t}. \]

  • Поэтому плотность вероятности \(T\) получается суммой вкладов двух ветвей:

Плотность вероятности \(f_T(t)\)

\[ \begin{aligned} f_T(t) &= f_Z(\sqrt t)\left|\frac{dz}{dt}\right|_{z=\sqrt t} + f_Z(-\sqrt t)\left|\frac{dz}{dt}\right|_{z=-\sqrt t}\\ &= \frac{f_Z(\sqrt t)+f_Z(-\sqrt t)}{2\sqrt t}. \end{aligned} \]

  • Так как \(f_Z(\sqrt t)=f_Z(-\sqrt t)\), \[ \boxed{ f_T(t) = \frac{1}{\sqrt{2\pi t}}e^{-t/2}, \qquad t>0. } \]

  • Это распределение \(\chi^2\) с одной степенью свободы.

Несколько независимых остатков

  • Пусть теперь есть \(\nu\) независимых стандартных гауссовых величин: \[ Z_1,\ldots,Z_\nu \sim \mathrm N(0,1). \]

  • Их совместная плотность пропорциональна \[ \exp\!\left[ -\frac12\sum_{i=1}^{\nu} z_i^2 \right]. \]

  • Обозначим \[ t = \sum_{i=1}^{\nu}z_i^2. \]

  • Это квадрат расстояния от центра в \(\nu\)-мерном пространстве нормированных отклонений.

Геометрия вывода

  • Совместная плотность зависит только от радиуса: \[ r^2=t. \]

  • Вероятность попасть в тонкий сферический слой радиуса \(r\) пропорциональна двум множителям: \[ \underbrace{e^{-r^2/2}}_{\text{плотность гаусса}} \times \underbrace{r^{\nu-1}\,dr}_{\text{объём слоя}}. \]

  • Переходим от \(r\) к \(t=r^2\): \[ r^{\nu-1}dr = \frac12 t^{\nu/2-1}dt. \]

  • Поэтому плотность по \(t\) имеет вид \[ f(t;\nu)\propto t^{\nu/2-1}e^{-t/2}. \]

Формула распределения

  • После нормировки получаем \[ \boxed{ f(t;\nu) = \frac{1}{2^{\nu/2}\Gamma(\nu/2)} t^{\nu/2-1}e^{-t/2}, \qquad t\ge 0. } \]

  • Это распределение \(\chi^2\) с \(\nu\) степенями свободы.

  • Здесь \[ \Gamma(a) = \int_0^\infty u^{a-1}e^{-u}\,du \] это гамма-функция.

  • Она обобщает факториал: \[ \Gamma(n)=(n-1)! \quad\text{для целого } n\ge 1, \qquad \Gamma(1/2)=\sqrt\pi. \]

Математическое ожидание

  • Используем плотность \(\chi^2_\nu\): \[ f(t;\nu) = \frac{1}{2^{\nu/2}\Gamma(\nu/2)} t^{\nu/2-1}e^{-t/2}. \]

  • Тогда \[ \begin{aligned} \mathbb E[T] &= \int_0^\infty t f(t;\nu)\,dt\\ &= \frac{1}{2^{\nu/2}\Gamma(\nu/2)} \int_0^\infty t^{\nu/2}e^{-t/2}\,dt. \end{aligned} \]

  • Подстановка \(t=2u\) даёт \[ \mathbb E[T] = 2\frac{\Gamma(\nu/2+1)}{\Gamma(\nu/2)} = \nu. \]

Дисперсия

  • Аналогично \[ \begin{aligned} \mathbb E[T^2] &= \frac{1}{2^{\nu/2}\Gamma(\nu/2)} \int_0^\infty t^{\nu/2+1}e^{-t/2}\,dt\\ &= 4\frac{\Gamma(\nu/2+2)}{\Gamma(\nu/2)} = \nu(\nu+2). \end{aligned} \]

  • Поэтому \[ \operatorname{Var}(T) = \mathbb E[T^2]-\mathbb E[T]^2 = \nu(\nu+2)-\nu^2 = 2\nu. \]

  • Итак, \[ \boxed{ \mathbb E[\chi^2_\nu]=\nu, \qquad \operatorname{Var}(\chi^2_\nu)=2\nu. } \]

  • Число степеней свободы задаёт не только среднее, но и ширину распределения.

Интерактив: распределение \(\chi^2\)

Если ошибки коррелированы

  • Если ошибки коррелированы, то сумма квадратов записывается через ковариационную матрицу: \[ \chi^2 = (\mathbf y-\boldsymbol\lambda)^{\mathsf T} V^{-1} (\mathbf y-\boldsymbol\lambda). \]

  • Вопрос:

    • Если между измерениями есть корреляции, не становится ли независимых величин меньше?

  • Ответ: нет, если ковариационная матрица невырождена.

  • Корреляции меняют координаты и масштабы, но не размерность пространства.

Диагонализация билинейной формы

  • Ковариационную матрицу можно диагонализовать: \[ V = Q\Lambda Q^{\mathsf T}, \qquad \Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_N). \]

  • Если все \(\lambda_i>0\), введём новые переменные \[ \mathbf z = \Lambda^{-1/2}Q^{\mathsf T} (\mathbf y-\boldsymbol\lambda). \]

  • Тогда \[ \chi^2 =\mathbf z^{\mathsf T}\mathbf z =\sum_{i=1}^{N}z_i^2. \]

  • Размерность осталась равной \(N\).

  • Точное уменьшение числа направлений возникает только при настоящем ограничении, когда одна из собственных дисперсий равна нулю и \(V\) становится вырожденной.

Почему после подгонки меньше степеней свободы

  • После перехода к независимым нормированным переменным данные живут в \(N\)-мерном пространстве.

  • Если параметры не подгоняются, то \(\chi^2\) — это квадрат длины вектора данных относительно заданной модели.

  • Если модель линейна по \(m\) подгоняемым параметрам, то эти параметры задают \(m\)-мерную поверхность возможных предсказаний.

  • Подгонка выбирает ближайшую точку этой поверхности.

  • Остаток — это часть вектора, ортогональная поверхности модели.

  • В линейной гауссовой задаче размерность пространства остатков равна \(N-m\).

  • Поэтому \[ \boxed{ \nu=N-m, } \] где \(m\) - число подгоняемых параметров.

  • Это не связано с корреляциями между измерениями. Это связано с тем, что сами данные использованы для выбора \(m\) параметров модели.

Пример: подгонка константы

  • Пусть измеряется одна и та же величина: \[ y_i=\mu+\varepsilon_i, \qquad \varepsilon_i\sim \mathrm N(0,\sigma^2). \]

  • Если \(\mu\) задана заранее, то \[ \chi^2(\mu) = \sum_{i=1}^{N} \frac{(y_i-\mu)^2}{\sigma^2} \sim \chi^2_N. \]

  • Если \(\mu\) оценивается по этим же данным, \[ \widehat\mu = \frac1N\sum_i y_i, \] то сумма остатков обязана быть равна нулю: \[ \sum_i (y_i-\widehat\mu) =0. \]

  • Это одно линейное условие на остатки.

Куда делась одна степень свободы

  • Обозначим вектор остатков \[ \mathbf r = (y_1-\widehat\mu,\ldots,y_N-\widehat\mu). \]

  • Условие \[ \sum_i r_i=0 \] означает, что \(\mathbf r\) лежит не во всём \(N\)-мерном пространстве, а в гиперплоскости размерности \(N-1\).

  • Геометрически подгонка константы вычитает из данных проекцию на направление \[ (1,1,\ldots,1). \]

  • Одно направление исчезло, поэтому \[ \chi^2_{\min}\sim \chi^2_{N-1}. \]

Всегда ли распределение \(\chi^2_{\mathrm{min}}\) это \(\chi^2_\nu\)?

  • Нет. Если

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

  • В таких случаях статистика может иметь распределение, отличающееся от \(\chi^2_\nu\).

  • В этих и похожих случаях следует получить распределение статистики \(\chi^2_{\mathrm{min}}\) при помощи псевдоэкспериментов.

\(p\)-уровень значимости

\(p\)-уровень значимости

  • Пусть после подгонки получено \[ \chi^2_{\mathrm{obs}}. \]

  • Если применимо распределение \(\chi^2_\nu\), то \(p\)-уровень равен \[ \boxed{ p = \int_{\chi^2_{\mathrm{obs}}}^{\infty} f(t;\nu)\,dt. } \]

  • Это вероятность получить такое же или большее несогласие с моделью при условии, что модель верна.

  • Чем меньше \(p\), тем менее типично наблюдённое значение \(\chi^2_{\mathrm{obs}}\).

\(p\)-уровень на графике

Алгоритм:

  1. Генерируем “теорию” \(T\) и псевдонаблюдения \(t_{obs}\).

  2. Вычисляем \(p\).

  3. Если псевдоданные и теория описываются одной моделью, значения \(p\) распределены равномерно.

import numpy as np

rng = np.random.default_rng(12345)

nu = 10
B_cal = 100_000
B_obs = 5_000

# Калибровка распределения статистики.
T_cal = rng.chisquare(df=nu, size=B_cal)
T_cal.sort()

# Независимые псевдонаблюдения.
T_obs = rng.chisquare(df=nu, size=B_obs)

# Для каждого наблюдённого t считаем p-уровень.
n_right = B_cal - np.searchsorted(T_cal, T_obs, side="left")
p = n_right / B_cal

hist, edges = np.histogram(p, bins=np.linspace(0, 1, 21), density=True)

Односторонний перевод \(p\to Z\)

  • Удобно переводить \(p\) в число сигма.

  • Односторонний перевод: \[ p = P(Z_{\mathrm G}\ge Z) = 1-\Phi(Z). \]

  • Значит, \[ \boxed{ Z=\Phi^{-1}(1-p). } \]

  • Здесь \(\Phi\) — функция распределения стандартной нормальной величины.

  • Ориентиры для одностороннего перевода:
\(Z\) \(p=1-\Phi(Z)\) \(q:\ P(\chi^2_1\ge q)=p\)
\(1\sigma\) \(1.59\cdot 10^{-1}\) \(1.99\)
\(2\sigma\) \(2.28\cdot 10^{-2}\) \(5.19\)
\(3\sigma\) \(1.35\cdot 10^{-3}\) \(10.27\)
\(4\sigma\) \(3.17\cdot 10^{-5}\) \(17.31\)
\(5\sigma\) \(2.87\cdot 10^{-7}\) \(26.34\)
  • Порог открытия в физике частиц обычно связывают с односторонним \(5\sigma\).

Двусторонний перевод

  • Для интервала у одного параметра нас обычно интересуют отклонения в обе стороны: \[ p = P(|Z_{\mathrm G}|\ge Z) = 2[1-\Phi(Z)]. \]

  • Если \[ \Delta\chi^2 = \chi^2(\theta)-\chi^2(\widehat\theta) \simeq \frac{(\theta-\widehat\theta)^2}{\sigma_\theta^2}, \] то \[ \Delta\chi^2 = Z^2. \]

  • Поэтому для одного параметра: \[ 1\sigma,\ 2\sigma,\ 3\sigma \quad\Longleftrightarrow\quad \Delta\chi^2=1,\ 4,\ 9. \]

  • Ориентиры для двустороннего перевода:
интервал двусторонний \(p\) \(\Delta\chi^2\)
\(1\sigma\) \(3.17\cdot 10^{-1}\) \(1\)
\(2\sigma\) \(4.55\cdot 10^{-2}\) \(4\)
\(3\sigma\) \(2.70\cdot 10^{-3}\) \(9\)
\(4\sigma\) \(6.33\cdot 10^{-5}\) \(16\)
\(5\sigma\) \(5.73\cdot 10^{-7}\) \(25\)
  • Не надо смешивать эти два соглашения: открытие обычно одностороннее, интервал для одного параметра обычно двусторонний.

Ответ дня. Почему нельзя писать \(\chi^2/\mathrm{ndf} = 1.2\)?

  • Вот почему
\(\chi^2_{\mathrm{obs}}\) \(\nu\) \(\chi^2_{\mathrm{obs}}/\nu\) \(z_{\chi^2}\) \(p\)
\(12\) \(10\) \(1.20\) \(0.45\) \(0.285\)
\(1200\) \(1000\) \(1.20\) \(4.47\) \(1.2\cdot 10^{-5}\)

  • где \[ z_{\chi^2} = \frac{\chi^2_{\mathrm{obs}}-\nu}{\sqrt{2\nu}}. \] отклонение от среднего в единицах “ошибки”

  • В обоих случаях отношение одинаковое.

  • Во втором случае значение \(\chi^2\) гораздо менее вероятно.

  • Правильно сообщать так:

\(\chi^2_{min}/ndf = 1200/1000\).

Как интерпретировать \(p\)?

  • Неправильно:

Модель неверна с вероятностью \(1-p\).

  • Правильно:

Вероятность получить такие или более несогласованные данные, при условии справедливости модели, равна p.

Пример: бозон Хиггса

  • По горизонтали перебирается масса гипотетического бозона.

  • В каждой точке считают локальный \(p\)-уровень фоновой гипотезы.

  • Это проверка гипотезы, а не \(p\)-уровень качества подгонки. Об этом в другой лекции.

  • Горизонтальные линии \(1\sigma,\ldots,5\sigma\) — это перевод \(p\leftrightarrow Z\).

  • Глубокий минимум около \(125\,\mathrm{GeV}\) означает: фон сам по себе редко даёт такие данные.

Как получить распределение статистики

Когда формулы недостаточно

Если известно, что \[ \chi^2_{\min}\sim\chi^2_{N-m}, \] то \(p\) находится аналитически.

Если это не гарантировано, не нужно угадывать новое число степеней свободы.

Нужно построить распределение той же статистики в той же процедуре анализа.

Идея:

Сколько раз верная модель сама рождает значение статистики не меньше наблюдённого?

Пример: малые пуассоновские бины

  • Пусть модель предсказывает ожидания \(\nu_i\) в \(N=20\) бинах.

  • Данные в каждом бине пуассоновские: \[ n_i \sim \operatorname{Pois}(\nu_i). \]

  • Статистика Пирсона: \[ \chi^2_{\mathrm P} = \sum_{i=1}^{N} \frac{(n_i-\nu_i)^2}{\nu_i}. \]

  • При больших \(\nu_i\) ожидаем приближение \[ \chi^2_{\mathrm P}\sim \chi^2_N. \]

  • Но если в части бинов ожидания малы, хвост распределения может отличаться.

Нужно моделировать ровно ту статистику, которую применяем к данным:

import numpy as np

nu = np.array([
    0.158, 0.295, 0.513, 0.833, 1.263,
    1.785, 2.355, 2.899, 3.330, 3.568,
    3.568, 3.330, 2.899, 2.355, 1.785,
    1.263, 0.833, 0.513, 0.295, 0.158,
])

rng = np.random.default_rng(12345)
B = 1_000_000

n = rng.poisson(nu, size=(B, len(nu)))
T = ((n - nu)**2 / nu).sum(axis=1)

T_obs = 29.8
p_model = np.mean(T >= T_obs)

print(p_model)

# Для T_obs = 29.8:
# p_model ≈ 0.11.
# Для сравнения: χ²_20-приближение даёт p ≈ 0.073.

  • \(\nu_i\) — не случайные числа, а ожидания модели в бинах.

  • Случайны только наблюдённые счёты: \[ n_i\sim\operatorname{Pois}(\nu_i). \]

  • В этом примере специально есть бины с малыми ожиданиями.

бин \(\nu_i\) бин \(\nu_i\) бин \(\nu_i\) бин \(\nu_i\) бин \(\nu_i\)
1 0.158 5 1.263 9 3.330 13 2.899 17 0.833
2 0.295 6 1.785 10 3.568 14 2.355 18 0.513
3 0.513 7 2.355 11 3.568 15 1.785 19 0.295
4 0.833 8 2.899 12 3.330 16 1.263 20 0.158

Интерактив: распределение статистики

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

  • В линейной гауссовой задаче с фиксированной ковариацией: \[ \chi^2_{\min}\sim\chi^2_{N-m}. \]

  • Подгоняемые параметры уменьшают размерность пространства остатков.

  • Если условия для формулы сомнительны, распределение статистики получают псевдоэкспериментами.

  • В результате всегда нужно сообщать: \[ \chi^2_{\min}, \qquad \nu\ \text{или способ калибровки}, \qquad p. \]