Двумерный гаусс

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

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

Двумерный гаусс

Независимые гауссовы величины

Для независимых \[ X\sim\mathrm N(\mu_X,\sigma_X^2), \qquad Y\sim\mathrm N(\mu_Y,\sigma_Y^2) \] совместная плотность является произведением: \[ \boxed{ 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]. } \]

Линии постоянной плотности — эллипсы с осями вдоль координат.

Догадка: наклонённый эллипс

Введём безразмерные отклонения \[ u=\frac{x-\mu_X}{\sigma_X}, \qquad v=\frac{y-\mu_Y}{\sigma_Y}. \]

Для независимых величин в экспоненте стоит \[ u^2+v^2. \]

Если величины коррелированы, эллипс должен повернуться. Самая простая квадратичная форма со смешанным членом: \[ \boxed{ Q_\rho(u,v) = \frac{u^2-2\rho uv+v^2}{1-\rho^2}, \qquad |\rho|<1. } \]

При \(\rho=0\) возвращается независимый случай.

Двумерная формула через \(\rho\)

Догадка для плотности: \[ \boxed{ 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]. } \]

Смешанный член \(-2\rho uv\) поворачивает эллипс.
Множитель \(\sqrt{1-\rho^2}\) нужен для нормировки.
Знак \(\rho\) задаёт направление наклона облака точек.

Это та же обратная ковариационная матрица

Пусть \[ V= \begin{pmatrix} \sigma_X^2&\rho\sigma_X\sigma_Y\\ \rho\sigma_X\sigma_Y&\sigma_Y^2 \end{pmatrix}. \]

Тогда \[ \det V = \sigma_X^2\sigma_Y^2(1-\rho^2), \] и \[ 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}. \]

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

Обозначим \[ \boldsymbol\delta = \begin{pmatrix} x-\mu_X\\ y-\mu_Y \end{pmatrix}. \]

Тогда \[ \boxed{ \boldsymbol\delta^\mathsf T V^{-1} \boldsymbol\delta = \frac{u^2-2\rho uv+v^2}{1-\rho^2} = Q_\rho(u,v). } \]

Значит формула через \(\rho\) и формула через \(V^{-1}\) — одна и та же запись.

Проверка догадки: восстановим ковариацию

Возьмём независимые стандартные гауссовы величины \[ Z_1,Z_2\sim\mathrm N(0,1), \qquad \operatorname{Cov}(Z_1,Z_2)=0. \]

Построим \[ X=\mu_X+\sigma_XZ_1, \qquad Y=\mu_Y+\sigma_Y\left(\rho Z_1+\sqrt{1-\rho^2}Z_2\right). \]

Ковариация действительно такая

Сразу получаем \[ \mathbb E[X]=\mu_X, \qquad \mathbb E[Y]=\mu_Y, \] \[ \operatorname{Var}(X)=\sigma_X^2, \qquad \operatorname{Var}(Y)=\sigma_Y^2 \left[\rho^2+(1-\rho^2)\right] = \sigma_Y^2, \] и \[ \boxed{ \operatorname{Cov}(X,Y) = \rho\sigma_X\sigma_Y. } \]

Значит \(\rho\) — именно коэффициент корреляции: \[ \rho = \frac{\operatorname{Cov}(X,Y)}{\sigma_X\sigma_Y}. \]

Шаг к произвольной размерности

В матричной форме двумерная формула сразу обобщается на \(d\) измерений: \[ \boxed{ f(\mathbf x) = \frac{1}{(2\pi)^{d/2}\sqrt{\det V}} \exp\!\left[ -\frac12 (\mathbf x-\boldsymbol\mu)^T V^{-1} (\mathbf x-\boldsymbol\mu) \right]. } \]

Для \(d=2\) это ровно предыдущая формула.

Величина \[ \boxed{ D^2 = (\mathbf x-\boldsymbol\mu)^T V^{-1} (\mathbf x-\boldsymbol\mu) } \] называется квадратом расстояния Махаланобиса.

Геометрия ковариационной матрицы

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

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

viewof rhoEllipse = Inputs.range([-0.98, 0.98], {value: 0.65, step: 0.02, label: "rho"})
viewof sigmaXEllipse = Inputs.range([0.4, 2.2], {value: 1.4, step: 0.05, label: "sigma_X"})
viewof sigmaYEllipse = Inputs.range([0.4, 2.2], {value: 0.8, step: 0.05, label: "sigma_Y"})
viewof regenerateEllipse = Inputs.button("Новое облако")

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

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

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

