7  Двумерное гауссово распределение

Аннотация

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

7.1 Физическая задача: два показания калориметра

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

\[ \boldsymbol\mu= \begin{pmatrix} 60\\ 40 \end{pmatrix}\ \text{ГэВ}. \tag{7.1}\]

Дисперсии измеренных энергий и их ковариация собраны в матрицу

\[ V= \begin{pmatrix} 1.44 & 0.36\\ 0.36 & 1.00 \end{pmatrix}\ \text{ГэВ}^2. \tag{7.2}\]

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

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

\[ \mathbf x= \begin{pmatrix} 61.2\\ 41.0 \end{pmatrix}\ \text{ГэВ}. \tag{7.3}\]

Каждая энергия отличается от среднего ровно на своё стандартное отклонение. Этого ещё недостаточно, чтобы оценить совместное отклонение. Для него нужна квадратичная форма

\[ D^2=(\mathbf x-\boldsymbol\mu)^{\mathsf T} V^{-1}(\mathbf x-\boldsymbol\mu)=1.54. \tag{7.4}\]

Число \(D^2\) безразмерно. Оно учитывает и масштабы двух ошибок, и направление, в котором они коррелируют. Контур, проходящий через нашу точку, охватывает около \(54\%\) двумерного гауссова распределения. Для области с вероятностью \(68.3\%\) потребуется \(D^2\leq2.30\).

Откуда взялись эти числа и почему привычное условие \(D^2\leq1\) здесь даёт не \(68.3\%\), а только \(39.3\%\)? Для ответа надо построить двумерный гаусс и посмотреть на геометрию его ковариационной матрицы.

7.2 От двух одномерных гауссов к двумерному

7.2.1 Независимые величины

Пусть

\[ X\sim\mathcal N(\mu_X,\sigma_X^2), \qquad Y\sim\mathcal N(\mu_Y,\sigma_Y^2) \tag{7.5}\]

и \(X\) и \(Y\) независимы. Совместная плотность равна произведению двух одномерных плотностей:

\[ f(x,y)= \frac{1}{2\pi\sigma_X\sigma_Y} \exp\!\left[-\frac12\left( \frac{(x-\mu_X)^2}{\sigma_X^2} +\frac{(y-\mu_Y)^2}{\sigma_Y^2} \right)\right]. \tag{7.6}\]

Введём безразмерные отклонения

\[ u=\frac{x-\mu_X}{\sigma_X}, \qquad v=\frac{y-\mu_Y}{\sigma_Y}. \tag{7.7}\]

Плотность постоянна при постоянном значении

\[ u^2+v^2=c. \tag{7.8}\]

В координатах \((u,v)\) это окружность. В исходных координатах получается эллипс с полуосями \(\sqrt c\,\sigma_X\) и \(\sqrt c\,\sigma_Y\). Его оси совпадают с координатными: отклонение \(X\) ничего не сообщает об отклонении \(Y\).

7.2.2 Корреляция поворачивает эллипс

Ковариационная матрица двух величин имеет вид

\[ V= \begin{pmatrix} \sigma_X^2 & \rho\sigma_X\sigma_Y\\ \rho\sigma_X\sigma_Y & \sigma_Y^2 \end{pmatrix}, \qquad |\rho|<1. \tag{7.9}\]

Смешанный член в совместной плотности появляется вместе с корреляцией:

\[ f(x,y)= \frac{1}{2\pi\sigma_X\sigma_Y\sqrt{1-\rho^2}} \exp\!\left[-\frac12 \frac{u^2-2\rho uv+v^2}{1-\rho^2} \right]. \tag{7.10}\]

При \(\rho=0\) уравнение 7.10 переходит в произведение двух одномерных гауссов. При положительной корреляции большие значения \(X\) чаще сопровождаются большими значениями \(Y\), поэтому эллипс наклонён вверх. Отрицательная корреляция меняет направление наклона.

