11  Метод наименьших квадратов и линейная подгонка

Аннотация

В этой главе мы применим метод максимального правдоподобия к измерениям с гауссовыми ошибками и получим метод наименьших квадратов. На примере прямой найдём оценки параметров, их ошибки и корреляцию. Затем разберём, с какой точностью известна сама подогнанная кривая и что можно узнать из остатков. Физическим примером послужит калибровка энергетической шкалы детектора.

11.1 Физическая задача: калибровка детектора

Детектор выдаёт амплитуду сигнала в каналах аналого-цифрового преобразователя, сокращённо АЦП. Для физического анализа нужна энергия. Чтобы связать одно с другим, измерим положение нескольких пиков с известными энергиями.

Рассмотрим учебный набор из трёх измерений в таблице 11.1. Энергии считаем известными точно. Ошибки положений пиков независимы, гауссовы и имеют известное стандартное отклонение \(2\) канала. Здесь речь об ошибке измеренного центра пика; ширина самого пика может быть значительно больше.

Таблица 11.1: Учебные данные для калибровки энергетической шкалы

Энергия \(x_i\), МэВ	Положение пика \(y_i\), канал	Ошибка \(\sigma_i\), канал
1	103	2
2	200	2
3	303	2

Предположим, что в этом диапазоне средняя амплитуда линейно зависит от энергии:

\[ f(x)=a+bx. \tag{11.1}\]

Параметр \(a\) задаёт смещение нуля в каналах, а \(b\) — коэффициент преобразования в каналах на МэВ. Подгонка, которую мы сейчас выведем, даёт

\[ \widehat a=(2.0\pm3.1)\ \text{канала}, \qquad \widehat b=(100.0\pm1.4)\ \text{канала}/\text{МэВ}. \tag{11.2}\]

Двух ошибок для описания результата ещё недостаточно. Оценки \(\widehat a\) и \(\widehat b\) сильно антикоррелированы: коэффициент корреляции равен примерно \(-0.93\). Если увеличить наклон, свободный член придётся уменьшить, чтобы прямая по-прежнему проходила вблизи точек.

Откуда взялись эти числа? Почему коррелируют параметры, хотя ошибки трёх измерений независимы? И как эту корреляцию учитывать, когда калибровка уже используется в анализе? Начнём с функции правдоподобия.

11.2 От правдоподобия к сумме квадратов

Пусть имеются \(N\) измерений \((x_i,y_i\pm\sigma_i)\) и модель \(y=f(x;\boldsymbol\theta)\). Во всей главе координаты \(x_i\) фиксированы и известны точно. Случайны результаты \(Y_i\):

\[ Y_i=f(x_i;\boldsymbol\theta)+\varepsilon_i, \qquad \varepsilon_i\sim\mathcal N(0,\sigma_i^2). \tag{11.3}\]

Ошибки \(\varepsilon_i\) независимы, а все \(\sigma_i>0\) известны и не зависят от параметров. Тогда правдоподобие имеет вид [1]

\[ L(\boldsymbol\theta)= \prod_{i=1}^{N}\frac{1}{\sqrt{2\pi}\sigma_i} \exp\!\left[-\frac{[y_i-f(x_i;\boldsymbol\theta)]^2}{2\sigma_i^2}\right]. \tag{11.4}\]

Возьмём логарифм и умножим его на \(-2\):

\[ -2\ln L(\boldsymbol\theta)= \sum_{i=1}^{N}\frac{[y_i-f(x_i;\boldsymbol\theta)]^2}{\sigma_i^2} +2\sum_{i=1}^{N}\ln(\sqrt{2\pi}\sigma_i). \tag{11.5}\]

Вторая сумма постоянна при изменении \(\boldsymbol\theta\). Максимум правдоподобия поэтому совпадает с минимумом функции

\[ \chi^2(\boldsymbol\theta)= \sum_{i=1}^{N}\frac{[y_i-f(x_i;\boldsymbol\theta)]^2}{\sigma_i^2}. \tag{11.6}\]

Это и есть взвешенный метод наименьших квадратов, или МНК.

Разность между измеренным значением и моделью называется остатком. Обозначим её через \(r_i\). Вес точки равен обратной дисперсии:

\[ r_i(\boldsymbol\theta)=y_i-f(x_i;\boldsymbol\theta), \qquad w_i=\frac1{\sigma_i^2}, \qquad \chi^2=\sum_iw_i r_i^2. \tag{11.7}\]

Каждый вклад в \(\chi^2\) безразмерен. Отклонение на два канала для точки с ошибкой один канал даёт вклад \(4\). Для точки с ошибкой четыре канала тот же остаток даёт всего \(1/4\). Так в подгонке учитывается разная точность измерений.

11.2.1 Что мы предположили об ошибках

Связь МНК с максимумом правдоподобия получена для конкретной модели 11.3. Если \(\sigma_i\) зависит от параметров, вторая сумма в формуле 11.5 тоже меняется, и отбрасывать её нельзя. При малом числе событий в бине следует использовать пуассоновское правдоподобие. Подстановка \(\sigma_i=\sqrt{y_i}\) особенно неудачна для пустого бина: его вес тогда вообще не определён.

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

Неопределённости по \(x\) и корреляции между измерениями требуют другого описания данных. В этой главе их пока нет. Физические ограничения на параметры тоже нужно учитывать отдельно: минимум может оказаться на границе, где обычные симметричные ошибки уже не описывают интервал.

