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

Аннотация

Подгонка даёт параметры модели и их ошибки. В этой главе проверим, насколько сама модель согласуется с данными. Выведем распределение χ², разберём число степеней свободы и научимся вычислять p-значение. Вернёмся к графикам поиска бозона Хиггса из первой главы и выясним, что означают линии «числа сигма». Для малых пуассоновских счётов получим распределение статистики псевдоэкспериментами.

13.1 Физическая задача: можно ли пользоваться калибровкой?

Мы откалибровали детектор по источникам с известными энергиями. В каждой точке измерили среднюю амплитуду сигнала и подогнали прямую

\[ y_i=a+bE_i+\varepsilon_i, \qquad \varepsilon_i\sim\mathcal N(0,\sigma_i^2). \tag{13.1}\]

Здесь \(E_i\) — энергия в МэВ, \(y_i\) — амплитуда в каналах, \(a\) — смещение нуля, \(b\) — коэффициент преобразования в каналах на МэВ. Энергии считаем известными точно, ошибки амплитуд — независимыми и гауссовыми, а \(\sigma_i\) — известными до подгонки.

Программа нашла \(\widehat a\), \(\widehat b\) и их ковариацию. Можно теперь пересчитывать амплитуды в энергии? Только если прямая описывает отклик детектора с нужной нам точностью. Электроника может давать нелинейность; при больших амплитудах возможно насыщение. Минимизатор об этом ничего не знает. Он найдёт наиболее подходящую прямую для тех точек, которые ему передали.

Возьмём два учебных результата. В первой калибровке было \(N=12\) точек, и после подгонки двух параметров получилось \(\chi^2_{\min}=12\). Во второй — \(N=1002\) точки и \(\chi^2_{\min}=1200\). Числа степеней свободы равны соответственно \(10\) и \(1000\). В обоих случаях

\[ \frac{\chi^2_{\min}}{\nu}=1.2. \tag{13.2}\]

В физическом фольклоре отношение, близкое к единице, часто принимают за достаточную проверку фита. Наши два примера показывают, чего в такой проверке не хватает. При верной модели вероятность получить значение \(\chi^2_{\min}\) не меньше указанного равна примерно \(0.285\) в первом случае и \(1.23\cdot10^{-5}\) во втором.

Разница больше четырёх порядков. Давайте получим эти числа и разберёмся, что именно они позволяют сказать о калибровке.

13.2 Остатки и статистика качества согласия

После подгонки остатки равны

\[ r_i=y_i-\widehat a-\widehat bE_i, \qquad u_i=\frac{r_i}{\sigma_i}, \qquad \chi^2_{\min}=\sum_{i=1}^N u_i^2. \tag{13.3}\]

Остаток \(r_i\) измеряется в каналах, нормированный остаток \(u_i\) безразмерен. Точка с \(u_i=3\) находится на расстоянии трёх своих исходных ошибок от подогнанной прямой и вносит в сумму квадратов девятку.

Но сами \(u_i\) после подгонки уже зависимы. Для прямой нормальные уравнения требуют

\[ \sum_i\frac{r_i}{\sigma_i^2}=0, \qquad \sum_i\frac{E_i r_i}{\sigma_i^2}=0. \tag{13.4}\]

Мы использовали эти же точки, чтобы выбрать прямую. Поэтому остатки связаны двумя условиями, а их дисперсии меньше дисперсий исходных ошибок. Нельзя считать все \(u_i\) независимыми величинами \(\mathcal N(0,1)\). В главе 11 это было видно из матрицы проекции.

Величина \(\chi^2_{\min}\) — статистика, то есть вычисляемая по данным числовая функция. При повторении калибровки изменятся амплитуды, оценки \(\widehat a\), \(\widehat b\) и минимум суммы квадратов. Полученное число нужно сравнить с распределением этой статистики в повторяемых опытах при верной модели [1].

При этом \(\chi^2_{\min}\) и \(\Delta\chi^2\) из прошлой главы отвечают на разные вопросы. Минимум характеризует остаточное расхождение данных с моделью после подгонки. Приращение относительно минимума показывает, как меняется согласие при смещении параметров. Распределения у этих двух величин тоже разные.

13.3 Откуда берётся распределение χ²

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

Начнём со стандартной гауссовой величины:

\[ G\sim\mathcal N(0,1), \qquad \varphi(g)=\frac{1}{\sqrt{2\pi}}e^{-g^2/2}, \qquad T=G^2. \tag{13.5}\]

Букву \(Z\) оставим для значимости, которая появится ниже. Для \(t>0\) уравнение \(g^2=t\) имеет два корня: \(\sqrt t\) и \(-\sqrt t\). При замене переменной нужно учесть оба, как мы делали в главе 5.

Для каждой ветви \(|dg/dt|=1/(2\sqrt t)\). Поэтому

\[ \begin{aligned} f_T(t) &=\frac{\varphi(\sqrt t)+\varphi(-\sqrt t)}{2\sqrt t}\\ &=\frac{1}{\sqrt{2\pi t}}e^{-t/2}, \qquad t>0. \end{aligned} \tag{13.6}\]

Это плотность распределения χ² с одной степенью свободы. Отрицательных значений у \(T\) нет. При приближении к нулю плотность растёт без границы, но её интеграл остаётся конечным. Бесконечная плотность в одной точке не означает бесконечной вероятности.

13.3.2 Сумма квадратов и геометрия

Теперь возьмём \(\nu\) независимых величин \(G_i\sim\mathcal N(0,1)\):

\[ T=\sum_{i=1}^{\nu}G_i^2, \qquad T\sim\chi^2_\nu. \tag{13.7}\]

Их совместная плотность пропорциональна \(\exp(-\sum_i g_i^2/2)\). Она зависит только от расстояния до начала координат. В двумерном случае мы уже интегрировали такую плотность по кругам; теперь число измерений произвольно.

Для сферического слоя радиуса \(r\) вероятность пропорциональна

\[ e^{-r^2/2}r^{\nu-1}\,dr. \]

Здесь \(r^{\nu-1}dr\) даёт зависимость объёма слоя от радиуса. При переходе к \(t=r^2\) получаем \(r^{\nu-1}dr=\tfrac12t^{\nu/2-1}dt\). Нормировка приводит к плотности

\[ f(t;\nu)= \frac{t^{\nu/2-1}e^{-t/2}}{2^{\nu/2}\Gamma(\nu/2)}, \qquad t>0. \tag{13.8}\]

Гамма-функция определяется интегралом

\[ \Gamma(a)=\int_0^\infty u^{a-1}e^{-u}\,du, \qquad a>0. \tag{13.9}\]

Из интегрирования по частям следует \(\Gamma(a+1)=a\Gamma(a)\). В частности, \(\Gamma(n)=(n-1)!\) для целого \(n\ge1\) и \(\Gamma(1/2)=\sqrt\pi\). Подстановка \(\nu=1\) в 13.8 возвращает результат 13.6.

13.3.3 Среднее, ширина и число степеней свободы