Предел \(|\rho|=1\) в формулу плотности не входит. В этом случае определитель ковариационной матрицы обращается в нуль, а случайные точки лежат на прямой. Двумерной плотности относительно площади уже нет.

7.2.3 Обратная ковариационная матрица

Для матрицы 7.9

\[ \det V=\sigma_X^2\sigma_Y^2(1-\rho^2), \tag{7.11}\]

а обратная матрица равна

\[ V^{-1}= \frac{1}{1-\rho^2} \begin{pmatrix} 1/\sigma_X^2 & -\rho/(\sigma_X\sigma_Y)\\ -\rho/(\sigma_X\sigma_Y) & 1/\sigma_Y^2 \end{pmatrix}. \tag{7.12}\]

Обозначим отклонение от среднего через

\[ \boldsymbol\delta= \mathbf x-\boldsymbol\mu= \begin{pmatrix} x-\mu_X\\ y-\mu_Y \end{pmatrix}. \tag{7.13}\]

Тогда показатель экспоненты можно записать одной строкой:

\[ \boldsymbol\delta^{\mathsf T}V^{-1}\boldsymbol\delta =\frac{u^2-2\rho uv+v^2}{1-\rho^2}. \tag{7.14}\]

Формула через \(\rho\) и матричная формула описывают одну плотность. Матрица \(V\) задаёт флуктуации, а \(V^{-1}\) оценивает отклонение с учётом этих флуктуаций.

7.2.4 Проверка коэффициента \(\rho\)

Построим коррелированные величины из двух независимых стандартных гауссов:

\[ Z_1,Z_2\sim\mathcal N(0,1), \qquad \operatorname{cov}(Z_1,Z_2)=0, \tag{7.15}\]

\[ X=\mu_X+\sigma_XZ_1, \qquad Y=\mu_Y+\sigma_Y \left(\rho Z_1+\sqrt{1-\rho^2}\,Z_2\right). \tag{7.16}\]

Такой способ генерации пригодится и в методе Монте-Карло. Пока проверим, что параметры в нём имеют заявленный смысл.

7.2.4.1 Средние, дисперсии и ковариация

Из уравнения 7.16 сразу следуют средние

\[ \mathbb E[X]=\mu_X, \qquad \mathbb E[Y]=\mu_Y. \tag{7.17}\]

Дисперсия \(X\) равна \(\sigma_X^2\). Для второй величины получаем

\[ \begin{aligned} \operatorname{Var}(Y) &=\sigma_Y^2\left[ \rho^2\operatorname{Var}(Z_1) +(1-\rho^2)\operatorname{Var}(Z_2) \right]\\ &=\sigma_Y^2\left[\rho^2+(1-\rho^2)\right] =\sigma_Y^2. \end{aligned} \tag{7.18}\]

Смешанное среднее даёт

\[ \begin{aligned} \operatorname{cov}(X,Y) &=\sigma_X\sigma_Y \operatorname{cov}\!\left( Z_1,\rho Z_1+\sqrt{1-\rho^2}\,Z_2 \right)\\ &=\rho\sigma_X\sigma_Y. \end{aligned} \tag{7.19}\]

Поэтому параметр \(\rho\) в уравнении 7.10 равен коэффициенту корреляции:

\[ \rho=\frac{\operatorname{cov}(X,Y)}{\sigma_X\sigma_Y}. \tag{7.20}\]

7.3 Многомерный гаусс и расстояние Махаланобиса

Для случайного вектора размерности \(d\) со средним \(\boldsymbol\mu\) и невырожденной ковариационной матрицей \(V\) плотность имеет вид [1]

\[ f(\mathbf x)= \frac{1}{(2\pi)^{d/2}\sqrt{\det V}} \exp\!\left[-\frac12 (\mathbf x-\boldsymbol\mu)^{\mathsf T} V^{-1} (\mathbf x-\boldsymbol\mu) \right]. \tag{7.21}\]

Квадратичная форма