11.3 Прямая по точкам

Вернёмся к модели \(f(x)=a+bx\). Нужно минимизировать

\[ \chi^2(a,b)=\sum_iw_i(y_i-a-bx_i)^2. \tag{11.8}\]

Слово «линейная» относится к зависимости от неизвестных параметров. Например, \(a+bx+cx^2\) тоже линейная модель по \(a,b,c\), хотя её график — парабола. Для начала достаточно двух параметров прямой.

11.3.1 Нормальные уравнения

Приравняем производные к нулю:

\[ \begin{aligned} \frac{\partial\chi^2}{\partial a} &=-2\sum_iw_i(y_i-a-bx_i)=0,\\ \frac{\partial\chi^2}{\partial b} &=-2\sum_iw_ix_i(y_i-a-bx_i)=0. \end{aligned} \tag{11.9}\]

Введём пять сумм. Все они вычисляются по данным и известным весам:

\[ \begin{aligned} S&=\sum_iw_i,& S_x&=\sum_iw_ix_i,& S_y&=\sum_iw_iy_i,\\ S_{xx}&=\sum_iw_ix_i^2,& S_{xy}&=\sum_iw_ix_iy_i. \end{aligned} \tag{11.10}\]

Тогда два условия 11.9 образуют систему

\[ \begin{pmatrix}S&S_x\\S_x&S_{xx}\end{pmatrix} \begin{pmatrix}\widehat a\\\widehat b\end{pmatrix} =\begin{pmatrix}S_y\\S_{xy}\end{pmatrix}. \tag{11.11}\]

Это нормальные уравнения МНК. Шляпки у параметров появились потому, что мы уже ищем их оценки. Обозначим матрицу слева через \(A\), вектор параметров через \(\boldsymbol\beta=(a,b)^{\mathsf T}\), а правую часть через \(\mathbf c\). Решение записывается как \(\widehat{\boldsymbol\beta}=A^{-1}\mathbf c\).

Для матрицы размера \(2\times2\) обратную матрицу можно выписать сразу. С определителем \(\Delta=SS_{xx}-S_x^2\) получаем

\[ \widehat a=\frac{S_{xx}S_y-S_xS_{xy}}{\Delta}, \qquad \widehat b=\frac{SS_{xy}-S_xS_y}{\Delta}. \tag{11.12}\]

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

11.3.2 Можно ли определить наклон?

Введём взвешенный центр координат и меру их разброса:

\[ \bar x_w=\frac{S_x}{S}, \qquad Q_x=\sum_iw_i(x_i-\bar x_w)^2 =S_{xx}-\frac{S_x^2}{S} =\frac{\Delta}{S}. \tag{11.13}\]

При положительных весах \(Q_x\) обращается в нуль только тогда, когда все \(x_i\) одинаковы. В этом случае мы измеряем значение прямой в одной точке. Через неё можно провести прямые с любым наклоном, меняя свободный член. Данные позволяют определить одну комбинацию \(a+bx_1\), а двух отдельных параметров не определяют.

Если есть хотя бы два различных \(x_i\), то \(Q_x>0\) и \(\Delta>0\). Матрица \(A\) положительно определена, поэтому найденная стационарная точка является единственным минимумом.

11.4 Ковариация оценок

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

Положим \(\delta a=a-\widehat a\) и \(\delta b=b-\widehat b\). Раскрывая квадраты в формуле 11.8, получаем

\[ \chi^2(a,b)-\chi^2_{\min} =S(\delta a)^2+2S_x\delta a\delta b+S_{xx}(\delta b)^2. \tag{11.14}\]

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

\[ \frac{L(\boldsymbol\beta)}{L(\widehat{\boldsymbol\beta})} =\exp\!\left[-\frac12 (\boldsymbol\beta-\widehat{\boldsymbol\beta})^{\mathsf T} A(\boldsymbol\beta-\widehat{\boldsymbol\beta})\right]. \tag{11.15}\]

Матрица \(A\) играет роль обратной ковариации. Проверим это непосредственно по распределению оценок, чтобы не приписывать самому правдоподобию смысл распределения параметров.

11.4.1 Дисперсии и ковариация из случайных ошибок

Вектор \(\mathbf c\) из нормальных уравнений содержит две случайные суммы: \(S_y=\sum_iw_iY_i\) и \(S_{xy}=\sum_iw_ix_iY_i\). Подстановка \(Y_i=a+bx_i+\varepsilon_i\) даёт

\[ \begin{aligned} \mathbb E[\mathbf c]&=A\boldsymbol\beta,\\ \operatorname{Cov}(\mathbf c)&= \begin{pmatrix} \sum_iw_i^2\sigma_i^2&\sum_iw_i^2x_i\sigma_i^2\\ \sum_iw_i^2x_i\sigma_i^2&\sum_iw_i^2x_i^2\sigma_i^2 \end{pmatrix} =A. \end{aligned} \tag{11.16}\]

Использовано равенство \(w_i\sigma_i^2=1\) и отсутствие ковариаций между ошибками разных точек. Поскольку \(\widehat{\boldsymbol\beta}=A^{-1}\mathbf c\),

\[ \mathbb E[\widehat{\boldsymbol\beta}]=\boldsymbol\beta, \qquad V_\beta=\operatorname{Cov}(\widehat{\boldsymbol\beta}) =A^{-1}A(A^{-1})^{\mathsf T}=A^{-1}. \tag{11.17}\]

Обращение матрицы даёт