Найдём первые два момента. В обоих интегралах достаточно сделать замену \(t=2u\):

\[ \begin{aligned} \mathbb E[T] &=2\frac{\Gamma(\nu/2+1)}{\Gamma(\nu/2)}=\nu,\\ \mathbb E[T^2] &=4\frac{\Gamma(\nu/2+2)}{\Gamma(\nu/2)}=\nu(\nu+2). \end{aligned} \tag{13.10}\]

Отсюда

\[ \operatorname{Var}(T)=2\nu, \qquad \mathbb E[T/\nu]=1, \qquad \sigma_{T/\nu}=\sqrt{\frac{2}{\nu}}. \tag{13.11}\]

Последняя формула объясняет проблему с правилом «χ² на степень свободы около единицы». При \(\nu=10\) стандартное отклонение отношения равно \(0.447\), при \(\nu=1000\) — \(0.0447\). То же превышение единицы на \(0.2\) во втором случае в десять раз больше в единицах ширины распределения.

13.3.4 Интерактив: форма распределения

Меняйте число степеней свободы \(\nu\). Сплошная вертикальная линия показывает среднее, пунктир — расстояние \(\sqrt{2\nu}\) от него. При \(\nu=1\) верхушка плотности у нуля выходит за пределы рисунка; сама формула остаётся той же.