\[ D^2=(\mathbf x-\boldsymbol\mu)^{\mathsf T} V^{-1} (\mathbf x-\boldsymbol\mu) \tag{7.22}\]

называется квадратом расстояния Махаланобиса. Обычное евклидово расстояние считает один ГэВ по каждой оси одинаковым отклонением. Расстояние 7.22 сначала делит отклонения на их характерные масштабы и одновременно учитывает корреляции.

7.3.1 Главные оси

Поскольку ковариационная матрица симметрична, её можно разложить как

\[ V=R\Lambda R^{\mathsf T}, \qquad \Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_d). \tag{7.23}\]

Столбцы ортогональной матрицы \(R\) являются собственными векторами \(V\). Перейдём к координатам

\[ \mathbf z=\Lambda^{-1/2}R^{\mathsf T} (\mathbf x-\boldsymbol\mu). \tag{7.24}\]

В них ковариационная матрица единична, а расстояние принимает обычный вид

\[ D^2=\mathbf z^{\mathsf T}\mathbf z =\sum_{i=1}^{d}z_i^2. \tag{7.25}\]

Отсюда сразу читается геометрия контура \(D^2=c\). В двух измерениях это эллипс. Его главные оси направлены вдоль собственных векторов \(V\), а длины полуосей равны

\[ a_i=\sqrt{c\lambda_i}. \tag{7.26}\]

Коэффициент корреляции поворачивает собственные векторы относительно исходных осей. Собственные значения определяют ширины распределения в новых направлениях.

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

viewof rhoEllipse = Inputs.range([-0.98, 0.98], {
  value: 0.65,
  step: 0.01,
  label: "Корреляция ρ"
})

viewof sigmaXEllipse = Inputs.range([0.4, 2.2], {
  value: 1.4,
  step: 0.05,
  label: "Ширина σ_X"
})

viewof sigmaYEllipse = Inputs.range([0.4, 2.2], {
  value: 0.8,
  step: 0.05,
  label: "Ширина σ_Y"
})

viewof regenerateEllipse = Inputs.button("Новое облако", {
  value: 0,
  reduce: value => value + 1
})

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

randnEllipse = random => {
  const u1 = Math.max(random(), 1e-12);
  const u2 = random();
  return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
}

ellipseSeed = (
  20260906
  + 1009 * regenerateEllipse
  + 37 * Math.round((rhoEllipse + 1) * 1000)
  + 101 * Math.round(sigmaXEllipse * 100)
  + 211 * Math.round(sigmaYEllipse * 100)
) >>> 0

ellipseCloud = {
  const random = seededUniformCh7(ellipseSeed);
  const points = [];
  const residualScale = Math.sqrt(1 - rhoEllipse * rhoEllipse);
  for (let i = 0; i < 1200; ++i) {
    const z1 = randnEllipse(random);
    const z2 = randnEllipse(random);
    points.push({
      x: sigmaXEllipse * z1,
      y: sigmaYEllipse * (rhoEllipse * z1 + residualScale * z2)
    });
  }
  return points;
}

ellipseContour = {
  const radius = Math.sqrt(2.30);
  const residualScale = Math.sqrt(1 - rhoEllipse * rhoEllipse);
  return Array.from({length: 241}, (_, i) => {
    const angle = 2 * Math.PI * i / 240;
    const z1 = radius * Math.cos(angle);
    const z2 = radius * Math.sin(angle);
    return {
      x: sigmaXEllipse * z1,
      y: sigmaYEllipse * (rhoEllipse * z1 + residualScale * z2)
    };
  });
}

ellipseExtentX = 3.5 * sigmaXEllipse
ellipseExtentY = 3.5 * sigmaYEllipse