\[ V_\beta=\frac1\Delta \begin{pmatrix}S_{xx}&-S_x\\-S_x&S\end{pmatrix}. \tag{11.18}\]

Итак, оценки несмещённые. Их дисперсии и ковариация — три элемента матрицы 11.18:

\[ \operatorname{Var}(\widehat a)=\frac{S_{xx}}\Delta, \qquad \operatorname{Var}(\widehat b)=\frac{S}\Delta, \qquad \operatorname{cov}(\widehat a,\widehat b)=-\frac{S_x}\Delta. \tag{11.19}\]

В этом выводе понадобились только средние и ковариации ошибок. Гауссова модель дополнительно гарантирует гауссово распределение вектора оценок: он является линейной функцией гауссовых данных.

Обратите внимание: в матрице \(V_\beta\) вообще нет значений \(y_i\). При фиксированных \(x_i\) и известных \(\sigma_i\) случайные значения \(y_i\) меняют положение минимума, но его кривизна остаётся той же. Это свойство именно линейной модели с заданными ошибками.

11.5 Геометрия ошибок

Через \(\bar x_w\) и \(Q_x\) формулы 11.19 переписываются в виде

\[ \sigma_b^2=\frac1{Q_x}, \qquad \sigma_a^2=\frac1S+\frac{\bar x_w^2}{Q_x}, \qquad \operatorname{cov}(\widehat a,\widehat b)=-\frac{\bar x_w}{Q_x}. \tag{11.20}\]

Здесь \(\sigma_a\) и \(\sigma_b\) обозначают стандартные отклонения оценок. Чем шире разнесены точки по \(x\) при тех же ошибках, тем больше \(Q_x\) и тем точнее определяется наклон. Для калибровки это означает: пики в узком энергетическом диапазоне плохо определяют коэффициент преобразования.

Со свободным членом ситуация другая. Параметр \(a\) — значение прямой при \(x=0\). Если калибровочные точки находятся далеко от нуля, приходится экстраполировать к нему прямую. Неопределённость наклона увеличивает неопределённость такого пересечения.

11.5.1 Через какую точку проходит прямая

Зафиксируем произвольный наклон \(b\) и найдём лучшее значение \(a\). Первое нормальное уравнение даёт

\[ \widehat{\widehat a}(b)=\bar y_w-b\bar x_w, \qquad \bar y_w=\frac{S_y}{S}. \tag{11.21}\]

Двойная шляпка имеет тот же смысл, что и в прошлой главе: свободный член подобран при фиксированном наклоне. Подставляя \(x=\bar x_w\), получаем \(f(\bar x_w)=\bar y_w\). Все такие прямые проходят через взвешенный центр данных \((\bar x_w,\bar y_w)\).

Теперь понятен знак ковариации. Если \(\bar x_w>0\), увеличение наклона компенсируется уменьшением свободного члена. При \(\bar x_w<0\) оба параметра меняются в одну сторону. Коэффициент корреляции равен

\[ \rho_{ab}= \frac{\operatorname{cov}(\widehat a,\widehat b)}{\sigma_a\sigma_b} =-\frac{S_x}{\sqrt{SS_{xx}}}. \tag{11.22}\]

11.5.2 Интерактив: центр точек и наклон

В апплете можно менять взвешенный центр \(\bar x_w\) и среднеквадратичный разброс координат \(s_x=\sqrt{Q_x/S}\). Суммарный вес \(S=18\) остаётся фиксированным. Координаты здесь безразмерные. Учебные точки построены так, чтобы их лучшая прямая была \(y=1+0.55x\).

Сплошная оранжевая линия показывает эту прямую. Две синие пунктирные линии соответствуют наклонам \(\widehat b\pm\sigma_b\) с заново подобранным свободным членом. На втором графике изображён контур \(\Delta\chi^2=2.30\) и направление профилирования \(\delta a=-\bar x_w\delta b\).

viewof geomXbar = Inputs.range([-3, 3], {value: 1.8, step: 0.1, label: "Центр x̄_w"})
viewof geomSpread = Inputs.range([0.45, 2.5], {value: 1.0, step: 0.05, label: "Разброс s_x"})

geomS = 18
geomQx = geomS * geomSpread * geomSpread
geomSigmaB = 1 / Math.sqrt(geomQx)
geomSigmaA = Math.sqrt(1 / geomS + geomXbar * geomXbar / geomQx)
geomCovAB = -geomXbar / geomQx
geomCorrAB = geomCovAB / (geomSigmaA * geomSigmaB)

geomA0 = 1.0
geomB0 = 0.55
geomYbar = geomA0 + geomB0 * geomXbar
geomSlopeSteps = [-1, 0, 1]
geomLineData = geomSlopeSteps.flatMap(k => {
  const b = geomB0 + k * geomSigmaB;
  const a = geomYbar - b * geomXbar;
  return Array.from({length: 121}, (_, i) => {
    const x = geomXLimits[0] + (geomXLimits[1] - geomXLimits[0]) * i / 120;
    return {x, y: a + b * x, k};
  });
})
geomFanLineData = geomLineData.filter(d => d.k !== 0)
geomBestLineData = geomLineData.filter(d => d.k === 0)

