6 Распространение ошибок

Аннотация

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

6.1 Физическая задача: измерение сечения

Эксперимент зарегистрировал $N=1000$ сигнальных событий. Эффективность отбора равна

\[ \varepsilon=0.80\pm0.02, \tag{6.1}\]

а интегральная светимость составляет

\[ L=(100\pm2)\ \text{фб}^{-1}. \tag{6.2}\]

В упрощённой модели фон отсутствует, число событий имеет распределение Пуассона, а эффективность и светимость измерены независимо от $N$ и друг от друга. Оценка сечения процесса равна

\[ \widehat\sigma_{\mathrm{proc}} =\frac{N}{\varepsilon L} =12.5\ \text{фб}. \tag{6.3}\]

Для пуассоновского счёта стандартное отклонение числа событий можно оценить как $\sigma_N=\sqrt N$. Относительная неопределённость сечения в линейном приближении тогда равна

\[ \begin{aligned} \left( \frac{\sigma_{\widehat\sigma_{\mathrm{proc}}}} {\widehat\sigma_{\mathrm{proc}}} \right)^2 &\simeq \left(\frac{\sigma_N}{N}\right)^2 +\left(\frac{\sigma_\varepsilon}{\varepsilon}\right)^2 +\left(\frac{\sigma_L}{L}\right)^2\\ &= \left(\frac{\sqrt{1000}}{1000}\right)^2 +\left(\frac{0.02}{0.80}\right)^2 +\left(\frac{2}{100}\right)^2 =0.045^2. \end{aligned} \tag{6.4}\]

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

\[ \widehat\sigma_{\mathrm{proc}} =(12.5\pm0.6)\ \text{фб}. \tag{6.5}\]

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

6.2 Что называется ошибкой

В названии главы используется привычное физикам выражение «распространение ошибок». В строгой метрологической терминологии ошибка — это разность между измеренным и опорным значениями, а неопределённость описывает разброс возможных результатов. В формулах этой главы речь идёт о стандартной неопределённости, которую мы обозначаем через $\sigma_X$.

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

6.3 Одна измеряемая величина

Пусть детектор измеряет случайную величину $X$, а результат вычисляется по формуле

\[ Y=f(X). \tag{6.6}\]

Обозначим

\[ \mathbb E[X]=\mu_X, \qquad \operatorname{Var}(X)=\sigma_X^2. \tag{6.7}\]

6.3.1 Линейное распространение ошибки

На масштабе $\sigma_X$ функцию $f$ можно заменить её касательной в точке $\mu_X$, если кривизна функции на этом интервале мала.

6.3.1.1 Формула для одной переменной

Первый порядок разложения Тейлора даёт