viewof chiDistNu = Inputs.range([1, 80], {value: 10, step: 1, label: "Степени свободы ν"})
chiDistLogGamma = (z) => {
  const coefficients = [676.5203681218851, -1259.1392167224028,
    771.32342877765313, -176.61502916214059, 12.507343278686905,
    -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
  z -= 1;
  let sum = 0.99999999999980993;
  for (let i = 0; i < coefficients.length; i++) sum += coefficients[i] / (z + i + 1);
  const t = z + 7.5;
  return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(sum);
}
chiDistPdf = (t, nu) => {
  if (t < 0) return 0;
  if (t === 0) return nu < 2 ? Infinity : nu === 2 ? 0.5 : 0;
  const a = nu / 2;
  return Math.exp((a - 1) * Math.log(t) - t / 2 - a * Math.log(2) - chiDistLogGamma(a));
}
chiDistSf = (t, nu) => {
  if (!(nu > 0) || Number.isNaN(t)) return NaN;
  if (t <= 0) return 1;
  if (t === Infinity) return 0;
  const a = nu / 2, x = t / 2;
  const scale = Math.exp(a * Math.log(x) - x - chiDistLogGamma(a));
  const eps = 3e-14, tiny = 1e-300;
  if (x < a + 1) {
    let term = 1 / a, sum = term;
    for (let n = 1; n < 10000; n++) {
      term *= x / (a + n);
      sum += term;
      if (Math.abs(term) < eps * Math.abs(sum)) return Math.max(0, Math.min(1, 1 - scale * sum));
    }
  } else {
    let b = x + 1 - a, c = 1 / tiny, d = 1 / b, h = d;
    for (let n = 1; n < 10000; n++) {
      const an = -n * (n - a);
      b += 2;
      d = an * d + b;
      if (Math.abs(d) < tiny) d = tiny;
      c = b + an / c;
      if (Math.abs(c) < tiny) c = tiny;
      d = 1 / d;
      const delta = d * c;
      h *= delta;
      if (Math.abs(delta - 1) < eps) return Math.max(0, Math.min(1, scale * h));
    }
  }
  throw new Error("Не сошлось вычисление хвоста χ²");
}
chiDistProbabilityText = p => p === 0 ? "0" : p < 0.001 ? p.toExponential(3) : p.toFixed(5)
chiDistPlotLimit = nu => nu === 1 ? 1.5 : 1.12 * chiDistPdf(Math.max(0, nu - 2), nu)
chiDistGrid = (nu, xmax, count = 1000) => Array.from({length: count}, (_, i) => {
  const t = Math.max(1e-8, xmax * (i / (count - 1)) ** 2);
  return {t, density: chiDistPdf(t, nu)};
})
chiDistXmax = Math.max(8, chiDistNu + 5 * Math.sqrt(2 * chiDistNu))
chiDistData = chiDistGrid(chiDistNu, chiDistXmax)
chiDistSigma = Math.sqrt(2 * chiDistNu)
chiDistPlot = Plot.plot({
  width: 760, height: 340, marginLeft: 58, marginBottom: 48,
  style: {fontSize: "14px"},
  x: {label: "t", domain: [0, chiDistXmax], grid: true},
  y: {label: "f(t; ν)", domain: [0, chiDistPlotLimit(chiDistNu)], grid: true},
  marks: [
    Plot.areaY(chiDistData, {x: "t", y: "density", fill: "#367ab5", fillOpacity: 0.14, clip: true}),
    Plot.line(chiDistData, {x: "t", y: "density", stroke: "#367ab5", strokeWidth: 2.5, clip: true}),
    Plot.ruleX([chiDistNu], {stroke: "#c44536", strokeWidth: 2}),
    Plot.ruleX([chiDistNu - chiDistSigma, chiDistNu + chiDistSigma].filter(t => t > 0),
      {stroke: "#536773", strokeDasharray: "5,4"})
  ]
})
html`<div class="book-applet-panel">
  ${chiDistPlot}
  <div class="book-applet-summary">
    <span>ν = ${chiDistNu}</span>
    <span>E[T] = ${chiDistNu}</span>
    <span>σT = ${chiDistSigma.toFixed(3)}</span>
    <span>σ(T/ν) = ${Math.sqrt(2 / chiDistNu).toFixed(3)}</span>
    <span>${chiDistNu === 1 ? "При t → 0 плотность неограниченно растёт; график обрезан сверху." : ""}</span>
  </div>
</div>`

При \(\nu=10\) среднее равно \(10\), стандартное отклонение — \(\sqrt{20}\). Плотность асимметрична. Отметки среднего и ширины не являются границами интервала с заданной вероятностью.

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

С ростом \(\nu\) распределение становится более симметричным. В центральной области работает гауссово приближение \((T-\nu)/\sqrt{2\nu}\approx\mathcal N(0,1)\). Для далёкого хвоста точность этого приближения нужно проверять: именно там нас могут интересовать вероятности порядка \(10^{-5}\) и меньше.

13.4 Сколько степеней свободы осталось после фита?

13.4.1 Коррелированные ошибки

Пусть данные гауссовы, их средние \(\boldsymbol\lambda\) заданы заранее, а ковариация \(V_y\) известна и положительно определена:

\[ \chi^2= (\mathbf y-\boldsymbol\lambda)^{\mathsf T} V_y^{-1}(\mathbf y-\boldsymbol\lambda). \tag{13.12}\]

Между измерениями могут быть сильные корреляции. Уменьшилось ли из-за этого число независимых направлений? Проверим заменой координат.

Диагонализуем ковариацию и нормируем каждую координату на её ошибку:

\[ \begin{aligned} V_y&=QDQ^{\mathsf T}, \qquad D=\operatorname{diag}(d_1,\ldots,d_N),\\ \mathbf g&=D^{-1/2}Q^{\mathsf T} (\mathbf y-\boldsymbol\lambda). \end{aligned} \tag{13.13}\]

Матрица \(Q\) ортогональна, все \(d_i>0\). В новых координатах ковариация единичная. Для совместно гауссовых данных это означает независимость компонент, и

\[ \chi^2=\sum_{i=1}^{N}g_i^2\sim\chi^2_N. \tag{13.14}\]

Осталось ровно \(N\) направлений. Корреляции изменили их ориентацию и масштабы. При вырожденной ковариации появляются точные ограничения: тогда сначала выделяют допустимое подпространство. Обычной обратной матрицы в формуле 13.12 уже нет.

13.4.2 Подгонка выбирает проекцию

Для линейной модели \(X\boldsymbol\theta\) с \(m\) независимыми параметрами предсказания образуют \(m\)-мерное линейное подпространство. Здесь \(m\) — ранг матрицы модели; два одинаковых столбца не дают двух независимых параметров.

После преобразования 13.13 подгонка выбирает ближайшую точку подпространства модели. Остаток ортогонален ему и имеет \(N-m\) независимых компонент. Если модель верна, все эти компоненты стандартные гауссовы. Поэтому [1]

\[ \chi^2_{\min}\sim\chi^2_\nu, \qquad \nu=N-m. \tag{13.15}\]

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

13.4.3 Куда уходит одна степень свободы

Для измерений одной величины с одинаковыми известными ошибками

\[ Y_i\sim\mathcal N(\mu,\sigma^2), \qquad \widehat\mu=\overline Y. \tag{13.16}\]

До подгонки сумма \(\sum_i(Y_i-\mu)^2/\sigma^2\) имеет \(N\) степеней свободы. После подгонки \(\sum_i(Y_i-\overline Y)=0\): свободных остатков стало на один меньше.

Это можно записать точным равенством

\[ \sum_i\frac{(Y_i-\mu)^2}{\sigma^2} = \sum_i\frac{(Y_i-\overline Y)^2}{\sigma^2} +\frac{N(\overline Y-\mu)^2}{\sigma^2}. \tag{13.17}\]

Последнее слагаемое — квадрат одной стандартной гауссовой величины. В гауссовой модели оно независимо от первой суммы справа. Подгонка убрала одно направление, оставив \(\chi^2_{N-1}\).

У калибровочной прямой два независимых параметра, если энергии источников не все одинаковы. Поэтому в наших двух калибровках \(\nu=12-2=10\) и \(\nu=1002-2=1000\).

13.4.4 А если есть калибровочные ограничения?

В прошлой главе мы включали в правдоподобие вспомогательные измерения. Их нужно учитывать и при проверке качества согласия. В линейной гауссовой задаче с \(N\) основными и \(K\) независимыми вспомогательными данными после подгонки \(r\) независимых параметров остаётся \(N+K-r\) степеней свободы.

Например, общая неизвестная поправка калибровки добавляет один параметр, но независимое измерение этой поправки добавляет одно наблюдение. Просто вычесть параметр из прежнего \(N\) было бы неверно. Этот подсчёт предполагает, что в повторяемом опыте флуктуируют и основные, и вспомогательные данные. Для другой процедуры распределение минимума проверяют отдельно.

13.5 p-значение: какую вероятность мы вычисляем?

Обозначим проверяемую гипотезу через \(H_0\). В калибровочном примере она утверждает, что отклик линейный, ошибки гауссовы и их ковариация задана верно. Параметры прямой подгоняются по данным.

Пусть получено \(t_{\mathrm{obs}}=\chi^2_{\min}\). Сумма квадратов тем больше, чем сильнее остаточное расхождение. Поэтому берём правый хвост:

\[ p=P_{H_0}(T\ge t_{\mathrm{obs}}) =\int_{t_{\mathrm{obs}}}^{\infty}f(t;\nu)\,dt. \tag{13.18}\]

Спросим словами: если повторять эксперимент при верной \(H_0\), в какой доле опытов получится сумма квадратов хотя бы такого размера?

Вопрос о вероятности самой гипотезы здесь не ставился. Поэтому из \(p=0.03\) нельзя получить утверждение «вероятность модели равна 3%» или «модель неверна с вероятностью 97%». Число относится к распределению возможных результатов при заданной гипотезе.

Если нужен формальный критерий отклонения, заранее выбирают уровень значимости \(\alpha\) и отклоняют \(H_0\) при \(p\le\alpha\). Для точной непрерывной калибровки доля ложных отклонений верной гипотезы будет равна \(\alpha\). Значение \(\alpha\) задаётся правилом проверки, а \(p\) вычисляется по конкретным данным.

13.5.1 Интерактив: площадь хвоста

В апплете можно менять \(\nu\) и полученное значение \(t_{\mathrm{obs}}\). При движении вертикальной линии вправо закрашенная площадь уменьшается. Вероятность вычисляется по функции распределения; точность числа не зависит от густоты точек на графике.

viewof pLevelNu = Inputs.range([1, 120], {value: 50, step: 1, label: "Степени свободы ν"})
viewof pLevelObs = Inputs.range([0, 180], {value: 70, step: 0.5, label: "Значение t_obs"})
pLevelXmax = Math.max(12, pLevelObs + 5, pLevelNu + 6 * Math.sqrt(2 * pLevelNu))
pLevelData = chiDistGrid(pLevelNu, pLevelXmax)
pLevelTailData = [{t: Math.max(pLevelObs, 1e-8), density: chiDistPdf(Math.max(pLevelObs, 1e-8), pLevelNu)},
  ...pLevelData.filter(d => d.t > pLevelObs)]
pLevelTailArea = chiDistSf(pLevelObs, pLevelNu)
pLevelPlot = Plot.plot({
  width: 760, height: 340, marginLeft: 58, marginBottom: 48,
  style: {fontSize: "14px"},
  x: {label: "t", domain: [0, pLevelXmax], grid: true},
  y: {label: "f(t; ν)", domain: [0, chiDistPlotLimit(pLevelNu)], grid: true},
  marks: [
    Plot.areaY(pLevelTailData, {x: "t", y: "density", fill: "#c44536", fillOpacity: 0.35, clip: true}),
    Plot.line(pLevelData, {x: "t", y: "density", stroke: "#367ab5", strokeWidth: 2.5, clip: true}),
    Plot.ruleX([pLevelObs], {stroke: "#c44536", strokeWidth: 2})
  ]
})
html`<div class="book-applet-panel">
  ${pLevelPlot}
  <div class="book-applet-summary">
    <span>ν = ${pLevelNu}</span>
    <span>t_obs = ${pLevelObs.toFixed(1)}</span>
    <span>p = ${chiDistProbabilityText(pLevelTailArea)}</span>
    <span>${pLevelNu === 1 ? "Плотность у нуля обрезана сверху. Вероятность вычислена по полному хвосту." : ""}</span>
  </div>
</div>`

Для \(\nu=50\) и \(t_{\mathrm{obs}}=70\) площадь правого хвоста равна \(p\simeq0.0324\). Закрашена вероятность получить \(T\ge70\) при этой модели.

p-значение как площадь хвоста

Большое \(p\) означает, что выбранная статистика не обнаружила сильного несогласия. Оно не устанавливает истинность модели: другая модель тоже может описывать эти данные, а статистики может не хватить для различения моделей.

Очень маленькое \(\chi^2\) тоже даёт повод проверить расчёт. Возможно, ошибки завышены или часть данных была учтена несколько раз как независимая. Но иногда малое значение возникает просто из-за флуктуации. Правохвостовое \(p\), близкое к единице, само по себе не является критерием для этой отдельной проверки.

13.5.2 Как распределены сами p-значения?

Представим много опытов, в которых \(H_0\) верна. Для каждого вычислим \(T\) и его p-значение. Если функция распределения \(F_0\) непрерывна, то \(p=1-F_0(T)\). Преобразование переменной даёт

\[ P_{H_0}(p\le u)=u, \qquad 0\le u\le1. \tag{13.19}\]

Плотность p-значений равна единице на отрезке \([0,1]\). Значения около \(0.5\) ничем не выделены. При пороге \(\alpha=0.05\) примерно каждый двадцатый опыт с верной моделью попадёт в область отклонения.

В апплете одна выборка из \(\chi^2_{10}\) задаёт численную калибровку, а независимая вторая выборка изображает повторные эксперименты. Число опытов во второй выборке можно менять. Для каждого результата подсчитывается правый хвост калибровочной выборки.

viewof pUniformCount = Inputs.range([500, 20000], {value: 5000, step: 500, label: "Число проверяемых опытов"})
pUniformRandomFactory = () => {
  let state = 246813579;
  return () => {state = (1664525 * state + 1013904223) >>> 0; return state / 4294967296;};
}
pUniformNormal = rand => Math.sqrt(-2 * Math.log(Math.max(rand(), 1e-12))) * Math.cos(2 * Math.PI * rand())
pUniformChi2Sample = (count, nu, rand) => Array.from({length: count}, () => {
  let sum = 0;
  for (let i = 0; i < nu; i++) sum += pUniformNormal(rand) ** 2;
  return sum;
})
pUniformSamples = {
  const rand = pUniformRandomFactory();
  return {calibration: pUniformChi2Sample(20000, 10, rand).sort((a,b) => a-b),
    observed: pUniformChi2Sample(20000, 10, rand)};
}
pUniformLowerBound = (values, t) => {
  let lo = 0, hi = values.length;
  while (lo < hi) {
    const mid = Math.floor((lo + hi) / 2);
    if (values[mid] < t) lo = mid + 1;
    else hi = mid;
  }
  return lo;
}
pUniformAllValues = pUniformSamples.observed.map(t => {
  const values = pUniformSamples.calibration, b = values.length;
  return (b - pUniformLowerBound(values, t) + 1) / (b + 1);
})
pUniformValues = pUniformAllValues.slice(0, pUniformCount)
pUniformHistogram = {
  const bins = new Array(20).fill(0);
  for (const p of pUniformValues) bins[Math.min(19, Math.floor(20 * p))]++;
  return bins.map((n, i) => ({x1:i/20, x2:(i+1)/20, density:20*n/pUniformCount}));
}
pUniformPlot = Plot.plot({
  width:760, height:340, marginLeft:58, marginBottom:48,
  style:{fontSize:"14px"},
  x:{label:"p",domain:[0,1],grid:true},
  y:{label:"плотность",domain:[0,Math.max(1.45,1.12*Math.max(...pUniformHistogram.map(d=>d.density)))],grid:true},
  marks:[
    Plot.rectY(pUniformHistogram,{x1:"x1",x2:"x2",y:"density",fill:"#367ab5",fillOpacity:0.6}),
    Plot.ruleY([1],{stroke:"#c44536",strokeWidth:2,strokeDasharray:"5,4"})
  ]
})
html`<div class="book-applet-panel">
  ${pUniformPlot}
  <div class="book-applet-summary">
    <span>20 000 опытов для калибровки</span>
    <span>${pUniformCount} проверяемых опытов</span>
    <span>Доля p ≤ 0.05: ${(pUniformValues.filter(p=>p<=0.05).length/pUniformCount).toFixed(4)}</span>
    <span>p = (k + 1)/(20 000 + 1)</span>
  </div>
</div>`

20 000 опытов задают калибровку; для 5000 независимых опытов построена гистограмма p-значений. Пунктир показывает единичную плотность. Отклонения столбцов связаны с конечным размером обеих выборок.

Распределение p-значений

При конечной калибровочной выборке равномерность лишь приближённая. Для дискретной статистики точное распределение p-значений также ступенчатое. При определении через \(P(T\ge t_{\mathrm{obs}})\) получается консервативное правило:

\[ P_{H_0}(p\le\alpha)\le\alpha. \tag{13.20}\]

Например, при малом пуассоновском счёте может вообще не существовать порога, дающего ровно 5% ложных отклонений. Это свойство дискретных данных, которое следует учитывать при чтении гистограммы p-значений.

13.6 Возвращаемся к двум калибровкам

Теперь вычислим интеграл 13.18 для результатов из начала главы. Для сравнения приведём и отклонение от среднего в единицах ширины:

\[ z_{\chi^2}=\frac{t_{\mathrm{obs}}-\nu}{\sqrt{2\nu}}. \tag{13.21}\]

Таблица 13.1: Два результата калибровки с одинаковым отношением χ²/ν

\(t_{\mathrm{obs}}\)	\(\nu\)	\(t_{\mathrm{obs}}/\nu\)	\(z_{\chi^2}\)	\(p\)
\(12\)	\(10\)	\(1.20\)	\(0.447\)	\(0.2851\)
\(1200\)	\(1000\)	\(1.20\)	\(4.472\)	\(1.226\cdot10^{-5}\)

Первая калибровка не вызывает возражений по этой проверке. Во второй остатки слишком велики для принятой модели ошибок. Нужно смотреть, где именно возникло расхождение: кривизна отклика, отдельная точка, зависимость ошибки от энергии, забытая общая систематика.

Само p-значение не выбирает причину. Оно относится ко всей проверяемой модели, включая ошибки. Если увеличить все \(\sigma_i\), то χ² уменьшится. Но подбирать ошибки так, чтобы отношение стало единицей, — значит менять исходную постановку по результату этой же проверки.

Поэтому запись результата должна сохранять оба числа:

\[ \chi^2_{\min}/\mathrm{ndf}=1200/1000, \qquad p=1.226\cdot10^{-5}. \tag{13.22}\]

Сокращение ndf означает число степеней свободы. Одна запись «χ²/ndf = 1.2» теряет информацию, от которой зависит вывод.

13.7 Что означает «число сигма»?

13.7.1 Один хвост

Малую хвостовую вероятность удобно перевести на знакомую гауссову шкалу. Для стандартной гауссовой величины \(G\) определим число \(Z\) условием

\[ p=P(G\ge Z)=1-\Phi(Z), \qquad Z=\Phi^{-1}(1-p). \tag{13.23}\]

Функция \(\Phi\) — функция распределения \(\mathcal N(0,1)\). Это соглашение о представлении вероятности. Исходные данные могут быть пуассоновскими, а статистика — иметь распределение χ². Гауссовой здесь служит шкала, на которую мы переводим уже найденное \(p\).

Для второго результата в таблице 13.1 односторонний перевод даёт \(Z\simeq4.22\). Число \(z_{\chi^2}=4.472\) в той же таблице другое: оно получено вычитанием среднего и делением на стандартное отклонение. У асимметричного распределения эта операция не воспроизводит точную вероятность хвоста.

viewof oneSideGaussZ = Inputs.range([-3, 6], {value: 2, step: 0.1, label: "Граница G: Z"})
oneSideGaussPdf = g => Math.exp(-g*g/2)/Math.sqrt(2*Math.PI)
oneSideGaussP = oneSideGaussZ >= 0 ? 0.5*chiDistSf(oneSideGaussZ**2,1) : 1-0.5*chiDistSf(oneSideGaussZ**2,1)
oneSideGaussXmax = Math.max(4.5, oneSideGaussZ+1)
oneSideGaussData = Array.from({length:1201},(_,i)=>{
  const g = -oneSideGaussXmax + 2*oneSideGaussXmax*i/1200;
  return {g, density:oneSideGaussPdf(g)};
})
oneSideGaussRight = [{g:oneSideGaussZ,density:oneSideGaussPdf(oneSideGaussZ)},...oneSideGaussData.filter(d=>d.g>oneSideGaussZ)]

oneSideGaussPlot = Plot.plot({
  width:760,height:340,marginLeft:58,marginBottom:48,
  style:{fontSize:"14px"},
  x:{label:"g",domain:[-oneSideGaussXmax,oneSideGaussXmax],grid:true},
  y:{label:"φ(g)",domain:[0,0.43],grid:true},
  marks:[
    Plot.areaY(oneSideGaussRight,{x:"g",y:"density",fill:"#c44536",fillOpacity:0.4}),

    Plot.line(oneSideGaussData,{x:"g",y:"density",stroke:"#367ab5",strokeWidth:2.5}),
    Plot.ruleX([oneSideGaussZ],{stroke:"#c44536",strokeWidth:2})
  ]
})
html`<div class="book-applet-panel">
  ${oneSideGaussPlot}
  <div class="book-applet-summary">
    <span>Z = ${oneSideGaussZ.toFixed(1)}</span>
    <span>p = P(G ≥ Z) = ${chiDistProbabilityText(oneSideGaussP)}</span>

  </div>
</div>`

При \(Z=2\) правый хвост содержит \(p=0.02275\). Ползунок в HTML позволяет сопоставить положение границы с площадью хвоста.

Односторонняя гауссова шкала

Таблица 13.2: Одностороннее соответствие p и Z

\(Z\)	\(p=1-\Phi(Z)\)
\(1\)	\(1.587\cdot10^{-1}\)
\(2\)	\(2.275\cdot10^{-2}\)
\(3\)	\(1.350\cdot10^{-3}\)
\(4\)	\(3.167\cdot10^{-5}\)
\(5\)	\(2.867\cdot10^{-7}\)

При поиске положительного сигнала обычно используют именно такое соглашение. Порог открытия \(5\sigma\) соответствует правому хвосту из таблицы 13.2. Если \(p>1/2\), формула 13.23 даёт отрицательное \(Z\). Такое число не означает свидетельства в пользу положительного сигнала.

13.7.2 Два хвоста

Если нас интересуют отклонения от заданного значения в обе стороны, складываются две хвостовые вероятности:

\[ p=P(|G|\ge Z)=2[1-\Phi(Z)], \qquad Z=\Phi^{-1}(1-p/2),\quad Z\ge0. \tag{13.24}\]

Для \(Z=2\) теперь получается \(p=0.04550\), вдвое больше одностороннего значения. Выбор определяется вопросом анализа. Его нельзя менять после того, как стал известен знак отклонения.

viewof twoSideGaussZ = Inputs.range([0, 6], {value: 2, step: 0.1, label: "Граница |G|: Z"})
twoSideGaussPdf = g => Math.exp(-g*g/2)/Math.sqrt(2*Math.PI)
twoSideGaussP = chiDistSf(twoSideGaussZ**2,1)
twoSideGaussXmax = Math.max(4.5, twoSideGaussZ+1)
twoSideGaussData = Array.from({length:1201},(_,i)=>{
  const g = -twoSideGaussXmax + 2*twoSideGaussXmax*i/1200;
  return {g, density:twoSideGaussPdf(g)};
})
twoSideGaussRight = [{g:twoSideGaussZ,density:twoSideGaussPdf(twoSideGaussZ)},...twoSideGaussData.filter(d=>d.g>twoSideGaussZ)]
twoSideGaussLeft = [...twoSideGaussData.filter(d=>d.g<-twoSideGaussZ),{g:-twoSideGaussZ,density:twoSideGaussPdf(twoSideGaussZ)}]
twoSideGaussPlot = Plot.plot({
  width:760,height:340,marginLeft:58,marginBottom:48,
  style:{fontSize:"14px"},
  x:{label:"g",domain:[-twoSideGaussXmax,twoSideGaussXmax],grid:true},
  y:{label:"φ(g)",domain:[0,0.43],grid:true},
  marks:[
    Plot.areaY(twoSideGaussRight,{x:"g",y:"density",fill:"#c44536",fillOpacity:0.4}),
    Plot.areaY(twoSideGaussLeft,{x:"g",y:"density",fill:"#c44536",fillOpacity:0.4}),
    Plot.line(twoSideGaussData,{x:"g",y:"density",stroke:"#367ab5",strokeWidth:2.5}),
    Plot.ruleX([-twoSideGaussZ,twoSideGaussZ],{stroke:"#c44536",strokeWidth:2})
  ]
})
html`<div class="book-applet-panel">
  ${twoSideGaussPlot}
  <div class="book-applet-summary">
    <span>Z = ${twoSideGaussZ.toFixed(1)}</span>
    <span>p = P(|G| ≥ Z) = ${chiDistProbabilityText(twoSideGaussP)}</span>
    <span>Каждый хвост: ${chiDistProbabilityText(twoSideGaussP/2)}</span>
  </div>
</div>`

За пределами \([-2,2]\) находится вероятность \(0.04550\). Каждый из двух хвостов содержит половину этой вероятности.

Два гауссовых хвоста

Уже знакомый интервал для одного параметра в квадратичном приближении строится из

\[ \Delta\chi^2(\theta) \simeq\frac{(\theta-\widehat\theta)^2}{\sigma_\theta^2}. \tag{13.25}\]

Его границы \(\widehat\theta\pm Z\sigma_\theta\) соответствуют \(\Delta\chi^2=Z^2\). В точной линейной гауссовой задаче покрытие такого интервала равно \(2\Phi(Z)-1\).

Таблица 13.3: Двусторонние гауссовы интервалы одного параметра

Границы интервала	Вероятность двух хвостов	\(\Delta\chi^2\) на границе
\(\widehat\theta\pm\sigma_\theta\)	\(3.173\cdot10^{-1}\)	\(1\)
\(\widehat\theta\pm2\sigma_\theta\)	\(4.550\cdot10^{-2}\)	\(4\)
\(\widehat\theta\pm3\sigma_\theta\)	\(2.700\cdot10^{-3}\)	\(9\)
\(\widehat\theta\pm4\sigma_\theta\)	\(6.334\cdot10^{-5}\)	\(16\)
\(\widehat\theta\pm5\sigma_\theta\)	\(5.733\cdot10^{-7}\)	\(25\)

Это утверждение относится к приращению функции для одного параметра. Подставлять вместо него полное \(\chi^2_{\min}\) нельзя.

13.8 Читаем графики поиска бозона Хиггса

Вернёмся к рисункам из первой главы. По горизонтали стоит пробная масса \(m_H\), по вертикали — локальное p-значение фоновой гипотезы. В каждой точке спрашивают: как часто один фон даёт избыток, настолько же похожий на сигнал с этой массой?

Рисунок 13.1: Локальное p-значение фоновой гипотезы в ATLAS. Сплошная чёрная кривая получена из данных; пунктир и голубая полоса относятся к ожидаемому результату при наличии бозона Хиггса Стандартной модели с пробной массой. Рисунок из статьи ATLAS [2].

На рисунке 13.1 минимум около \(126\) ГэВ достигает локальной значимости \(5.9\sigma\). Вертикальная шкала логарифмическая. Красные линии переводят её в односторонние \(Z\): чем глубже минимум p-значения, тем сильнее несогласие с фоном.

Пунктирная кривая отвечает другому расчёту: какой результат ожидается, если сигнал Стандартной модели действительно есть при каждой проверяемой массе. Она не является подгонкой чёрной кривой.

Рисунок 13.2: Локальные p-значения в предварительном результате CMS с данными при 7 и 8 ТэВ и светимостью до 12.2 фб\(^{-1}\) при 8 ТэВ. Цветные кривые показывают отдельные каналы распада, чёрная — их объединение. Источник: CMS PAS HIG-12-045 [3].

В первой статье об открытии CMS сообщала локальную значимость \(5.0\sigma\) [4]. На рисунке 13.2 используется уже больший набор данных, и минимум стал глубже. Сравнивая результаты, нужно проверять накопленную светимость и версию анализа: число сигма относится к конкретной выборке.

Обе иллюстрации относятся к проверке фоновой гипотезы. Качество согласия модели «сигнал плюс фон» с данными проверяется отдельно. Малое p-значение фона само по себе не сообщает, насколько подробно выбранная модель сигнала описывает все свойства событий.

13.8.1 Локальный результат и поиск по массе

При фиксированной заранее массе получается локальное p-значение. Если массу выбирали по самому глубокому минимуму в диапазоне, в повторяемом опыте нужно каждый раз просматривать весь этот диапазон. Так определяется глобальное p-значение: вероятность найти где-либо в нём столь же сильное отклонение [5].

Соседние массы проверяются по одним и тем же данным, поэтому их результаты коррелированы. Число точек на нарисованной кривой не равно числу независимых попыток.

13.8.2 Откуда возникает правило \(Z=\sqrt{q_0}\)

Покажем его на простой гауссовой модели. Пусть измеряется амплитуда сигнала \(Y\sim\mathcal N(s,\sigma^2)\), где \(\sigma\) известна, а физически \(s\ge0\). Проверяется фон \(s=0\). Подгонка с ограничением даёт \(\widehat s=\max(0,y)\). Улучшение суммы квадратов равно

\[ q_0=\chi^2(0)-\chi^2(\widehat s) = \begin{cases} y^2/\sigma^2,&y>0,\\ 0,&y\le0. \end{cases} \tag{13.26}\]

При верном фоне половина опытов имеет \(Y\le0\) и даёт \(q_0=0\). Положительная часть соответствует половине распределения \(\chi^2_1\). Для полученного \(q_0>0\) имеем

\[ p_0=P_{s=0}(q_0(Y)\ge q_0) =1-\Phi(\sqrt{q_0}), \qquad Z=\sqrt{q_0}. \tag{13.27}\]

Поэтому \(q_0=25\) даёт односторонние \(5\sigma\). Для чистого \(\chi^2_1\) вероятность правее 25 вдвое больше: квадрат учитывает оба знака исходного гауссова отклонения.

В задачах поиска частиц аналогичная формула возникает как асимптотика профильной статистики при одном неотрицательном параметре сигнала и выполнении соответствующих условий регулярности [5]. К произвольному минимуму χ² она не относится.

13.9 Когда распределение нужно получить моделированием

13.9.1 Малые пуассоновские счёты

Рассмотрим пример из лекции: двадцать независимых пуассоновских бинов с заранее заданными ожиданиями \(\mu_i\). Здесь нет подгоняемых параметров и не фиксируется суммарное число событий:

\[ N_i\sim\operatorname{Pois}(\mu_i), \qquad T_{\mathrm P}=\sum_{i=1}^{20}\frac{(N_i-\mu_i)^2}{\mu_i}. \tag{13.28}\]

Это статистика Пирсона. При больших ожиданиях каждое \((N_i-\mu_i)/\sqrt{\mu_i}\) приближается к стандартной гауссовой величине, и сумма квадратов приближённо имеет распределение \(\chi^2_{20}\). Но в примере ожидания малы:

Таблица 13.4: Ожидаемые числа событий в двадцати бинах

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

У статистики 13.28 среднее всё равно равно 20. Для каждого бина \(\mathbb E[(N_i-\mu_i)^2]=\mu_i\). Но совпадения среднего недостаточно, чтобы совпало распределение. Используя центральный четвёртый момент Пуассона \(\mathbb E[(N_i-\mu_i)^4]=\mu_i+3\mu_i^2\), получаем

\[ \operatorname{Var}(T_{\mathrm P}) =\sum_{i=1}^{20}\left(2+\frac{1}{\mu_i}\right) \simeq71.14. \tag{13.29}\]

У \(\chi^2_{20}\) дисперсия равна 40. Уже по ширине видно, что аналитическая замена неточна. Малые \(\mu_i\) дают особенно большой добавочный вклад.

13.9.2 Интерактив: распределение из опытов

Для каждого опыта генерируем двадцать счётов из заданных пуассоновских распределений и вычисляем ту же сумму Пирсона. Полученные значения задают распределение \(T_{\mathrm P}\) при нашей модели.

Ползунок задаёт проверяемое значение \(t_{\mathrm{obs}}\). При \(t_{\mathrm{obs}}=29.8\) независимая проверка на миллионе опытов даёт \(p\simeq0.113\) со статистической ошибкой моделирования около \(0.0003\). Формула для \(\chi^2_{20}\) даёт \(0.07316\). Апплет использует 80 000 опытов, поэтому последняя цифра его результата может отличаться.

viewof cowanObs = Inputs.range([0, 60], {value: 29.8, step: 0.1, label: "Значение T_obs"})

cowanNuModel = [
  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
]

cowanRandomFactory = () => {
  let state = 123456789;
  return () => {
    state = (1664525 * state + 1013904223) >>> 0;
    return state / 4294967296;
  };
}

cowanPoisson = (lambda, rand) => {
  const l = Math.exp(-lambda);
  let k = 0;
  let p = 1;
  do {
    k += 1;
    p *= rand();
  } while (p > l);
  return k - 1;
}

cowanToyStats = {
  const rand = cowanRandomFactory();
  const out = [];
  const bCount = 80000;
  for (let b = 0; b < bCount; b++) {
    let t = 0;
    for (const mu of cowanNuModel) {
      const n = cowanPoisson(mu, rand);
      const d = n - mu;
      t += d * d / mu;
    }
    out.push(t);
  }
  return out;
}

cowanHistogram = {
  const binWidth = 1;
  const bins = Array.from({length: 60}, () => 0);
  for (const t of cowanToyStats) {
    if (t >= 0 && t < 60) bins[Math.floor(t / binWidth)] += 1;
  }
  return bins.map((n, i) => ({
    x1: i * binWidth,
    x2: (i + 1) * binWidth,
    y: n / (cowanToyStats.length * binWidth)
  }));
}

cowanChi2Data = Array.from({length: 1200}, (_, i) => {
  const x = 1e-4 + 60 * i / 1199;
  return {x, y: chiDistPdf(x, 20)};
})

cowanChi2Tail = chiDistSf(cowanObs, 20)

cowanModelP = cowanToyStats.filter(t => t >= cowanObs).length / cowanToyStats.length;
cowanModelError = Math.sqrt(cowanModelP * (1 - cowanModelP) / cowanToyStats.length)
cowanExceedances = cowanToyStats.filter(t => t >= cowanObs).length

cowanPlot = Plot.plot({
  width: 760,
  height: 340,
  marginLeft: 58,
  marginBottom: 48,
  style: {fontSize: "14px"},
  x: {label: "T", domain: [0, 60], grid: true},
  y: {label: "плотность", domain: [0, 0.10], grid: true},
  marks: [
    Plot.rectY(cowanHistogram, {
      x1: "x1",
      x2: "x2",
      y: "y",
      fill: "#8d9ea8",
      fillOpacity: 0.62
    }),
    Plot.line(cowanChi2Data, {
      x: "x",
      y: "y",
      stroke: "#367ab5",
      strokeWidth: 2.5,
      strokeDasharray: "6,4"
    }),
    Plot.ruleX([cowanObs], {stroke: "#c44536", strokeWidth: 2}),
    Plot.text([{x: Math.min(cowanObs + 5, 52), y: 0.075, label: "T_obs"}],
      {x: "x", y: "y", text: "label", fill: "#263640", fontSize: 14})
  ]
})

html`
<div class="book-applet-panel">
  <div class="book-applet-center">${cowanPlot}</div>
  <div class="book-applet-summary">
    <span>80 000 опытов; 20 независимых бинов</span>
    <span>T_obs = ${cowanObs.toFixed(1)}</span>
    <span>Превышений: ${cowanExceedances}</span>
    <span>p из опытов = ${chiDistProbabilityText(cowanModelP)}; ошибка ≈ ${cowanModelError.toFixed(5)}</span>
    <span>p по χ²₂₀ = ${chiDistProbabilityText(cowanChi2Tail)}</span>
  </div>
  <p>Гистограмма — псевдоэксперименты; синий пунктир — χ²₂₀. Нормировка учитывает все опыты, в том числе T ≥ 60.</p>
</div>`

Гистограмма построена по 80 000 опытов для ожиданий из таблицы 13.4. Пунктир — плотность \(\chi^2_{20}\), вертикальная линия — \(t_{\mathrm{obs}}=29.8\). Оба p-значения вычисляются по всему хвосту, включая значения за правой границей рисунка.

Проверка χ²-приближения при малых счётах

Можно выбрать другую статистику, например отношение пуассоновских правдоподобий к насыщенной модели. В насыщенной модели для каждого бина разрешено отдельное среднее, и его оценка равна \(n_i\):

\[ T_{\mathrm L} =2\sum_i\left[ \mu_i-n_i+n_i\ln\frac{n_i}{\mu_i} \right]. \tag{13.30}\]

При \(n_i=0\) логарифмическое слагаемое по пределу равно нулю. Если \(n_i\) близко к \(\mu_i\), разложение логарифма даёт статистику Пирсона в ведущем квадратичном порядке. При малых счётах обе статистики определены, но имеют разные распределения. Само использование правдоподобия не делает распределение 13.30 точно равным χ² при любом числе событий.

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

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

При полностью заданной гипотезе \(H_0\) доля превышений оценивается как

\[ \widehat p=\frac{k}{B}, \qquad \sigma_{\widehat p}\simeq \sqrt{\frac{\widehat p(1-\widehat p)}{B}}, \tag{13.31}\]

где \(B\) — число независимых опытов, \(k\) — число значений \(T\ge t_{\mathrm{obs}}\). Формула ошибки годится, когда превышений и непревышений достаточно много. Если \(k=0\), нулевая оценка ошибки по этой формуле бессодержательна: редкий хвост просто не разрешён данной выборкой. При \(p\simeq2.9\cdot10^{-7}\) даже миллион опытов даёт в среднем меньше одного превышения.

Для численного рангового теста с конечной калибровочной выборкой часто используют \((k+1)/(B+1)\), включая сам проверяемый результат в набор из \(B+1\) равноценных реализаций. Именно так устроен апплет распределения p-значений выше. Этот приём не заменяет проверку модели и не устраняет неопределённость редкого хвоста.

Если \(H_0\) содержит неизвестные параметры, остаётся вопрос, при каких значениях генерировать опыты. Подстановка найденных оценок даёт распространённую приближённую калибровку. Её частотные свойства нужно проверять по допустимым значениям мешающих параметров. Одна генерация при наилучшем фите не гарантирует правильный уровень для всей составной гипотезы.

13.10 Итоги главы

• В линейной гауссовой задаче с известной ковариацией минимум χ² имеет \(N-m\) степеней свободы. Корреляции учитываются ковариацией; подгонка уменьшает размерность пространства остатков.
• p-значение — хвостовая вероятность выбранной статистики при проверяемой гипотезе. Отношение χ²/ν без самого ν не определяет эту вероятность.
• «Число сигма» задаёт гауссово представление p-значения. Один и два хвоста, проверка фона и качество фита требуют разных постановок; их нельзя смешивать.
• При неточной аналитической калибровке распределение статистики получают псевдоэкспериментами с той же процедурой анализа. Число опытов ограничивает точность вычисления редких хвостов.

13.11 Задачи

Задача 1. Две статистики для пуассоновских бинов

В двадцати независимых бинах \(N_i\sim\operatorname{Pois}(\mu_i)\). Первые десять ожиданий равны \((0.158,0.295,0.513,0.833,1.263,1.785,2.355,2.899,3.330,3.568)\), остальные заданы симметрией \(\mu_{21-i}=\mu_i\). Все ожидания известны; параметры по данным не подгоняются.

1. Запишите статистику Пирсона и выведите статистику отношения правдоподобий к насыщенной модели, где каждому бину разрешено своё среднее.
2. Разложите логарифм при \(n_i=\mu_i+\delta_i\), \(|\delta_i|\ll\mu_i\). Сравните ведущие квадратичные члены.
3. Найдите вклады пустого бина в обе статистики. Почему малое число событий само по себе не мешает вычислить эти вклады?
4. Выведите среднее и дисперсию статистики Пирсона, используя \(\mathbb E[(N_i-\mu_i)^4]=\mu_i+3\mu_i^2\).
5. Сгенерируйте 50 000 опытов. Сравните распределения обеих статистик с \(\chi^2_{20}\). Повторите при увеличении всех \(\mu_i\) в 100 раз.

Задача 2. Степени свободы до и после подгонки

Даны \(N=12\) независимых измерений \(Y_i\sim\mathcal N(\mu,\sigma^2)\) при \(\mu=10\) каналах и известной \(\sigma=2\) канала. Проведите 20 000 псевдоэкспериментов.

1. В каждом опыте вычислите среднее \(\overline Y\), сумму квадратов относительно истинного \(\mu\) и минимальную сумму квадратов относительно \(\overline Y\); ошибки в знаменателях равны \(\sigma\).
2. Проверьте точное разложение суммы квадратов на два слагаемых: минимум и \(N(\overline Y-\mu)^2/\sigma^2\).
3. Сравните три распределения с \(\chi^2_{12}\), \(\chi^2_{11}\) и \(\chi^2_1\). Проверьте их средние и дисперсии.
4. Объясните, почему остатки после подгонки не являются двенадцатью независимыми стандартными гауссовыми величинами.
5. Найдите p-значение для \(\chi^2_{\min}=12\). Что получится, если ошибочно оставить двенадцать степеней свободы?

Задача 3. Одинаковый χ² и разная форма остатков

Построим два учебных набора калибровочных данных. Энергии \(E_i=i\) МэВ, \(i=0,\ldots,11\), ошибки амплитуды — по 5 каналов. Матрица \(X\) содержит столбцы \(1\) и \(E_i\).

Для первого набора сгенерируйте двенадцать независимых чисел \(w_i\sim\mathcal N(0,1)\) и удалите их проекцию на прямую: \(\mathbf v=\mathbf w-X(X^{\mathsf T}X)^{-1}X^{\mathsf T}\mathbf w\). Для второго положите \(v_i=(E_i-5.5)^2-\tfrac1{12}\sum_j(E_j-5.5)^2\), используя числовые значения энергий в МэВ. В каждом случае нормируйте остаток: \(\mathbf u=\sqrt{12}\,\mathbf v/\|\mathbf v\|\), и задайте амплитуды \(y_i=2+100E_i+5u_i\) в каналах.

1. Подгоните оба набора прямой с двумя свободными параметрами.
2. Проверьте, что в обоих случаях \(\chi^2_{\min}=12\) и p-значения совпадают.
3. Постройте нормированные остатки по энергии. Какую особенность второго набора скрывает одно число χ²?
4. Объясните, почему эти специально сконструированные наборы нельзя использовать как случайную выборку для проверки распределения \(\chi^2_{\min}\).

Задача 4. Когда p-значения равномерны

Сгенерируйте 20 000 независимых значений \(T\sim\chi^2_{10}\).

1. Вычислите для каждого \(p=P(\chi^2_{10}\ge T)\), постройте гистограмму и найдите долю \(p\le0.05\).
2. Повторите расчёт, заменив каждое \(T\) на \(T/0.8^2\), но сохранив калибровочное распределение \(\chi^2_{10}\). Какой ошибке в задании стандартных отклонений соответствует такая замена?
3. Сгенерируйте отдельную калибровочную выборку из 20 000 значений \(\chi^2_{10}\). Для исходных \(T\) вычислите ранговые p-значения \((k+1)/(B+1)\), где \(k\) — число калибровочных значений не меньше \(T\), \(B=20000\).
4. Повторите опыт для \(N\sim\operatorname{Pois}(2)\), используя точное \(p(n)=P(N\ge n)\). Найдите наименьшее целое \(n\), при котором \(p(n)\le0.05\), и точную вероятность отклонения верной гипотезы.
5. Почему в опыте с Пуассоном доля отклонений не обязана быть ровно 5% даже при неограниченном числе псевдоэкспериментов?

Задача 5. Один хвост, два хвоста и открытие

1. Вычислите \(1-\Phi(Z)\) и \(2[1-\Phi(Z)]\) для \(Z=1,\ldots,5\). Проверьте таблицы главы.
2. Найдите оба перевода в \(Z\) для \(p=0.05\), \(2.7\cdot10^{-3}\) и \(2.87\cdot10^{-7}\).
3. Для \(T=1200\) при \(T\sim\chi^2_{1000}\) сравните точный односторонний перевод p-значения с \((T-1000)/\sqrt{2000}\).
4. В модели \(Y\sim\mathcal N(s,\sigma^2)\) при известной \(\sigma\), \(s\ge0\), выведите \(q_0=\max(0,Y/\sigma)^2\) для проверки \(s=0\). Найдите вероятность \(q_0\ge25\) при фоне. Сравните с \(P(\chi^2_1\ge25)\).
5. Исследователь решил считать отклонение односторонним, но сторону выбирает по знаку результата. Если в каждой выбранной стороне он использует порог 5%, какова полная вероятность ложного отклонения верной гауссовой гипотезы?

Литература

[1] Cowan, G. Statistical Data Analysis. Oxford University Press. 1998.

[2] ATLAS Collaboration. Observation of a New Particle in the Search for the Standard Model Higgs Boson with the ATLAS Detector at the LHC. Physics Letters B, 716, 1–29. 2012. doi:10.1016/j.physletb.2012.08.020.

[3] CMS Collaboration. Combination of Standard Model Higgs Boson Searches and Measurements of the Properties of the New Boson with a Mass near 125 GeV. 2012. https://twiki.cern.ch/twiki/bin/view/CMSPublic/Hig12045TWiki.

[4] CMS Collaboration. Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC. Physics Letters B, 716, 30–61. 2012. doi:10.1016/j.physletb.2012.08.021.

[5] Cowan, G., Cranmer, K., Gross, E., и Vitells, O. Asymptotic Formulae for Likelihood-Based Tests of New Physics. The European Physical Journal C, 71, 1554. 2011. doi:10.1140/epjc/s10052-011-1554-0. https://arxiv.org/abs/1007.1727.