geomOffsets = {
  const raw = [-1.55, -1.05, -0.55, -0.12, 0.35, 0.82, 1.28, 1.62];
  const mean = raw.reduce((a, b) => a + b, 0) / raw.length;
  const rms = Math.sqrt(raw.reduce((s, u) => s + (u - mean) ** 2, 0) / raw.length);
  return raw.map(u => (u - mean) / rms);
}
geomPointSigma = Math.sqrt(geomOffsets.length / geomS)
geomResidualPattern = {
  const raw = [0.09, -0.04, 0.06, -0.08, 0.02, 0.07, -0.05, -0.07];
  const mean = raw.reduce((a, b) => a + b, 0) / raw.length;
  const slope = raw.reduce((s, r, i) => s + r * geomOffsets[i], 0) / raw.length;
  return raw.map((r, i) => r - mean - slope * geomOffsets[i]);
}
geomPoints = geomOffsets.map((u, i) => {
  const x = geomXbar + geomSpread * u;
  return {x, y: geomA0 + geomB0 * x + geomResidualPattern[i]};
})

geomCenterPoint = [{x: geomXbar, y: geomYbar}]
geomXLimits = [Math.min(0, ...geomPoints.map(d => d.x)) - 0.5, Math.max(0, ...geomPoints.map(d => d.x)) + 0.5]

geomDataPlot = Plot.plot({
  width: 760,
  height: 320,
  style: {fontSize: "14px"},
  marginLeft: 58,
  marginBottom: 46,
  x: {label: "x", domain: geomXLimits, grid: true},
  y: {label: "y", grid: true},
  marks: [
    Plot.ruleX([0], {stroke: "#aaa", strokeDasharray: "6,5"}),
    Plot.ruleX([geomXbar], {stroke: "#ffb347", strokeWidth: 3}),
    Plot.line(geomFanLineData, {x: "x", y: "y", z: "k", stroke: "#367ab5", strokeWidth: 2.4, strokeDasharray: "7,5"}),
    Plot.line(geomBestLineData, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 4}),
    Plot.ruleX(geomPoints, {x: "x", y1: d => d.y - geomPointSigma, y2: d => d.y + geomPointSigma, stroke: "#aaa", strokeOpacity: 0.65}),
    Plot.dot(geomPoints, {x: "x", y: "y", r: 6, fill: "#4ea1ff", fillOpacity: 0.82}),
    Plot.dot(geomCenterPoint, {x: "x", y: "y", r: 8, fill: "#ffb347", stroke: "#111", strokeWidth: 1})
  ]
})

geomSa2 = geomSigmaA * geomSigmaA
geomSb2 = geomSigmaB * geomSigmaB
geomL11 = geomSigmaA
geomL21 = geomCovAB / geomL11
geomL22 = Math.sqrt(Math.max(geomSb2 - geomL21 * geomL21, 0))
geomLevel = Math.sqrt(2.30)
geomEllipse = Array.from({length: 241}, (_, i) => {
  const t = 2 * Math.PI * i / 240;
  const z1 = geomLevel * Math.cos(t);
  const z2 = geomLevel * Math.sin(t);
  return {
    da: geomL11 * z1,
    db: geomL21 * z1 + geomL22 * z2
  };
})
geomDaMax = Math.max(0.55, 1.18 * Math.max(...geomEllipse.map(d => Math.abs(d.da))))
geomDbMax = Math.max(0.42, 1.18 * Math.max(...geomEllipse.map(d => Math.abs(d.db))))
geomCompLine = [
  {db: -geomDbMax, da: geomXbar * geomDbMax},
  {db: geomDbMax, da: -geomXbar * geomDbMax}
]

geomParamPlot = Plot.plot({
  width: 760,
  height: 320,
  style: {fontSize: "14px"},
  marginLeft: 66,
  marginBottom: 46,
  x: {label: "δa", domain: [-geomDaMax, geomDaMax], grid: true},
  y: {label: "δb", domain: [-geomDbMax, geomDbMax], grid: true},
  marks: [
    Plot.ruleX([0], {stroke: "#888"}),
    Plot.ruleY([0], {stroke: "#888"}),
    Plot.line(geomEllipse, {x: "da", y: "db", stroke: "#ffb347", strokeWidth: 4}),
    Plot.line(geomCompLine, {x: "da", y: "db", stroke: "#367ab5", strokeWidth: 3, strokeDasharray: "7,5", clip: true}),
    Plot.dot([{da: 0, db: 0}], {x: "da", y: "db", r: 6, fill: "#263640"})
  ]
})

html`
<div class="book-applet-panel">
  <div class="book-applet-column">
    ${geomDataPlot}
    ${geomParamPlot}
  </div>
  <div class="book-applet-summary">
    <span>S = ${geomS}; Q_x = ${geomQx.toFixed(2)}</span>
    <span>σ_y фиксирована: ${geomPointSigma.toFixed(2)}</span>
    <span>σ_b = ${geomSigmaB.toFixed(3)}</span>
    <span>σ_a = ${geomSigmaA.toFixed(3)}</span>
    <span>ρ(â,b̂) = ${geomCorrAB.toFixed(2)}</span>
    <span>Контур: Δχ² = 2.30</span>
  </div>
</div>
`

При \(\bar x_w=1.8\) и \(s_x=1\) корреляция оценок равна \(-0.87\). Перенос центра к нулю устраняет наклон эллипса. Увеличение разброса точек уменьшает ошибку наклона прямой.

Геометрия корреляции параметров

11.5.3 Центрированная параметризация

Можно сразу измерять \(x\) от центра данных:

\[ x'=x-\bar x_w, \qquad f(x)=a'+bx', \qquad a'=a+b\bar x_w. \tag{11.23}\]