\[ f(X)\simeq f(\mu_X)+f'(\mu_X)(X-\mu_X). \tag{6.8}\]

В этом приближении

\[ \mathbb E[Y]\simeq f(\mu_X), \qquad Y-\mathbb E[Y] \simeq f'(\mu_X)(X-\mu_X). \tag{6.9}\]

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

\[ \begin{aligned} \operatorname{Var}(Y) &\simeq \mathbb E\!\left[ f'(\mu_X)^2(X-\mu_X)^2 \right]\\ &=f'(\mu_X)^2\sigma_X^2. \end{aligned} \tag{6.10}\]

Поэтому

\[ \boxed{ \sigma_Y\simeq |f'(\mu_X)|\sigma_X. } \tag{6.11}\]

Производная в уравнении 6.11 показывает, во сколько раз небольшое отклонение входной величины растягивается или сжимается при переходе к $Y$. Для линейной функции формула точна. Для нелинейной её точность определяется размером области, занятой распределением $X$.

6.3.2 Поперечный импульс в магнитном поле

Радиус $R$ траектории заряженной частицы в однородном магнитном поле связан с поперечным импульсом соотношением

\[ p_T=|q|BR. \tag{6.12}\]

В единицах, обычных для физики частиц,

\[ p_T[\text{ГэВ}/c] \simeq 0.2998\,|q/e|\,B[\text{Тл}]\,R[\text{м}]. \tag{6.13}\]

Если $q$ и $B$ считаются известными, неопределённость создаёт измерение радиуса. Поскольку зависимость линейна,

\[ \sigma_{p_T}=|q|B\sigma_R, \qquad \frac{\sigma_{p_T}}{p_T}=\frac{\sigma_R}{R}. \tag{6.14}\]

Ошибка радиуса переходит в такую же относительную ошибку импульса. Если неопределённость магнитного поля существенна, $B$ становится второй входной величиной. Такой случай разберём в разделе о ковариациях.

6.3.3 Степенная зависимость

Для функции

\[ Y=X^a \tag{6.15}\]

производная равна

\[ f'(x)=ax^{a-1}. \tag{6.16}\]

В точке $x=\mu_X$ уравнение 6.11 принимает вид

\[ \frac{\sigma_Y}{|f(\mu_X)|} \simeq |a|\frac{\sigma_X}{|\mu_X|}. \tag{6.17}\]

Несколько часто встречающихся функций собраны в таблице 6.1. Здесь $x_0=\mu_X$ и $y_0=f(x_0)$.

Таблица 6.1: Линейное распространение ошибки для часто встречающихся функций

Функция $f(x)$	Линейная неопределённость	Условие на точку разложения
$ax+b$	$\sigma_Y=\|a\|\sigma_X$	формула точна
$x^a$	$\sigma_Y/\|y_0\|\simeq \|a\|\sigma_X/\|x_0\|$	функция определена около $x_0$
$1/x$	$\sigma_Y/\|y_0\|\simeq\sigma_X/\|x_0\|$	$x_0$ далёк от нуля
$\ln x$	$\sigma_Y\simeq\sigma_X/x_0$	$x_0>0$, вероятность $X\le0$ мала
$e^x$	$\sigma_Y/y_0\simeq\sigma_X$	$X$ безразмерна

Таблица даёт локальные формулы. У отношения, логарифма и других функций с границей одной малости $\sigma_X$ недостаточно. Распределение $X$ должно с большой вероятностью оставаться в области, где функция гладкая.

6.4 Где линейное приближение перестаёт работать

Касательная описывает функцию в малой окрестности точки разложения. Размер этой окрестности задаёт $\sigma_X$. Линейная формула требует проверки в четырёх типичных ситуациях:

производная заметно меняется на интервале порядка $\sigma_X$;
$f'(\mu_X)$ близка к нулю и первым ненулевым оказывается квадратичный член;
рядом находится полюс, порог или физическая граница;
преобразованное распределение заметно асимметрично.

Возьмём $Y=X^2$ и $\mu_X=0$. Касательная к параболе в нуле горизонтальна, и линейная формула даёт $\sigma_Y\simeq0$. Сама величина $X^2$ при этом флуктуирует. Нулевой ответ сообщает о непригодности первого порядка разложения.

6.4.1 Нелинейное преобразование $Y=X+aX^2$

Рассмотрим модель из лекции:

\[ X\sim\mathcal N(1,\sigma_X^2), \qquad Y=X+aX^2. \tag{6.18}\]

Касательная в точке $\mu_X=1$ даёт

\[ \mathbb E[Y]\simeq1+a, \qquad \sigma_Y\simeq|1+2a|\sigma_X. \tag{6.19}\]

Для гауссовой величины моменты можно вычислить без приближения:

\[ \mathbb E[Y]=1+a+a\sigma_X^2, \tag{6.20}\]

\[ \operatorname{Var}(Y) =(1+2a)^2\sigma_X^2+2a^2\sigma_X^4. \tag{6.21}\]

Член $a\sigma_X^2$ сдвигает среднее, а $2a^2\sigma_X^4$ добавляет дисперсию. Оба исчезают в линейном приближении.

viewof propSigma = Inputs.range([0.05, 1.5], {value: 0.35, step: 0.05, label: "стандартное отклонение σ_X"})
viewof propCurvature = Inputs.range([0, 1.5], {value: 0.35, step: 0.05, label: "коэффициент нелинейности a"})
viewof regenerateProp = Inputs.button("Новая выборка", {
  value: 0,
  reduce: value => value + 1
})

mulberry32Prop = function(seed) {
  let state = seed >>> 0;
  return function() {
    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;
  };
}

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

propSeed = (
  20260905
  + Math.round(propSigma * 100) * 1009
  + Math.round(propCurvature * 100) * 9176
  + regenerateProp * 104729
) >>> 0

propSample = {
  const rng = mulberry32Prop(propSeed);
  const n = 18000;
  const values = new Array(n);
  for (let i = 0; i < n; ++i) {
    const x = 1 + propSigma * normal01Prop(rng);
    values[i] = x + propCurvature * x * x;
  }
  return values;
}

histProp = {
  const bins = 80;
  const lo = Math.min(...propSample);
  const hi = Math.max(...propSample);
  const width = (hi - lo) / bins;
  const counts = Array(bins).fill(0);
  for (const y of propSample) {
    const j = Math.min(bins - 1, Math.max(0, Math.floor((y - lo) / width)));
    counts[j] += 1;
  }
  return counts.map((count, j) => ({
    x0: lo + j * width,
    x1: lo + (j + 1) * width,
    density: count / (propSample.length * width)
  }));
}

propMeanLinear = 1 + propCurvature
propSigmaLinear = Math.abs(1 + 2 * propCurvature) * propSigma
propMeanExact = 1 + propCurvature + propCurvature * propSigma * propSigma
propSigmaExact = Math.sqrt(
  (1 + 2 * propCurvature) ** 2 * propSigma ** 2
  + 2 * propCurvature ** 2 * propSigma ** 4
)
normalProp = Array.from({length: 401}, (_, i) => {
  const lo = histProp[0].x0;
  const hi = histProp[histProp.length - 1].x1;
  const y = lo + (hi - lo) * i / 400;
  const z = (y - propMeanLinear) / propSigmaLinear;
  return {
    y,
    density: Math.exp(-0.5 * z * z) / (propSigmaLinear * Math.sqrt(2 * Math.PI))
  };
})

propActualMean = propSample.reduce((a, b) => a + b, 0) / propSample.length
propActualSigma = Math.sqrt(propSample.reduce((s, y) => s + (y - propActualMean) ** 2, 0) / propSample.length)

propPlot = Plot.plot({
  width: 940,
  height: 410,
  marginLeft: 70,
  x: {label: "Y = X + aX^2", grid: true},
  y: {label: "Плотность", grid: true},
  marks: [
    Plot.rectY(histProp, {x1: "x0", x2: "x1", y: "density", fill: "#4ea1ff", fillOpacity: 0.7}),
    Plot.line(normalProp, {x: "y", y: "density", stroke: "#ffb347", strokeWidth: 3})
  ]
})

html`
<div class="book-applet-panel">
  ${propPlot}
  <div class="book-applet-summary">
    <span> Линейное приближение: среднее ${propMeanLinear.toFixed(3)}, σ ${propSigmaLinear.toFixed(3)}; </span>
    <span> Точный расчёт: среднее ${propMeanExact.toFixed(3)}, σ ${propSigmaExact.toFixed(3)}; </span>
    <span> Монте-Карло: среднее ${propActualMean.toFixed(3)}, σ ${propActualSigma.toFixed(3)}; </span>
    <span> Синим показан результат Монте-Карло, оранжевым — линейное приближение.</span>
  </div>
</div>
`

Нелинейное преобразование $Y=X+aX^2$ — Нелинейное преобразование Y = X + aX²

На рисунке 6.1 показан режим, в котором кривизна уже существенна на масштабе входного распределения.

Рисунок 6.1: Распределение $Y=X+aX^2$ при $\sigma_X=0.9$ и $a=0.9$. Синяя гистограмма получена методом Монте-Карло, оранжевая кривая показывает гауссово распределение с параметрами линейного приближения.

Здесь расхождение видно по положению максимума, ширине и длинному правому хвосту. Запись $y\pm\sigma_Y$ сохраняет два числа, но уже не описывает форму распределения. В такой ситуации используют преобразование полного распределения или численное моделирование.

6.5 Разброс оценки стандартного отклонения

До сих пор $\mu$ и $\sigma$ были известными параметрами распределения. В эксперименте их оценивают по выборке

\[ X_1,\ldots,X_N\sim\mathcal N(\mu,\sigma^2), \tag{6.22}\]

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

\[ \widehat\mu=\overline X =\frac1N\sum_{k=1}^{N}X_k, \tag{6.23}\]

\[ \widehat{\sigma^2} =\frac1{N-1}\sum_{k=1}^{N}(X_k-\overline X)^2, \qquad \widehat\sigma=\sqrt{\widehat{\sigma^2}}. \tag{6.24}\]

Если повторить весь эксперимент, выборка изменится вместе со всеми тремя оценками. Поэтому у $\widehat\sigma$ есть собственное распределение и собственная дисперсия.

Для выборочного среднего

\[ \operatorname{Var}(\widehat\mu)=\frac{\sigma^2}{N}. \tag{6.25}\]

Параметр $\sigma^2$ неизвестен, и в оценке дисперсии среднего его заменяют выборочной дисперсией:

\[ \widehat{\operatorname{Var}}(\widehat\mu) =\frac{\widehat{\sigma^2}}{N}, \qquad \widehat\sigma_{\widehat\mu} =\frac{\widehat\sigma}{\sqrt N}. \tag{6.26}\]

Возникает следующий вопрос: насколько устойчиво значение $\widehat\sigma$, найденное по конечной выборке?

6.5.1 Почему в выборочной дисперсии стоит $N-1$

Введём отклонения от истинного среднего

\[ Y_k=X_k-\mu. \tag{6.27}\]

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

\[ \mathbb E[Y_k]=0, \qquad \mathbb E[Y_k^2]=\sigma^2. \tag{6.28}\]

Гауссовость для доказательства несмещённости выборочной дисперсии не нужна.

6.5.1.1 Несмещённость выборочной дисперсии

Обозначим две суммы

\[ A=\sum_{k=1}^{N}Y_k^2, \qquad B=\sum_{k=1}^{N}Y_k. \tag{6.29}\]

Разность выборочного и истинного среднего равна

\[ \overline X-\mu=\frac{B}{N}. \tag{6.30}\]

Сумма квадратов отклонений от выборочного среднего принимает вид

\[ \begin{aligned} T &=\sum_{k=1}^{N}(X_k-\overline X)^2\\ &=\sum_{k=1}^{N}\left(Y_k-\frac BN\right)^2\\ &=A-\frac{B^2}{N}. \end{aligned} \tag{6.31}\]

Независимость и центрирование дают

\[ \mathbb E[A]=N\sigma^2, \qquad \mathbb E[B^2]=N\sigma^2. \tag{6.32}\]

Отсюда

\[ \mathbb E[T] =N\sigma^2-\frac1N N\sigma^2 =(N-1)\sigma^2. \tag{6.33}\]

После деления на $N-1$

\[ \boxed{ \mathbb E[\widehat{\sigma^2}]=\sigma^2. } \tag{6.34}\]

Выборочные отклонения связаны условием

\[ \sum_{k=1}^{N}(X_k-\overline X)=0. \tag{6.35}\]

Одно из $N$ отклонений определяется через остальные. Поэтому сумма их квадратов содержит $N-1$ степеней свободы, и такой же делитель появляется в несмещённой оценке дисперсии.

6.5.2 Дисперсия оценки дисперсии

Теперь понадобится гауссовость. Для центрированной гауссовой величины

\[ \mathbb E[Y_k^4]=3\sigma^4. \tag{6.36}\]

Поскольку $\widehat{\sigma^2}=T/(N-1)$, достаточно найти дисперсию $T$.

6.5.2.1 Вычисление $\operatorname{Var}(\widehat{\sigma^2})$

Из уравнения 6.31

\[ T^2=A^2-\frac2N AB^2+\frac1{N^2}B^4. \tag{6.37}\]

Первый момент вычисляется раскрытием квадрата суммы:

\[ \begin{aligned} \mathbb E[A^2] &=\sum_{k=1}^{N}\mathbb E[Y_k^4] +2\sum_{i<j}\mathbb E[Y_i^2Y_j^2]\\ &=3N\sigma^4+N(N-1)\sigma^4\\ &=N(N+2)\sigma^4. \end{aligned} \tag{6.38}\]

Для смешанного момента выделим одно слагаемое, $B=Y_i+C_i$, где $C_i=\sum_{j\ne i}Y_j$. Величины $Y_i$ и $C_i$ независимы, $\mathbb E[C_i]=0$ и $\mathbb E[C_i^2]=(N-1)\sigma^2$. Поэтому

\[ \begin{aligned} \mathbb E[Y_i^2B^2] &=\mathbb E[Y_i^2(Y_i+C_i)^2]\\ &=3\sigma^4+(N-1)\sigma^4\\ &=(N+2)\sigma^4, \end{aligned} \tag{6.39}\]

и после суммирования по $i$

\[ \mathbb E[AB^2]=N(N+2)\sigma^4. \tag{6.40}\]

Сумма $B$ также гауссова, причём $\operatorname{Var}(B)=N\sigma^2$. Её четвёртый момент равен

\[ \mathbb E[B^4]=3N^2\sigma^4. \tag{6.41}\]

Подстановка трёх средних в уравнение 6.37 даёт

\[ \mathbb E[T^2]=(N^2-1)\sigma^4. \tag{6.42}\]

Вместе с уравнением 6.33 это приводит к

\[ \operatorname{Var}(T) =2(N-1)\sigma^4. \tag{6.43}\]

Искомая дисперсия оценки равна

\[ \boxed{ \operatorname{Var}(\widehat{\sigma^2}) =\frac{2\sigma^4}{N-1}. } \tag{6.44}\]

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

6.5.3 Разброс оценки стандартного отклонения

Оценка $\widehat\sigma$ связана с $\widehat{\sigma^2}$ функцией квадратного корня. Применим к ней формулу 6.11 в точке $\sigma^2$:

\[ f(y)=\sqrt y, \qquad f'(\sigma^2)=\frac1{2\sigma}. \tag{6.45}\]

Получаем

\[ \operatorname{Var}(\widehat\sigma) \simeq \frac1{4\sigma^2} \operatorname{Var}(\widehat{\sigma^2}) =\frac{\sigma^2}{2(N-1)}. \tag{6.46}\]

Относительный разброс оценки стандартного отклонения равен

\[ \boxed{ \frac{\sqrt{\operatorname{Var}(\widehat\sigma)}}{\sigma} \simeq \frac1{\sqrt{2(N-1)}}. } \tag{6.47}\]

Здесь снова использовано линейное приближение, теперь уже для квадратного корня. При конечном $N$ оценка $\widehat\sigma$ немного смещена вниз; с ростом выборки смещение исчезает, и формула становится точнее.

Для относительного разброса около $10\%$

\[ \frac1{\sqrt{2(N-1)}}\simeq0.1, \qquad N\simeq51. \tag{6.48}\]

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

6.5.4 Десять гауссовских измерений

viewof pars = Inputs.form({
  mu: Inputs.range([-5, 5], {value: 0.8, step: 0.1, label: "μ"}),
  sigma: Inputs.range([0.2, 5], {value: 1.5, step: 0.1, label: "σ"}),
  seed: Inputs.range([1, 1000], {value: 123, step: 1, label: "начальное состояние"})
})

muTrue = pars.mu
sigmaTrue = pars.sigma
seed = pars.seed
N = 10

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

normal01 = function(rng) {
  const u1 = Math.max(rng(), 1e-12)
  const u2 = rng()
  return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2)
}

rng = mulberry32(seed)
X = Array.from({length: N}, () => muTrue + sigmaTrue * normal01(rng))
xbar = X.reduce((a, b) => a + b, 0) / N
shat2 = X.map(x => (x - xbar) ** 2).reduce((a, b) => a + b, 0) / (N - 1)
shat = Math.sqrt(shat2)
varMeanHat = shat2 / N
sigmaMeanHat = shat / Math.sqrt(N)
varSHat = shat2 / (2 * (N - 1))
sigmaSHat = shat / Math.sqrt(2 * (N - 1))
fmt = x => Number(x).toFixed(3)

html`
<div class="book-applet-data">
  <div class="book-applet-grid">
    <div class="book-applet-table-wrap">
      <table class="book-applet-table">
        <thead><tr><th>k</th><th>X_k</th><th>X_k - X̄</th><th>(X_k - X̄)^2</th></tr></thead>
        <tbody>
          ${X.map((x, i) => html`<tr><td>${i + 1}</td><td>${fmt(x)}</td><td>${fmt(x - xbar)}</td><td>${fmt((x - xbar) ** 2)}</td></tr>`)}
        </tbody>
      </table>
    </div>
    <div class="book-applet-results">
      <div>Истинные параметры: μ = ${fmt(muTrue)}, σ = ${fmt(sigmaTrue)}, N = ${N}</div>
      <div class="book-applet-result"><div>Оценка среднего:</div><div>X̄ = ${fmt(xbar)}</div></div>
      <div class="book-applet-result"><div>Оценки дисперсии и стандартного отклонения:</div><div>σ̂² = ${fmt(shat2)}</div><div>σ̂ = ${fmt(shat)}</div></div>
      <div class="book-applet-result"><div>Неопределённость среднего:</div><div>σ̂ / √N = ${fmt(sigmaMeanHat)}</div><div>оценка Var(X̄) = ${fmt(varMeanHat)}</div></div>
      <div class="book-applet-result"><div>Неопределённость оценки σ̂:</div><div>σ̂ / √[2(N−1)] = ${fmt(sigmaSHat)}</div><div>оценка Var(σ̂) ≈ ${fmt(varSHat)}</div></div>
    </div>
  </div>
</div>
`

Десять чисел выглядят вполне убедительно, пока рядом не написаны истинные параметры. Новая выборка даёт другие $\widehat\mu$ и $\widehat\sigma$. Формулы выше описывают разброс этих результатов в серии повторений.

6.6 Несколько измеряемых величин

В задаче о сечении использовались три входные величины. Теперь выведем формулу, которая учитывает их совместные флуктуации.

6.6.1 Ковариация и ошибка суммы

Ковариация $X$ и $Y$ равна

\[ \operatorname{cov}(X,Y) =\mathbb E[(X-\mu_X)(Y-\mu_Y)]. \tag{6.49}\]

Для суммы $S=X+Y$

\[ \operatorname{Var}(S) =\sigma_X^2+\sigma_Y^2 +2\operatorname{cov}(X,Y). \tag{6.50}\]

Коэффициент корреляции убирает размерность ковариации:

\[ \rho_{XY} =\frac{\operatorname{cov}(X,Y)}{\sigma_X\sigma_Y}, \qquad -1\le\rho_{XY}\le1. \tag{6.51}\]

Для суммы и разности

\[ \sigma_{X\pm Y}^2 =\sigma_X^2+\sigma_Y^2 \pm2\rho_{XY}\sigma_X\sigma_Y. \tag{6.52}\]

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

Независимость влечёт $\rho_{XY}=0$. Обратное утверждение для произвольного распределения неверно: при нулевой корреляции нелинейная зависимость может сохраняться.

6.6.2 Оценка корреляции по выборке

Пусть получены $N$ независимых пар $(x_k,y_k)$ из одного совместного распределения. Выборочная ковариация равна

\[ \widehat{\operatorname{cov}}(X,Y) =\frac1{N-1}\sum_{k=1}^{N} (x_k-\overline x)(y_k-\overline y), \tag{6.53}\]

а выборочный коэффициент корреляции —

\[ \widehat\rho =\frac{\widehat{\operatorname{cov}}(X,Y)}{s_Xs_Y}. \tag{6.54}\]

При небольшом $N$ значение $\widehat\rho$ заметно меняется от выборки к выборке. В детекторных задачах ковариации получают также из моделирования, калибровочных данных и вариаций систематических параметров.

6.6.3 Общая формула распространения

Пусть

\[ F=f(X_1,\ldots,X_n), \tag{6.55}\]

а средние входных величин образуют вектор $\boldsymbol\mu=(\mu_1,\ldots,\mu_n)^{\mathsf T}$. Первый порядок разложения имеет вид

\[ F-\mathbb E[F] \simeq \sum_{i=1}^{n} \left. \frac{\partial f}{\partial x_i} \right|_{\boldsymbol\mu} (X_i-\mu_i). \tag{6.56}\]

6.6.3.1 Дисперсия функции нескольких переменных

Обозначим производные в точке $\boldsymbol\mu$ через

\[ g_i= \left. \frac{\partial f}{\partial x_i} \right|_{\boldsymbol\mu}. \tag{6.57}\]

После возведения линейного разложения в квадрат

\[ \begin{aligned} \operatorname{Var}(F) &\simeq \mathbb E\!\left[ \left(\sum_i g_i(X_i-\mu_i)\right) \left(\sum_j g_j(X_j-\mu_j)\right) \right]\\ &=\sum_{i,j}g_i \operatorname{cov}(X_i,X_j)g_j. \end{aligned} \tag{6.58}\]

Итак,

\[ \boxed{ \operatorname{Var}(F) \simeq \sum_{i,j} \frac{\partial f}{\partial x_i} \operatorname{cov}(X_i,X_j) \frac{\partial f}{\partial x_j}. } \tag{6.59}\]

Все производные в уравнении 6.59 вычисляются в одной точке $\boldsymbol\mu$. На практике неизвестные средние заменяют полученными центральными значениями.

Для двух переменных формула разворачивается в три слагаемых:

\[ \begin{aligned} \sigma_F^2\simeq{} &\left(\frac{\partial f}{\partial x}\right)^2\sigma_X^2 +\left(\frac{\partial f}{\partial y}\right)^2\sigma_Y^2\\ &+2\rho_{XY} \frac{\partial f}{\partial x} \frac{\partial f}{\partial y} \sigma_X\sigma_Y. \end{aligned} \tag{6.60}\]

Знак корреляционного вклада определяется одновременно знаком $\rho_{XY}$ и знаками двух производных.

6.6.4 Корреляционный вклад в дисперсию

viewof gradX = Inputs.range([-3, 3], {value: 1, step: 0.1, label: "∂f/∂x"})
viewof gradY = Inputs.range([-3, 3], {value: 1, step: 0.1, label: "∂f/∂y"})
viewof sigmaXCalc = Inputs.range([0.1, 2], {value: 1, step: 0.05, label: "σ_X"})
viewof sigmaYCalc = Inputs.range([0.1, 2], {value: 0.8, step: 0.05, label: "σ_Y"})
viewof rhoCalc = Inputs.range([-0.99, 0.99], {value: 0, step: 0.01, label: "ρ"})

termX = gradX * gradX * sigmaXCalc * sigmaXCalc
termY = gradY * gradY * sigmaYCalc * sigmaYCalc
termCorr = 2 * rhoCalc * gradX * gradY * sigmaXCalc * sigmaYCalc
varianceCalc = Math.max(0, termX + termY + termCorr)

contribData = [
  {term: "вклад X", value: termX, color: "#4ea1ff"},
  {term: "вклад Y", value: termY, color: "#ffb347"},
  {term: "корреляция", value: termCorr, color: "#b388ff"}
]

contribPlot = Plot.plot({
  width: 900,
  height: 350,
  marginLeft: 100,
  x: {label: "Вклад в дисперсию", grid: true},
  y: {label: null},
  marks: [
    Plot.ruleX([0], {stroke: "#ddd"}),
    Plot.barX(contribData, {x: "value", y: "term", fill: d => d.color})
  ]
})

html`
<div class="book-applet-panel">
  ${contribPlot}
  <div class="book-applet-summary">
    <span>σ_F² = ${varianceCalc.toFixed(3)}</span>
    <span>σ_F = ${Math.sqrt(varianceCalc).toFixed(3)}</span>
    <span>корреляционный вклад: ${termCorr.toFixed(3)}</span>
  </div>
</div>
`

На рисунке 6.2 производные имеют одинаковый знак, а корреляция отрицательна. Корреляционный член вычитается из двух положительных вкладов.

Рисунок 6.2: Вклады в дисперсию функции двух переменных при $\partial f/\partial x=\partial f/\partial y=1$, $\sigma_X=1$, $\sigma_Y=0.8$ и $\rho=-0.8$. Отрицательный корреляционный член уменьшает полную дисперсию до $0.36$.

6.6.5 Отношение двух величин

Для

\[ R=\frac XY \tag{6.61}\]

две производные равны

\[ \frac{\partial R}{\partial X}=\frac1Y, \qquad \frac{\partial R}{\partial Y}=-\frac{X}{Y^2}. \tag{6.62}\]

Линейная относительная неопределённость составляет

\[ \boxed{ \left(\frac{\sigma_R}{R}\right)^2 \simeq \left(\frac{\sigma_X}{X}\right)^2 +\left(\frac{\sigma_Y}{Y}\right)^2 -2\rho_{XY} \frac{\sigma_X}{X} \frac{\sigma_Y}{Y}. } \tag{6.63}\]

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

6.6.6 Возвращение к сечению

Для функции

\[ \widehat\sigma_{\mathrm{proc}}(N,\varepsilon,L) =\frac{N}{\varepsilon L} \tag{6.64}\]

производные удобно записать через само сечение:

\[ \frac{\partial\widehat\sigma_{\mathrm{proc}}}{\partial N} =\frac{\widehat\sigma_{\mathrm{proc}}}{N}, \qquad \frac{\partial\widehat\sigma_{\mathrm{proc}}}{\partial\varepsilon} =-\frac{\widehat\sigma_{\mathrm{proc}}}{\varepsilon}, \qquad \frac{\partial\widehat\sigma_{\mathrm{proc}}}{\partial L} =-\frac{\widehat\sigma_{\mathrm{proc}}}{L}. \tag{6.65}\]

Общая относительная дисперсия равна

\[ \begin{aligned} \left( \frac{\sigma_{\widehat\sigma_{\mathrm{proc}}}} {\widehat\sigma_{\mathrm{proc}}} \right)^2 \simeq{} &\frac{\sigma_N^2}{N^2} +\frac{\sigma_\varepsilon^2}{\varepsilon^2} +\frac{\sigma_L^2}{L^2}\\ &-2\frac{\operatorname{cov}(N,\varepsilon)}{N\varepsilon} -2\frac{\operatorname{cov}(N,L)}{NL} +2\frac{\operatorname{cov}(\varepsilon,L)}{\varepsilon L}. \end{aligned} \tag{6.66}\]

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

6.7 Ковариационная матрица

Для случайного вектора

\[ \mathbf X=(X_1,\ldots,X_n)^{\mathsf T} \tag{6.67}\]

определим матрицу

\[ (\mathbf V_X)_{ij} =\operatorname{cov}(X_i,X_j). \tag{6.68}\]

В двумерном случае

\[ \mathbf V_X= \begin{pmatrix} \sigma_X^2 & \rho\sigma_X\sigma_Y\\ \rho\sigma_X\sigma_Y & \sigma_Y^2 \end{pmatrix}. \tag{6.69}\]

Ковариационная матрица симметрична. На диагонали стоят дисперсии, а для любого вектора $\mathbf a$

\[ \mathbf a^{\mathsf T}\mathbf V_X\mathbf a =\operatorname{Var}(\mathbf a^{\mathsf T}\mathbf X) \ge0. \tag{6.70}\]

Последнее равенство объясняет требование положительной полуопределённости: никакая линейная комбинация измерений не может иметь отрицательную дисперсию. Матрица с нарушенным условием не описывает совместное распределение случайных величин.

6.7.1 Матричная запись распространения ошибки

Соберём производные скалярной функции в градиент

\[ \mathbf g= \left. \nabla f \right|_{\boldsymbol\mu} = \begin{pmatrix} \partial f/\partial x_1\\ \vdots\\ \partial f/\partial x_n \end{pmatrix}_{\boldsymbol\mu}. \tag{6.71}\]

Уравнение 6.59 записывается одной строкой:

\[ \boxed{ \sigma_F^2\simeq \mathbf g^{\mathsf T}\mathbf V_X\mathbf g. } \tag{6.72}\]

Пусть выходных величин несколько,

\[ \mathbf Y=\mathbf f(\mathbf X), \tag{6.73}\]

а матрица Якоби имеет элементы

\[ J_{ai}= \left. \frac{\partial f_a}{\partial x_i} \right|_{\boldsymbol\mu}. \tag{6.74}\]

Тогда

\[ \boxed{ \mathbf V_Y \simeq \mathbf J\mathbf V_X\mathbf J^{\mathsf T}. } \tag{6.75}\]

Для линейного преобразования матрица $\mathbf J$ постоянна и уравнение 6.75 является точным. Для нелинейной функции это первый порядок разложения около $\boldsymbol\mu$.

6.7.2 Сумма и разность

Возьмём

\[ \mathbf X= \begin{pmatrix}X\\Y\end{pmatrix}, \qquad \mathbf U= \begin{pmatrix}S\\D\end{pmatrix} = \begin{pmatrix}X+Y\\X-Y\end{pmatrix}. \tag{6.76}\]

Матрица Якоби равна

\[ \mathbf J= \begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix}. \tag{6.77}\]

Для входной матрицы

\[ \mathbf V_X= \begin{pmatrix} \sigma_X^2&\sigma_{XY}\\ \sigma_{XY}&\sigma_Y^2 \end{pmatrix} \tag{6.78}\]

получаем

\[ \mathbf V_U =\mathbf J\mathbf V_X\mathbf J^{\mathsf T} = \begin{pmatrix} \sigma_X^2+\sigma_Y^2+2\sigma_{XY} &\sigma_X^2-\sigma_Y^2\\ \sigma_X^2-\sigma_Y^2 &\sigma_X^2+\sigma_Y^2-2\sigma_{XY} \end{pmatrix}. \tag{6.79}\]

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

\[ \operatorname{cov}(S,D)=\sigma_X^2-\sigma_Y^2. \tag{6.80}\]

При $\sigma_X=\sigma_Y$ сумма и разность некоррелированы даже при $\sigma_{XY}\ne0$. В следующей главе мы увидим геометрический смысл этого результата на двумерном гауссовом распределении.

6.8 Итоги главы

При малой входной неопределённости функция заменяется касательной, и для одной переменной $\sigma_Y\simeq|f'(\mu_X)|\sigma_X$.
Выборочная дисперсия с делителем $N-1$ является несмещённой. Для гауссовой выборки её дисперсия равна $2\sigma^4/(N-1)$.
Корреляции входят в ошибку результата со знаками, заданными производными. Квадратурная сумма применима при нулевых ковариациях.
Ковариационная матрица преобразуется по закону $\mathbf V_Y\simeq\mathbf J\mathbf V_X\mathbf J^{\mathsf T}$; для линейной функции это равенство точное.

6.9 Задачи

Задача 1. Радиус и площадь круга

Измеренный радиус круга равен \[ R=(10.0\pm0.3)\ \text{см}. \]

Найдите площадь $S=\pi R^2$ и её ошибку в линейном приближении.
Сравните относительные ошибки $S$ и $R$.
Оцените поправку к математическому ожиданию площади, возникающую из-за квадратичности функции.

Задача 2. Сумма и разность коррелированных величин

Пусть \[ \sigma_X=2, \qquad \sigma_Y=3, \qquad \rho_{XY}=0.8. \]

Вычислите ошибки: \[ S=X+Y, \qquad D=X-Y. \]

Повторите вычисление для $\rho_{XY}=-0.8$ и объясните изменение результата.

Задача 3. Ошибка отношения

Эффективность определяется отношением \[ \varepsilon=\frac{N_{\mathrm{selected}}}{N_{\mathrm{total}}}. \]

Выведите линейную формулу для $\sigma_\varepsilon$, предполагая известные $\sigma_{\mathrm{selected}}$, $\sigma_{\mathrm{total}}$ и корреляцию $\rho$.
Объясните, почему предположение о независимости числителя и знаменателя обычно неверно.
Сравните ответы для $\rho=0$ и $\rho>0$.

Задача 4. Преобразование ковариационной матрицы

Пусть \[ \mathbf X= \begin{pmatrix}X\\Y\end{pmatrix}, \qquad V_X= \begin{pmatrix} 4&1.2\\ 1.2&1 \end{pmatrix}. \]

Для новых величин \[ S=X+Y, \qquad D=X-Y \] найдите:

матрицу Якоби;
ковариационную матрицу $V_{SD}$;
коэффициент корреляции между $S$ и $D$.

Задача 5. Проверка линейного приближения

Пусть \[ X\sim\mathrm N(1,\sigma_X^2), \qquad Y=X+aX^2. \]

Получите линейную оценку $\sigma_Y$.
Выведите точные $\mathbb E[Y]$ и $\operatorname{Var}(Y)$.
Сгенерируйте не менее $10^5$ значений $X$ и сравните выборочные моменты $Y$ с точным и линейным результатами.
На сетке $0\le a\le1$ и $0.05\le\sigma_X\le1$ найдите область, в которой линейная оценка стандартного отклонения отличается от точной менее чем на $5\%$.

Литература

[1] Joint Committee for Guides in Metrology. Evaluation of Measurement Data—Guide to the Expression of Uncertainty in Measurement. 2008. doi:10.59161/JCGM100-2008E.

Функция \(f(x)\)	Линейная неопределённость	Условие на точку разложения
\(ax+b\)	\(\sigma_Y=\|a\|\sigma_X\)	формула точна
\(x^a\)	\(\sigma_Y/\|y_0\|\simeq \|a\|\sigma_X/\|x_0\|\)	функция определена около \(x_0\)
\(1/x\)	\(\sigma_Y/\|y_0\|\simeq\sigma_X/\|x_0\|\)	\(x_0\) далёк от нуля
\(\ln x\)	\(\sigma_Y\simeq\sigma_X/x_0\)	\(x_0>0\), вероятность \(X\le0\) мала
\(e^x\)	\(\sigma_Y/y_0\simeq\sigma_X\)	\(X\) безразмерна