ellipsePlot = Plot.plot({
  width: 700,
  aspectRatio: 1,
  marginLeft: 62,
  x: {label: "X", domain: [-ellipseExtentX, ellipseExtentX], grid: true},
  y: {label: "Y", domain: [-ellipseExtentY, ellipseExtentY], grid: true},
  marks: [
    Plot.frame(),
    Plot.ruleX([0], {stroke: "#777"}),
    Plot.ruleY([0], {stroke: "#777"}),
    Plot.dot(ellipseCloud, {
      x: "x",
      y: "y",
      r: 2.1,
      fill: "#4ea1ff",
      fillOpacity: 0.34
    }),
    Plot.line(ellipseContour, {
      x: "x",
      y: "y",
      stroke: "#ff8c42",
      strokeWidth: 3.5
    })
  ]
})

html`
<div class="book-applet-panel">
  ${ellipsePlot}
  <div class="book-applet-summary">
    <span>ρ = ${rhoEllipse.toFixed(2)};</span>
    <span>σ_X = ${sigmaXEllipse.toFixed(2)};</span>
    <span>σ_Y = ${sigmaYEllipse.toFixed(2)};</span>
    <span>линия: D² = 2.30, внутри 68.3% распределения.</span>
  </div>
</div>
`

Меняйте \(\rho\), \(\sigma_X\) и \(\sigma_Y\). Точки показывают выборку из двумерного гаусса, линия — эллипс \(D^2=2.30\).

корреляция и эллипс

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

7.4 Сколько вероятности находится внутри эллипса

В одном измерении фраза «в пределах одной сигмы» понятна:

\[ P(|Z|\leq1)=0.683, \qquad Z\sim\mathcal N(0,1). \tag{7.27}\]

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

\[ D^2=Z_1^2+Z_2^2\sim\chi_2^2. \tag{7.28}\]

Плотность \(\chi_2^2\) особенно проста:

\[ f_{\chi_2^2}(q)=\frac12e^{-q/2}, \qquad q\geq0. \tag{7.29}\]

Интегрируя её от нуля до \(c\), получаем вероятность внутри эллипса:

\[ P(D^2\leq c)=1-e^{-c/2}. \tag{7.30}\]

Для заданной вероятности \(P\) граница находится из обратной формулы

\[ c=-2\ln(1-P). \tag{7.31}\]

Значение \(c=1\) в двух измерениях даёт

\[ P(D^2\leq1)=1-e^{-1/2}=0.393. \tag{7.32}\]

В таблице 7.1 приведены уровни, которые часто встречаются на двумерных графиках.

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

Вероятность внутри контура	Граница \(c\)
\(39.3\%\)	\(1.00\)
\(68.3\%\)	\(2.30\)
\(90\%\)	\(4.61\)
\(95\%\)	\(5.99\)
\(99\%\)	\(9.21\)

7.4.1 Интерактив: вероятность и размер контура

viewof contourProb = Inputs.range([0.1, 0.99], {
  value: 0.683,
  step: 0.001,
  label: "Вероятность внутри P"
})

viewof contourRho = Inputs.range([-0.95, 0.95], {
  value: -0.55,
  step: 0.05,
  label: "Корреляция ρ"
})

contourThreshold = -2 * Math.log(1 - contourProb)

contourCurveCh7 = {
  const radius = Math.sqrt(contourThreshold);
  const sigmaX = 1.5;
  const sigmaY = 0.8;
  const residualScale = Math.sqrt(1 - contourRho * contourRho);
  return Array.from({length: 241}, (_, i) => {
    const angle = 2 * Math.PI * i / 240;
    const z1 = radius * Math.cos(angle);
    const z2 = radius * Math.sin(angle);
    return {
      x: sigmaX * z1,
      y: sigmaY * (contourRho * z1 + residualScale * z2)
    };
  });
}

unitCurveCh7 = {
  const sigmaX = 1.5;
  const sigmaY = 0.8;
  const residualScale = Math.sqrt(1 - contourRho * contourRho);
  return Array.from({length: 241}, (_, i) => {
    const angle = 2 * Math.PI * i / 240;
    const z1 = Math.cos(angle);
    const z2 = Math.sin(angle);
    return {
      x: sigmaX * z1,
      y: sigmaY * (contourRho * z1 + residualScale * z2)
    };
  });
}