ellipsePlot = Plot.plot({
  width: 860,
  height: 500,
  marginLeft: 70,
  x: {label: "X", domain: [-4.5, 4.5], grid: true},
  y: {label: "Y", domain: [-4, 4], grid: true},
  marks: [
    Plot.ruleX([0], {stroke: "#777"}),
    Plot.ruleY([0], {stroke: "#777"}),
    Plot.dot(ellipseCloud, {x: "x", y: "y", r: 2.2, fill: "#4ea1ff", fillOpacity: 0.35}),
    Plot.line(ellipseContour, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 4})
  ]
})

html`
<div class="clt-panel">
  ${ellipsePlot}
  <div class="clt-summary">
    <span>rho = ${rhoEllipse.toFixed(2)}</span>
    <span>оранжевый эллипс: вероятность 68.3%</span>
    <span>Delta chi^2 = 2.30</span>
  </div>
</div>
`

Эллипсы и вероятностные контуры

Одномерный и двумерный случаи различаются

Для стандартного одномерного гаусса: \[ P(|Z|\le1)=0.683. \]

Но для двумерного стандартного гаусса область \[ X^2+Y^2\le1 \] содержит только \[ \boxed{ P(D^2\le1)=1-e^{-1/2}\approx0.393. } \]

Контур «одна сигма по радиусу» в двух измерениях не является областью с вероятностью 68.3%.

Вероятность внутри двумерного контура

Для двумерного гаусса \[ D^2 = (\mathbf x-\boldsymbol\mu)^TV^{-1}(\mathbf x-\boldsymbol\mu) \sim\chi^2_2. \]

Поэтому \[ \boxed{ P(D^2\le c)=1-e^{-c/2}. } \]

Обратно: \[ \boxed{ c=-2\ln(1-P). } \]

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

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

Числа зависят от размерности пространства параметров.

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

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: "Корр. ρ"})

contourC = -2 * Math.log(1 - contourProb)
contourCurve = {
  const q = Math.sqrt(contourC);
  const sx = 1.5;
  const sy = 0.8;
  const s = Math.sqrt(1 - contourRho * contourRho);
  return Array.from({length: 241}, (_, i) => {
    const a = 2 * Math.PI * i / 240;
    const z1 = q * Math.cos(a);
    const z2 = q * Math.sin(a);
    return {
      x: sx * z1,
      y: sy * (contourRho * z1 + s * z2)
    };
  });
}

unitCurve = {
  const q = 1;
  const sx = 1.5;
  const sy = 0.8;
  const s = Math.sqrt(1 - contourRho * contourRho);
  return Array.from({length: 241}, (_, i) => {
    const a = 2 * Math.PI * i / 240;
    const z1 = q * Math.cos(a);
    const z2 = q * Math.sin(a);
    return {
      x: sx * z1,
      y: sy * (contourRho * z1 + s * z2)
    };
  });
}

contourPlot = Plot.plot({
  width: 850,
  height: 500,
  marginLeft: 70,
  x: {label: "X", domain: [-5, 5], grid: true},
  y: {label: "Y", domain: [-4, 4], grid: true},
  marks: [
    Plot.ruleX([0], {stroke: "#777"}),
    Plot.ruleY([0], {stroke: "#777"}),
    Plot.line(unitCurve, {x: "x", y: "y", stroke: "#8ec5ff", strokeWidth: 3, strokeDasharray: "8,6"}),
    Plot.line(contourCurve, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 5})
  ]
})

html`
<div class="clt-panel">
  ${contourPlot}
  <div class="clt-summary">
    <span>оранжевый: P = ${(100 * contourProb).toFixed(1)}%</span>
    <span>Delta chi^2 = ${contourC.toFixed(2)}</span>
    <span>голубой пунктир: Delta chi^2 = 1, P = 39.3%</span>
  </div>
</div>
`

Контур распределения и контур параметров

В этой лекции эллипс описывает распределение известного случайного вектора \(\mathbf X\).
В задаче оценивания параметров контур строят вокруг best fit в пространстве параметров.
Формулы часто похожи, потому что в обоих случаях возникает квадратичная форма.
Но вероятностная интерпретация контуров параметров требует дополнительных условий.

Доверительные контуры параметров подробно рассмотрим после likelihood и метода \(\chi^2\).

Практические применения

Усреднение независимых измерений

Пусть независимые измерения одной величины равны \[ x_i\pm\sigma_i. \]

Оптимальное взвешенное среднее: \[ \boxed{ \widehat\mu = \frac{\sum_i x_i/\sigma_i^2} {\sum_i1/\sigma_i^2}, \qquad \sigma_{\widehat\mu}^2 = \frac{1}{\sum_i1/\sigma_i^2}. } \]

Более точные измерения получают больший вес.

Усреднение коррелированных измерений