Теперь \(\sum_iw_ix_i'=0\), поэтому

\[ \widehat a'=\bar y_w, \qquad \operatorname{Var}(\widehat a')=\frac1S, \qquad \operatorname{cov}(\widehat a',\widehat b)=0. \tag{11.24}\]

Наклон и его ошибка сохранились. Параметр \(a'\) теперь описывает значение прямой в центре данных. Сама калибровка не изменилась, но её параметры стали некоррелированными. В гауссовой модели оценки \(\widehat a'\) и \(\widehat b\) также независимы.

11.5.4 Нормальные координаты и профиль

В новых координатах квадратичная форма диагональна:

\[ \Delta\chi^2=S(\delta a')^2+Q_x(\delta b)^2 =z_a^2+z_b^2, \qquad z_a=\sqrt S\,\delta a',\quad z_b=\sqrt{Q_x}\,\delta b. \tag{11.25}\]

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

При профилировании по \(a'\) первое слагаемое обращается в нуль. Остаётся

\[ \chi^2_{\mathrm p}(b)-\chi^2_{\min} =Q_x(b-\widehat b)^2 =\frac{(b-\widehat b)^2}{\sigma_b^2}. \tag{11.26}\]

Поэтому пределы \(\widehat b\pm\sigma_b\) соответствуют росту профильного \(\chi^2\) на единицу. Это одномерный интервал с покрытием около \(68.3\%\) в нашей гауссовой модели. Для совместной области двух параметров с тем же покрытием нужен порог \(\Delta\chi^2=2.30\). Именно он показан в апплете.

11.6 Возвращаемся к калибровке

Для данных из таблицы 11.1 все веса равны \(1/4\) в обратных квадратах канала. Получаем

\[ \bar x_w=2\ \text{МэВ}, \qquad \bar y_w=202\ \text{канала}, \qquad Q_x=0.5\ \text{МэВ}^2/\text{канал}^2. \tag{11.27}\]

Из формул 11.12 и 11.20 следуют \(\widehat b=100\) каналов/МэВ и \(\widehat a=2\) канала, \(\sigma_b=\sqrt2\) канала/МэВ и \(\sigma_a=\sqrt{28/3}\) канала. Ковариация равна \(-4\) канал\(^2\)/МэВ. Это и даёт результат 11.2 с корреляцией \(-\sqrt{6/7}\simeq-0.93\).

В центрированной записи та же прямая выглядит проще:

\[ \widehat f(x)=202\ \text{канала} +\left(100\ \frac{\text{каналов}}{\text{МэВ}}\right)(x-2\ \text{МэВ}). \tag{11.28}\]

Значение в центре равно \(202\pm1.15\) канала и не коррелирует с наклоном. Ошибка \(1.15\) заметно меньше ошибки свободного члена \(3.1\): центр данных находится внутри диапазона измерений, а нулевая энергия — за его пределами.

11.6.1 С какой точностью известна прямая?

В произвольной фиксированной точке \(x\) оценка среднего сигнала равна \(\widehat f(x)=\widehat a+\widehat bx\). Закон распространения ошибок даёт

\[ \begin{aligned} \operatorname{Var}[\widehat f(x)] &=\sigma_a^2+2x\operatorname{cov}(\widehat a,\widehat b)+x^2\sigma_b^2\\ &=\frac1S+\frac{(x-\bar x_w)^2}{Q_x}. \end{aligned} \tag{11.29}\]

Это точное равенство: \(\widehat f\) линейна по оценкам. Минимальная ошибка достигается в центре данных. По мере удаления от центра полоса вокруг подогнанной прямой расширяется.

Если забыть ковариацию, уже при \(x=2\) МэВ получится дисперсия \(28/3+8=52/3\) канал\(^2\) вместо \(4/3\) канал\(^2\). Ошибка вырастет с \(1.15\) до \(4.16\) канала. Сильная корреляция здесь имеет вполне измеримое последствие.

11.6.2 Средний сигнал и новое измерение

Теперь проведём ещё одно независимое измерение при фиксированной энергии \(x_*\), с гауссовой ошибкой известной ширины \(\sigma_*\):

\[ Y_*=a+bx_*+\varepsilon_*, \qquad \operatorname{Var}[Y_*-\widehat f(x_*)] =\sigma_*^2+\frac1S+\frac{(x_*-\bar x_w)^2}{Q_x}. \tag{11.30}\]

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

В центре нашей калибровки средний сигнал известен с ошибкой \(1.15\) канала. Если новое измерение имеет \(\sigma_*=2\) канала, стандартное отклонение разности \(Y_*-\widehat f(x_*)\) равно \(\sqrt{4+4/3}=2.31\) канала. Обе полосы показаны на рисунке 11.1.

Рисунок 11.1: Учебная калибровка из таблицы: точки с ошибками и подогнанная прямая. На нижнем графике из каждой величины вычтена подогнанная прямая, чтобы различить полосы. Узкая полоса соответствует одной стандартной ошибке среднего сигнала; широкая учитывает также ошибку нового независимого измерения \(\sigma_*=2\) канала.

Границы обеих полос имеют точечное покрытие около \(68.3\%\) при каждом заранее фиксированном \(x\). Утверждение о всей кривой сразу требует одновременной доверительной полосы с другим порогом. Кроме того, за пределами калибровочного диапазона эти формулы учитывают лишь ошибки параметров прямой. Возможную нелинейность реального отклика они не описывают.

11.7 Что остаётся после подгонки

Остатки при найденных параметрах равны

\[ r_i=y_i-\widehat a-\widehat bx_i, \qquad u_i=\frac{r_i}{\sigma_i}. \tag{11.31}\]

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

\[ \sum_iw_i r_i=0, \qquad \sum_iw_ix_i r_i=0. \tag{11.32}\]

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

11.7.1 Матрица влияния

Разделим данные и модель в каждой строке на \(\sigma_i\):

\[ z_i=\frac{y_i}{\sigma_i}, \qquad B_{i1}=\frac1{\sigma_i}, \qquad B_{i2}=\frac{x_i}{\sigma_i}. \tag{11.33}\]

Матрица \(B\) имеет \(N\) строк и два столбца. Ошибки вектора \(\mathbf z\) независимы и имеют единичную дисперсию. Нормальные уравнения теперь записываются как \(B^{\mathsf T}B\widehat{\boldsymbol\beta}=B^{\mathsf T}\mathbf z\). Значения подогнанной модели в этих координатах равны

\[ \widehat{\mathbf z}=H\mathbf z, \qquad H=B(B^{\mathsf T}B)^{-1}B^{\mathsf T}. \tag{11.34}\]

Матрица \(H\) называется матрицей влияния. Она симметрична и удовлетворяет \(H^2=H\): если ещё раз подогнать уже лежащие на прямой значения, они не изменятся. Вектор нормированных остатков равен \(\mathbf u=(I-H)\mathbf z\), где \(I\) — единичная матрица. Поэтому

\[ \operatorname{Cov}(\mathbf u) =(I-H)(I-H)^{\mathsf T}=I-H. \tag{11.35}\]

В частности,

\[ \operatorname{Var}(u_i)=1-h_{ii}, \qquad h_{ii}=w_i\left[\frac1S+\frac{(x_i-\bar x_w)^2}{Q_x}\right]. \tag{11.36}\]

Точка с большим \(h_{ii}\) сильнее подтягивает к себе прямую, поэтому её остаток имеет меньшую дисперсию. При прочих равных большое влияние имеют точки далеко от центра \(x\). Малый остаток такой точки сам по себе ещё не означает, что с измерением всё в порядке.

Для известной дисперсии ошибки и \(h_{ii}<1\) величина \(u_i/\sqrt{1-h_{ii}}\) имеет единичную дисперсию, а в нашей гауссовой модели и стандартное нормальное распределение. Эта нормировка не устраняет корреляций между разными остатками.

11.7.2 Почему остаётся \(N-2\) степеней свободы

Матрица \(H\) проектирует данные на двумерное пространство моделей, заданное двумя столбцами \(B\). У неё два собственных значения \(1\), остальные равны \(0\), и \(\operatorname{tr}H=2\). Для остатков остаётся \(N-2\) независимых направлений. При верной линейной модели и известных гауссовых ошибках

\[ \chi^2_{\min}=\mathbf u^{\mathsf T}\mathbf u \sim\chi^2_{N-2}, \qquad \mathbb E[\chi^2_{\min}]=N-2. \tag{11.37}\]

Распределение справа обозначено тем же символом \(\chi^2\), что и функция, которую мы минимизировали. Здесь \(N>2\). Для нашей калибровки остатки равны \(1,-2,1\) канала, поэтому \(\chi^2_{\min}=1.5\) при одной степени свободы. Из одного такого числа нельзя заключить, что модель детектора верна. Формальную проверку согласия разберём в главе 13.

Если ошибки всех точек имеют одинаковую, но неизвестную дисперсию \(\sigma^2\), её можно оценить по остаткам:

\[ s^2=\frac1{N-2}\sum_i r_i^2, \qquad \mathbb E[s^2]=\sigma^2. \tag{11.38}\]

Для интервалов тогда требуется учитывать случайность \(s\), что приводит к распределению Стьюдента. В нашей основной задаче все \(\sigma_i\) заданы заранее. Перемасштабировать их до \(\chi^2_{\min}/(N-2)=1\) оснований нет.

11.8 Остатки и выбросы

Сгенерируем \(16\) точек около прямой \(y=1+0.8x\) с независимыми гауссовыми ошибками. Обе координаты здесь безразмерные. В апплете можно изменить масштаб ошибок и добавить сдвиг \(d\) к одной отмеченной точке. При изменении сдвига остальные точки сохраняются. Кнопка «Новые данные» меняет всю случайную выборку.

viewof lsSigma = Inputs.range([0.15, 1.5], {value: 0.55, step: 0.05, label: "Масштаб ошибок σ"})
viewof lsOutlier = Inputs.range([-4, 6], {value: 0, step: 0.25, label: "Сдвиг точки d"})
viewof regenerateLs = Inputs.button("Новые данные", {value: 0, reduce: n => n + 1})

lsRandom = seed => {
  let a = seed >>> 0;
  return () => {
    a = (a + 0x6D2B79F5) >>> 0;
    let t = a;
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}

lsNoise = {
  const random = lsRandom(20260911 + 104729 * regenerateLs);
  return Array.from({length: 16}, () => {
    const u1 = Math.max(random(), 1e-12);
    return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * random());
  });
}

lsData = {
  const data = [];
  for (let i = 0; i < 16; ++i) {
    const x = -4.5 + 9 * i / 15;
    const sigma = lsSigma * (0.72 + 0.56 * (i % 4) / 3);
    const y = 1 + 0.8 * x + sigma * lsNoise[i] + (i === 12 ? lsOutlier : 0);
    data.push({x, y, sigma, outlier: i === 12});
  }
  return data;
}

weightedLineFit = data => {
  let S = 0, Sx = 0, Sy = 0, Sxx = 0, Sxy = 0;
  for (const d of data) {
    const w = 1 / (d.sigma * d.sigma);
    S += w;
    Sx += w * d.x;
    Sy += w * d.y;
    Sxx += w * d.x * d.x;
    Sxy += w * d.x * d.y;
  }
  const delta = S * Sxx - Sx * Sx;
  const a = (Sxx * Sy - Sx * Sxy) / delta;
  const b = (S * Sxy - Sx * Sy) / delta;
  const chi2 = data.reduce((sum, d) => sum + ((d.y - a - b * d.x) / d.sigma) ** 2, 0);
  return {a, b, chi2, sigmaA: Math.sqrt(Sxx / delta), sigmaB: Math.sqrt(S / delta), covAB: -Sx / delta};
}

lsFit = weightedLineFit(lsData)
lsFitLine = Array.from({length: 201}, (_, i) => {
  const x = -5 + 10 * i / 200;
  return {x, y: lsFit.a + lsFit.b * x};
})
lsTrueLine = Array.from({length: 201}, (_, i) => {
  const x = -5 + 10 * i / 200;
  return {x, y: 1 + 0.8 * x};
})
lsResiduals = lsData.map(d => ({
  x: d.x,
  pull: (d.y - lsFit.a - lsFit.b * d.x) / d.sigma,
  outlier: d.outlier
}))

lsMainPlot = Plot.plot({
  width: 760,
  height: 340,
  style: {fontSize: "14px"},
  marginLeft: 65,
  x: {label: "x", domain: [-5, 5], grid: true},
  y: {label: "y", grid: true},
  marks: [
    Plot.ruleX(lsData, {x: "x", y1: d => d.y - d.sigma, y2: d => d.y + d.sigma, stroke: "#aaa"}),
    Plot.dot(lsData.filter(d => !d.outlier), {x: "x", y: "y", r: 6, fill: "#4ea1ff"}),
    Plot.dot(lsData.filter(d => d.outlier), {x: "x", y: "y", r: 6, symbol: "square", fill: "#c44536"}),
    Plot.line(lsTrueLine, {x: "x", y: "y", stroke: "#999", strokeDasharray: "7,5", strokeWidth: 3}),
    Plot.line(lsFitLine, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 4})
  ]
})

lsResidualPlot = Plot.plot({
  width: 760,
  height: 260,
  style: {fontSize: "14px"},
  marginLeft: 65,
  x: {label: "x", domain: [-5, 5], grid: true},
  y: {label: "u = остаток / σᵢ", domain: [-lsResidualLimit, lsResidualLimit], grid: true},
  marks: [
    Plot.ruleY([0], {stroke: "#263640"}),
    Plot.ruleY([-2, 2], {stroke: "#999", strokeDasharray: "6,5"}),
    Plot.dot(lsResiduals.filter(d => !d.outlier), {x: "x", y: "pull", r: 6, fill: "#4ea1ff"}),
    Plot.dot(lsResiduals.filter(d => d.outlier), {x: "x", y: "pull", r: 6, symbol: "square", fill: "#c44536"})
  ]
})

lsResidualLimit = Math.max(3, 1.15 * Math.max(...lsResiduals.map(d => Math.abs(d.pull))))

html`
<div class="book-applet-panel">
  <div class="book-applet-column">
    ${lsMainPlot}
    ${lsResidualPlot}
  </div>
  <div class="book-applet-summary">
    <span>â = ${lsFit.a.toFixed(3)} ± ${lsFit.sigmaA.toFixed(3)}</span>
    <span>b̂ = ${lsFit.b.toFixed(3)} ± ${lsFit.sigmaB.toFixed(3)}</span>
    <span>χ²_min = ${lsFit.chi2.toFixed(1)}</span>
    <span>N − 2 = 14</span>
    <span>Пунктир: истинная прямая; сплошная: подгонка</span>
  </div>
</div>
`

Подгонка \(16\) псевдоизмерений при масштабе ошибок \(0.55\) и без дополнительного сдвига. Пунктиром показана истинная прямая, сплошной линией — результат подгонки. Нижний график содержит остатки, делённые на ошибки точек.

Прямая, остатки и выброс

На рисунке 11.2 к той же выборке добавлен сдвиг \(d=4\) в одной точке. Её заявленная ошибка не изменилась. Подгонка учитывает эту точку с прежним весом, и прямая смещается.

Рисунок 11.2: Те же псевдоизмерения после сдвига отмеченной точки на \(d=4\). Пунктирная прямая задаёт исходную модель, сплошная получена из подгонки. Сдвиг одной точки изменил и её собственный остаток, и остатки остальных измерений.

Этот эффект можно вычислить без повторной минимизации. Если заменить \(y_j\) на \(y_j+d\) при неизменных координатах и весах, то

\[ \delta\widehat b=\frac{w_j(x_j-\bar x_w)}{Q_x}\,d, \qquad \delta\widehat a=\frac{w_j}{S}\,d-\bar x_w\delta\widehat b. \tag{11.39}\]

Смещение наклона растёт линейно с \(d\) и с расстоянием точки от центра. Обычный МНК не ограничивает влияние сколь угодно большого отклонения.

Что делать, если такая точка встретилась в эксперименте? Нужно проверить причину: сбой электроники, неверную ошибку измерения, примесь другого процесса, насыщение детектора или недостаточность самой модели. Удаление точки требует физического основания и явно описанного правила отбора. Если большие отклонения составляют часть распределения ошибок, это распределение следует включить в правдоподобие.

На графике остатков полезно смотреть также на согласованное поведение многих точек. Изгиб может указывать на нелинейность отклика. Рост разброса с \(x\) может означать, что зависимость ошибок от энергии описана неверно. Такой график помогает поставить вопрос к модели; сам по себе он ещё не выбирает причину отклонения.

11.9 Итоги главы

• Для независимых гауссовых измерений с известными, не зависящими от параметров дисперсиями МНК совпадает с методом максимального правдоподобия.
• Подгонка прямой сводится к двум нормальным уравнениям. Ковариация оценок определяется координатами точек и их ошибками.
• Центрирование координаты устраняет ковариацию наклона и значения прямой в центре данных. При вычислении ошибки кривой в исходных параметрах ковариационный член необходимо сохранять.
• После подгонки остатки связаны между собой и имеют уменьшенные дисперсии. Их структура и влияние отдельных точек позволяют проверять описание измерений.

11.10 Задачи

Задача 1. Среднее как результат МНК

Имеются независимые измерения одной величины: \[ Y_i\sim\mathrm N(\mu,\sigma_i^2). \]

Все \(\sigma_i>0\) известны и не зависят от \(\mu\).

1. Запишите функцию \(\chi^2(\mu)\).
2. Найдите её минимум.
3. Покажите, что полученная оценка является взвешенным средним.
4. Найдите дисперсию этой оценки.
5. Рассмотрите частный случай одинаковых \(\sigma_i\).

Задача 2. Аналитическая подгонка прямой

Для данных \((x_i,y_i,\sigma_i)\) рассматривается модель \[ y=a+bx. \]

Координаты \(x_i\) известны точно, ошибки по \(y\) независимы и гауссовы, их стандартные отклонения \(\sigma_i>0\) известны. Среди \(x_i\) есть хотя бы два различных значения.

1. Выведите нормальные уравнения.
2. Получите аналитические формулы для \(\widehat a\) и \(\widehat b\).
3. Выведите ковариационную матрицу параметров.
4. Покажите, что центрирование \(x_i\) около взвешенного среднего устраняет ковариацию параметров.
5. Проверьте формулы численно на искусственной выборке.

Задача 3. Ошибка подогнанной кривой

Три независимых гауссовых измерения при точно известных энергиях \(x=(1,2,3)\) МэВ дали положения пиков \(y=(103,200,303)\) канала. Ошибка каждого положения равна \(2\) каналам. Для модели \(f(x)=a+bx\) найдите:

1. дисперсию оценки среднего \(\widehat y(x)=\widehat a+\widehat bx\);
2. точку, в которой эта дисперсия минимальна;
3. дисперсию разности \(Y_*-\widehat y(x_*)\) для нового независимого измерения при фиксированном \(x_*\) с известной ошибкой \(\sigma_*\);
4. различие между полосой среднего сигнала и полосой для нового измерения.

Постройте обе полосы для этих данных при \(\sigma_*=2\) канала. Объясните, почему точечное покрытие не означает покрытия всей кривой сразу.

Задача 4. Остатки и диагональные элементы матрицы влияния

Пусть \(N>2\), \(y_i=a+bx_i+\varepsilon_i\), а \(\varepsilon_i\) независимы, имеют нулевые средние и одну известную дисперсию \(\sigma^2\). Координаты \(x_i\) фиксированы и не все одинаковы. В этом случае \(\widehat{\mathbf y}=H\mathbf y\), где \(H=X(X^{\mathsf T}X)^{-1}X^{\mathsf T}\), \(X_{i1}=1\), \(X_{i2}=x_i\). Обозначьте остатки через \(\mathbf r=\mathbf y-\widehat{\mathbf y}\), а единичную матрицу через \(I\).

1. покажите, что \(H^2=H\);
2. покажите, что \(\operatorname{tr}H=2\);
3. выведите ковариацию остатков \[ \operatorname{Cov}(\mathbf r)=\sigma^2(I-H); \]
4. найдите диагональные элементы \(h_{ii}\), характеризующие влияние точек при подгонке прямой;
5. добавьте одну далёкую по \(x\) точку и исследуйте её влияние.

Задача 5. Влияние одной точки

Сгенерируйте \(16\) независимых гауссовых измерений около прямой \(y=1+0.8x\) при \(x_i=-4.5+9(i-1)/15\), \(i=1,\ldots,16\), и известных одинаковых ошибках \(\sigma_i=0.55\). Зафиксируйте начальное состояние генератора.

1. Выполните подгонку обычным методом наименьших квадратов.
2. Добавьте к одной точке сдвиг \(d\), сохранив её заявленную ошибку. Выведите изменения \(\widehat a\) и \(\widehat b\) из нормальных уравнений.
3. Меняйте \(d\) от \(-4\) до \(6\) и сравните формулы с повторной подгонкой.
4. Повторите опыт для точки вблизи центра и для крайней точки. Сравните влияние на наклон и остатки.
5. Проверьте, меняется ли вычисленная ковариация параметров при этих сдвигах. Объясните, почему неизменная ошибка параметра не гарантирует отсутствия смещения из-за испорченного измерения.

Литература

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