contourScaleCh7 = 1.25 * Math.sqrt(Math.max(contourThreshold, 1))
contourExtentXCh7 = contourScaleCh7 * 1.5
contourExtentYCh7 = contourScaleCh7 * 0.8

contourPlotCh7 = Plot.plot({
  width: 700,
  aspectRatio: 1,
  marginLeft: 62,
  x: {label: "X", domain: [-contourExtentXCh7, contourExtentXCh7], grid: true},
  y: {label: "Y", domain: [-contourExtentYCh7, contourExtentYCh7], grid: true},
  marks: [
    Plot.frame(),
    Plot.ruleX([0], {stroke: "#777"}),
    Plot.ruleY([0], {stroke: "#777"}),
    Plot.line(unitCurveCh7, {
      x: "x",
      y: "y",
      stroke: "#4ea1ff",
      strokeWidth: 3,
      strokeDasharray: "8,6"
    }),
    Plot.line(contourCurveCh7, {
      x: "x",
      y: "y",
      stroke: "#ff8c42",
      strokeWidth: 4
    })
  ]
})

html`
<div class="book-applet-panel">
  ${contourPlotCh7}
  <div class="book-applet-summary">
    <span>сплошная линия: P = ${(100 * contourProb).toFixed(1)}%;</span>
    <span>D² = ${contourThreshold.toFixed(2)};</span>
    <span>пунктир: D² = 1, внутри 39.3% распределения;</span>
    <span>ρ = ${contourRho.toFixed(2)}.</span>
  </div>
</div>
`

Сплошная линия ограничивает выбранную долю распределения. Пунктир показывает \(D^2=1\), который в двух измерениях содержит \(39.3\%\) вероятности.

вероятность и размер контура

7.4.2 Возвращение к калориметру

Для измерения 7.3 мы нашли \(D^2=1.54\). Уравнение 7.30 даёт

\[ P(D^2\leq1.54)=1-e^{-1.54/2}=0.537. \tag{7.33}\]

Иными словами, эллипс через эту точку содержит около \(54\%\) распределения. Точка лежит внутри области \(68.3\%\), поскольку \(1.54<2.30\). Положительная корреляция здесь существенна: оба измерения отклонились в одну сторону, то есть вдоль вытянутой оси облака.

7.4.3 Распределение данных и пространство параметров

Эллипсы этой главы описывают случайный вектор \(\mathbf X\) при известных \(\boldsymbol\mu\) и \(V\). В статьях похожие контуры часто рисуют в пространстве оцениваемых параметров вокруг полученного минимума. Одинаковая квадратичная форма ещё не делает их одним объектом.

Для параметров вероятностная интерпретация зависит от того, как построен интервал: от распределения оценки в частотном подходе или от апостериорной плотности в байесовском. Значения \(2.30\), \(4.61\) и \(5.99\) применимы к двумерной квадратичной гауссовой задаче. Условия, при которых ими можно пользоваться для параметров, мы обсудим вместе с правдоподобием и \(\chi^2\).

7.5 Объединение измерений

Пусть \(n\) приборов измеряют одну величину \(\mu\). Результаты собраны в вектор \(\mathbf x\), а их ковариационная матрица равна \(V\). Если совместное распределение гауссово, значение \(\mu\), при котором полученные данные имеют наибольшую плотность, минимизирует

\[ Q(\mu)= (\mathbf x-\mu\mathbf1)^{\mathsf T} V^{-1} (\mathbf x-\mu\mathbf1). \tag{7.34}\]

Здесь \(\mathbf1\) — вектор из единиц. Дифференцирование по \(\mu\) даёт

\[ \widehat\mu= \frac{\mathbf1^{\mathsf T}V^{-1}\mathbf x} {\mathbf1^{\mathsf T}V^{-1}\mathbf1}, \qquad \operatorname{Var}(\widehat\mu)= \frac{1}{\mathbf1^{\mathsf T}V^{-1}\mathbf1}. \tag{7.35}\]

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