Для вектора измерений \(\mathbf x\) с ковариационной матрицей \(V\): \[ \boxed{ \widehat\mu = \frac{\mathbf 1^TV^{-1}\mathbf x} {\mathbf 1^TV^{-1}\mathbf 1}, \qquad \sigma_{\widehat\mu}^2 = \frac{1}{\mathbf 1^TV^{-1}\mathbf 1}. } \]

Это обобщение взвешенного среднего.
Коррелированная информация не должна учитываться несколько раз.
При сильных корреляциях некоторые веса могут оказаться отрицательными.

Пример: инвариантная масса

Для двух безмассовых частиц \[ m^2=2E_1E_2(1-\cos\theta). \]

Удобно сначала распространить ошибку на \(m^2\): \[ \sigma_{m^2}^2 \approx (\nabla m^2)^TV(\nabla m^2), \] где \[ \nabla m^2 = \begin{pmatrix} 2E_2(1-\cos\theta)\\ 2E_1(1-\cos\theta)\\ 2E_1E_2\sin\theta \end{pmatrix}. \]

Затем \[ \sigma_m\approx\frac{\sigma_{m^2}}{2m}. \]

Пример: инвариантная масса

Пусть измерены энергии и угол двух для безмассовых частиц: \[ E_1 = 60.0~\mathrm{GeV}, \qquad E_2 = 40.0~\mathrm{GeV}, \qquad \theta = 0.200~\mathrm{rad}. \]
Пусть ковариационная матрица величин \[ \mathbf x = \begin{pmatrix} E_1\\ E_2\\ \theta \end{pmatrix} \]
равна \[ V_x = \begin{pmatrix} 1.44 & 0.36 & 0\\ 0.36 & 1.00 & 0\\ 0 & 0 & 2.5\cdot 10^{-5} \end{pmatrix}. \]
Вычислить ошибку величин \(m^2\) и \(m\).

Пример: инвариантная масса

Сначала вычислим саму массу: \[ m^2 = 2\cdot 60\cdot 40\cdot (1-\cos 0.200) = 95.680~\mathrm{GeV}^2, \] и \[ m = \sqrt{m^2}=9.782~\mathrm{GeV}. \]
Градиент величины \(m^2\): \[ \nabla m^2 = \begin{pmatrix} 2E_2(1-\cos\theta)\\ 2E_1(1-\cos\theta)\\ 2E_1E_2\sin\theta \end{pmatrix} = \begin{pmatrix} 1.595\\ 2.392\\ 953.613 \end{pmatrix}. \]
Тогда \[ \operatorname{Var}(m^2) \approx (\nabla m^2)^T V_x(\nabla m^2)=34.864~\mathrm{GeV}^4. \]

Пример: инвариантная масса

Следовательно, \[ \sigma_{m^2} = \sqrt{34.864}~\mathrm{GeV}^2 = 5.905~\mathrm{GeV}^2. \]
Для массы: \[ \sigma_m \approx \frac{\sigma_{m^2}}{2m} = \frac{5.905}{2\cdot 9.782} = 0.302~\mathrm{GeV}. \]
Итог: \[ \boxed{ m^2 = 95.68 \pm 5.91~\mathrm{GeV}^2 } \] или \[ \boxed{ m = 9.78 \pm 0.30~\mathrm{GeV}. } \]

Аналитика или Монте-Карло?

Линейное распространение

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

Монте-Карло

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

Практическая проверка результата

Записать функцию результата \(\mathbf Y=\mathbf f(\mathbf X)\).
Перечислить все источники неопределённости.
Построить полную ковариационную матрицу \(V_X\).
Вычислить Якобиан \(J\).
Получить \(V_Y\approx JV_XJ^T\).
Проверить единицы измерения и предельные случаи.
При заметной нелинейности сравнить с Монте-Карло.
Сохранять и публиковать ковариационную матрицу вместе с результатом.

Типичные ошибки анализа

Складывать все неопределённости в квадратуре, игнорируя корреляции.
Путать нулевую корреляцию с независимостью.
Использовать относительные ошибки рядом с нулём.
Линеаризовать функцию около неправильной точки.
Считать симметричной ошибку после сильно нелинейного преобразования.
Называть контур \(\Delta\chi^2=1\) областью \(68\%\) независимо от размерности.
Публиковать только диагональные ошибки без ковариаций.

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

Резюме

Ошибки производных величин возникают из локальной чувствительности функции к исходным измерениям.
Ковариация необходима для корректного объединения коррелированных ошибок.
Общий закон линейного преобразования: \[ \boxed{V_Y\approx JV_XJ^T.} \]
Ковариационная матрица задаёт геометрию многомерного гаусса.
В двумерном случае корреляции проявляются как поворот эллипса.
Вероятность внутри контура зависит от размерности.
Нелинейные случаи необходимо проверять методом Монте-Карло.