\[ \widehat\mu= \frac{\displaystyle\sum_{i=1}^{n}x_i/\sigma_i^2} {\displaystyle\sum_{i=1}^{n}1/\sigma_i^2}, \qquad \operatorname{Var}(\widehat\mu)= \frac{1}{\displaystyle\sum_{i=1}^{n}1/\sigma_i^2}. \tag{7.36}\]

Для двух измерений с ковариацией \(c_{12}\) веса можно выписать явно:

\[ \widehat\mu=w_1x_1+w_2x_2, \qquad w_1=\frac{\sigma_2^2-c_{12}} {\sigma_1^2+\sigma_2^2-2c_{12}}, \qquad w_2=\frac{\sigma_1^2-c_{12}} {\sigma_1^2+\sigma_2^2-2c_{12}}. \tag{7.37}\]

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

7.6 Инвариантная масса двух частиц

Возьмём ещё один пример из лекции. Для двух безмассовых частиц с энергиями \(E_1\), \(E_2\) и углом \(\theta\) между импульсами

\[ m^2=2E_1E_2(1-\cos\theta). \tag{7.38}\]

Вектор измеряемых величин и его ковариационная матрица равны

\[ \mathbf x= \begin{pmatrix} E_1\\ E_2\\ \theta \end{pmatrix}, \qquad V_{\mathbf x}= \begin{pmatrix} 1.44 & 0.36 & 0\\ 0.36 & 1.00 & 0\\ 0 & 0 & 2.5\cdot10^{-5} \end{pmatrix}. \tag{7.39}\]

В энергетическом блоке числа указаны в \(\text{ГэВ}^2\), угловая дисперсия — в \(\text{рад}^2\). Возьмём

\[ E_1=60.0\ \text{ГэВ}, \qquad E_2=40.0\ \text{ГэВ}, \qquad \theta=0.200\ \text{рад}. \tag{7.40}\]

Градиент квадрата массы равен

\[ \nabla m^2= \begin{pmatrix} 2E_2(1-\cos\theta)\\ 2E_1(1-\cos\theta)\\ 2E_1E_2\sin\theta \end{pmatrix}. \tag{7.41}\]

7.6.0.1 Численный расчёт

Сначала найдём саму массу:

\[ m^2=2\cdot60\cdot40(1-\cos0.200) =95.680\ \text{ГэВ}^2, \qquad m=9.782\ \text{ГэВ}. \tag{7.42}\]

В заданной точке численный градиент равен

\[ \nabla m^2= \begin{pmatrix} 1.595\\ 2.392\\ 953.613 \end{pmatrix}. \tag{7.43}\]

По формуле распространения ковариации из главы 6

\[ \operatorname{Var}(m^2) \simeq (\nabla m^2)^{\mathsf T} V_{\mathbf x} (\nabla m^2) =34.864\ \text{ГэВ}^4. \tag{7.44}\]

Отсюда

\[ \sigma_{m^2}=5.905\ \text{ГэВ}^2, \qquad \sigma_m\simeq\frac{\sigma_{m^2}}{2m} =0.302\ \text{ГэВ}. \tag{7.45}\]

Результат можно записать как

\[ m=(9.78\pm0.30)\ \text{ГэВ}. \tag{7.46}\]

Ковариация двух энергий даёт в дисперсию \(m^2\) вклад

\[ 2\frac{\partial m^2}{\partial E_1} \frac{\partial m^2}{\partial E_2} \operatorname{cov}(E_1,E_2) =2.75\ \text{ГэВ}^4. \tag{7.47}\]

Если этот член отбросить, получится \(\sigma_m=0.290\) ГэВ вместо \(0.302\) ГэВ. Та же положительная корреляция, которая наклоняла эллипс двух энергий, увеличивает ошибку массы: обе производные по энергиям положительны.

7.7 Где заканчивается гауссова геометрия

Многомерный гаусс полностью определяется средним и ковариационной матрицей. Для произвольного распределения этих двух объектов уже недостаточно. Два распределения могут иметь одинаковые \(\boldsymbol\mu\) и \(V\), но разные хвосты и разную вероятность далёких отклонений.

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

Отдельная численная трудность возникает при \(|\rho|\simeq1\). Матрица \(V\) становится почти вырожденной, и результат обращения чувствителен к малым изменениям её элементов. Это означает, что два измерения содержат почти одну и ту же информацию. Печатать больше знаков в \(V^{-1}\) здесь бесполезно: следует проверить модель общей ошибки.

7.8 Итоги главы

• Ковариационная матрица задаёт масштабы и направления флуктуаций многомерного гаусса.
• Квадрат расстояния Махаланобиса учитывает единицы, дисперсии и корреляции.
• В двух измерениях область \(D^2\leq1\) содержит \(39.3\%\) вероятности; для \(68.3\%\) нужна граница \(D^2=2.30\).
• Объединение измерений и распространение ошибок требуют полной ковариационной матрицы.

7.9 Задачи

Задача 1. Вероятность внутри эллипса

Для двумерного гауссова распределения

\[ D^2=(\mathbf X-\boldsymbol\mu)^{\mathsf T} V^{-1}(\mathbf X-\boldsymbol\mu) \sim\chi_2^2. \]

1. Получите плотность \(D^2\) и покажите, что \[ P(D^2\leq c)=1-e^{-c/2}. \]
2. Найдите \(c\) для областей с вероятностями \(68.3\%\), \(90\%\) и \(95\%\).
3. Вычислите вероятность внутри контура \(D^2=1\).
4. Сравните этот результат с одномерной вероятностью \(P(|Z|\leq1)\) и объясните различие.

Задача 2. Геометрия ковариационной матрицы

Даны средний вектор и ковариационная матрица двух измеренных энергий:

\[ \boldsymbol\mu= \begin{pmatrix}60\\40\end{pmatrix}\ \text{ГэВ}, \qquad V= \begin{pmatrix} 1.44&0.36\\ 0.36&1.00 \end{pmatrix}\ \text{ГэВ}^2. \]

1. Найдите стандартные отклонения и коэффициент корреляции.
2. Вычислите собственные значения и собственные векторы \(V\).
3. Найдите длины и направления полуосей эллипса \(D^2=2.30\).
4. Для точки \(\mathbf x=(61.2,41.0)^{\mathsf T}\) ГэВ вычислите \(D^2\) и определите, лежит ли она внутри области с вероятностью \(68.3\%\).

Задача 3. Объединение двух измерений

Одна величина измерена двумя способами:

\[ x_1=10.2\pm0.8, \qquad x_2=11.1\pm0.4. \]

1. Найдите взвешенное среднее и его стандартное отклонение, считая измерения независимыми.
2. Повторите расчёт для коэффициента корреляции \(\rho=0.6\).
3. Выпишите оба веса во втором случае и объясните знак каждого из них.
4. Проверьте, что ковариационная матрица положительно определена.

Задача 4. Инвариантная масса

Для двух безмассовых частиц

\[ m^2=2E_1E_2(1-\cos\theta). \]

Пусть \(E_1=60\) ГэВ, \(E_2=40\) ГэВ, \(\theta=0.200\) рад, а ковариационная матрица величин \((E_1,E_2,\theta)\) равна

\[ V= \begin{pmatrix} 1.44&0.36&0\\ 0.36&1.00&0\\ 0&0&2.5\cdot10^{-5} \end{pmatrix}. \]

1. Найдите градиент \(m^2\) и вычислите \(\sigma_{m^2}\).
2. Получите \(\sigma_m\) линейным распространением ошибки.
3. Повторите расчёт при \(\operatorname{cov}(E_1,E_2)=0\) и сравните ответы.
4. Исследуйте предел \(\theta\ll1\): получите приближённые выражения для \(m\) и его производных.

Литература

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