5  Когда гаусс не возникает

Аннотация

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

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

5.1 Физическая задача: ионизационные потери в трековом детекторе

Заряженная частица проходит через трековый детектор и оставляет сигнал в \(n\) последовательных слоях. По измеренным зарядам

\[ Q_1,Q_2,\ldots,Q_n \tag{5.1}\]

нужно оценить её удельные потери энергии \(dE/dx\). Первая идея состоит в том, чтобы взять среднее:

\[ \overline Q_n=\frac1n\sum_{i=1}^{n}Q_i. \tag{5.2}\]

Если \(Q_i\) независимы и имеют конечную дисперсию \(\sigma_Q^2\), то

\[ \operatorname{Var}(\overline Q_n)=\frac{\sigma_Q^2}{n}. \tag{5.3}\]

Можно ожидать, что с ростом числа слоёв распределение среднего станет узким и гауссовым. Однако распределение заряда в одном тонком слое имеет длинный правый хвост. Иногда частица передаёт большую энергию одному электрону, и один \(Q_i\) оказывается намного больше остальных. Обычное среднее полностью сохраняет этот вклад.

В газовых трековых детекторах часто используют усечённое среднее. Заряды сначала располагают по возрастанию,

\[ Q_{(1)}\le Q_{(2)}\le\cdots\le Q_{(n)}, \tag{5.4}\]

а затем усредняют только \(m\) меньших значений:

\[ \overline Q_{\mathrm{tr}} =\frac1m\sum_{i=1}^{m}Q_{(i)}, \qquad m<n. \tag{5.5}\]

Так редкие большие потери не определяют оценку \(dE/dx\). Выбор доли отбрасываемых измерений зависит от детектора и оптимизируется по его разрешению. Причину длинного хвоста мы получим ниже из спектра передачи энергии электрону [1].

Этот пример сразу ставит правильный вопрос. Число слагаемых само по себе ещё не определяет форму суммы. Нужно знать распределение отдельных вкладов и связи между ними.

5.2 Что именно утверждает центральная предельная теорема

Пусть

\[ S_n=X_1+X_2+\cdots+X_n. \tag{5.6}\]

Для независимых одинаково распределённых величин с параметрами

\[ \mathbb E[X_i]=\mu, \qquad 0<\operatorname{Var}(X_i)=\sigma^2<\infty \tag{5.7}\]

центральная предельная теорема даёт

\[ Z_n= \frac{S_n-n\mu}{\sigma\sqrt n} =\frac{S_n/n-\mu}{\sigma/\sqrt n} \xrightarrow[n\to\infty]{d}\mathcal N(0,1). \tag{5.8}\]

В этой формулировке содержатся четыре условия, которые понадобится проверять.

• Рассматривается сумма случайных величин.
• Слагаемые независимы и имеют одно распределение.
• Их дисперсия конечна.
• Переход к пределу выполняется после центрирования и нормировки суммы.

Для независимых величин с разными распределениями существуют более общие варианты теоремы. Условия Линдеберга и Ляпунова запрещают одному слагаемому сохранять конечную долю полной дисперсии. Для некоторых зависимых последовательностей также доказаны свои предельные теоремы [2]. Но из одной фразы «слагаемых много» ни одна из них не следует.

Общая формула для дисперсии суммы имеет вид

\[ \operatorname{Var}(S_n) =\sum_{i=1}^{n}\operatorname{Var}(X_i) +2\sum_{i<j}\operatorname{cov}(X_i,X_j). \tag{5.9}\]

При независимых слагаемых ковариации равны нулю. Корреляции оставляют вторую сумму и могут полностью изменить зависимость ширины от \(n\). Если дисперсия \(X_i\) не существует, нормировка через \(\sqrt n\) лишается основания.

5.3 Распределение Коши

Стандартное распределение Коши имеет плотность

\[ f_X(x)=\frac{1}{\pi(1+x^2)}. \tag{5.10}\]

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

\[ F_X(x)=\frac12+\frac1\pi\arctan x. \tag{5.11}\]

Поэтому величину с распределением Коши можно получить из равномерной \(U\sim U(0,1)\) преобразованием

\[ X=F_X^{-1}(U) =\tan\!\left[\pi\left(U-\frac12\right)\right]. \tag{5.12}\]

Хвост плотности убывает как \(1/x^2\). Из-за этого

\[ \mathbb E[|X|]=\infty, \qquad \mathbb E[X^2]=\infty. \tag{5.13}\]

Математическое ожидание распределения Коши не существует. Симметричное главное значение интеграла равно нулю, но математическим ожиданием оно не является.

5.3.1 Среднее Коши не сужается

Пусть \(X_1,\ldots,X_n\) — независимые стандартные величины Коши. В главе 3 мы нашли характеристическую функцию их суммы:

\[ \varphi_{S_n}(t)=e^{-n|t|}, \qquad S_n=\sum_{i=1}^{n}X_i. \tag{5.14}\]

Это характеристическая функция распределения Коши с масштабом \(n\). Деление суммы на \(n\) возвращает исходный масштаб:

\[ \frac{S_n}{n}\sim\operatorname{Cauchy}(0,1). \tag{5.15}\]

Распределение среднего при \(n=1\) и \(n=200\) совпадает. Для гауссовых слагаемых ширина среднего уменьшилась бы в \(\sqrt{200}\) раз.

viewof nCauchy = Inputs.range([1, 200], {
  value: 1,
  step: 1,
  label: "Число слагаемых n"
})

viewof sampleCauchy = Inputs.range([1000, 50000], {
  value: 15000,
  step: 1000,
  label: "Число опытов N"
})

viewof regenerateCauchy = Inputs.button("Новая статистика", {
  value: 0,
  reduce: value => value + 1
})

seededUniformCh5 = (seed) => {
  let state = seed >>> 0;
  return () => {
    state += 0x6D2B79F5;
    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;
  };
}

randnCauchy = (random) => {
  let u = 0;
  let v = 0;
  while (u === 0) u = random();
  while (v === 0) v = random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

cauchyRandom = (random) => Math.tan(Math.PI * (random() - 0.5))

cauchyMeans = {
  const random = seededUniformCh5(
    51001 + 1009 * regenerateCauchy + 17 * nCauchy + sampleCauchy
  );
  const values = new Array(sampleCauchy);
  for (let i = 0; i < sampleCauchy; ++i) {
    let sum = 0;
    for (let j = 0; j < nCauchy; ++j) sum += cauchyRandom(random);
    values[i] = sum / nCauchy;
  }
  return values;
}

gaussMeans = {
  const random = seededUniformCh5(
    71011 + 1013 * regenerateCauchy + 19 * nCauchy + sampleCauchy
  );
  const values = new Array(sampleCauchy);
  for (let i = 0; i < sampleCauchy; ++i) {
    let sum = 0;
    for (let j = 0; j < nCauchy; ++j) sum += randnCauchy(random);
    values[i] = sum / Math.sqrt(nCauchy);
  }
  return values;
}

histogramCauchy = (values, bins, min, max) => {
  const width = (max - min) / bins;
  const counts = Array(bins).fill(0);
  for (const value of values) {
    if (value < min || value >= max) continue;
    counts[Math.floor((value - min) / width)] += 1;
  }
  return counts.map((count, i) => ({
    x0: min + i * width,
    x1: min + (i + 1) * width,
    density: count / (values.length * width)
  }));
}

histCauchy = histogramCauchy(cauchyMeans, 100, -10, 10)
histGaussMean = histogramCauchy(gaussMeans, 100, -10, 10)

cauchyTheory = Array.from({length: 501}, (_, i) => {
  const x = -10 + 20 * i / 500;
  return {x, y: 1 / (Math.PI * (1 + x * x))};
})

cauchyPlot = Plot.plot({
  width: 760,
  height: 390,
  marginLeft: 55,
  x: {label: "Значение", domain: [-10, 10]},
  y: {label: "Плотность", domain: [0, 0.38]},
  grid: true,
  marks: [
    Plot.rectY(histGaussMean, {
      x1: "x0", x2: "x1", y: "density",
      fill: "#9be564", fillOpacity: 0.40
    }),
    Plot.rectY(histCauchy, {
      x1: "x0", x2: "x1", y: "density",
      fill: "#4ea1ff", fillOpacity: 0.68
    }),
    Plot.line(cauchyTheory, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 2.7})
  ]
})

html`
<div class="book-applet-panel">
  ${cauchyPlot}
  <div class="book-applet-summary">
    <span>Синий: среднее величин Коши;</span>
    <span>зелёный: нормированная сумма гауссовых величин;</span>
    <span>n = ${nCauchy}.</span>
  </div>
</div>
`

Распределение Коши встречается и в физике. Нерелятивистская линия Брейта—Вигнера имеет ту же форму и служит простой моделью резонансного пика. Для релятивистского резонанса форма зависит от выбранной кинематической переменной и фазового пространства.

Из поведения среднего не следует, что накопление данных для распределения Коши бесполезно. Параметр положения можно оценивать по медиане или методом максимального правдоподобия. Эти оценки используют форму распределения и не требуют существования \(\mathbb E[X]\).

5.4 Энергетические потери в тонком слое вещества

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

\[ \Delta=\omega_1+\omega_2+\cdots+\omega_N, \tag{5.16}\]

где \(\omega_a\) — энергия, переданная в одном столкновении. Формально перед нами сумма. Форма распределения определяется хвостом отдельных \(\omega_a\).

5.4.1 Редкие большие передачи энергии

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

\[ \frac{d\sigma}{d\omega}\propto\frac1{\omega^2}, \qquad \omega\le W_{\max}. \tag{5.17}\]

Здесь \(W_{\max}\) — максимальная передаваемая энергия, заданная кинематикой столкновения. Вероятность передачи больше порога \(W\) пропорциональна

\[ P(\omega>W) \propto \int_W^{W_{\max}}\frac{d\omega}{\omega^2} =\frac1W-\frac1{W_{\max}}. \tag{5.18}\]

При \(W\ll W_{\max}\) получаем

\[ P(\omega>W)\propto\frac1W. \tag{5.19}\]

Хвост убывает медленно. Один \(\delta\)-электрон может унести заметную часть энергии, доступной в слое, и дать заряд намного больше наиболее вероятного. Именно эти события образуют правый хвост распределения \(Q_i\) в трековом детекторе.

5.4.2 Распределение Ландау

В идеализированном пределе тонкого слоя распределение потери энергии имеет ландау-подобную форму

\[ p(\Delta) =\frac1\xi L\!\left(\frac{\Delta-\Delta_0}{\xi}\right). \tag{5.20}\]

Параметр \(\Delta_0\) задаёт положение распределения, а \(\xi\) — масштаб флуктуаций. При используемом здесь определении стандартной функции Ландау \(\Delta_0\) не совпадает с наиболее вероятным значением. Для частицы с зарядом \(ze\) в веществе

\[ \xi= \frac K2\frac ZA\frac{z^2}{\beta^2}\rho x. \tag{5.21}\]

Здесь \(Z\) и \(A\) — атомный номер и молярная масса вещества, \(\rho x\) — массовая толщина слоя, \(\beta=v/c\), а \(K=4\pi N_A r_e^2m_ec^2\simeq0.307\, \text{МэВ}\,\text{см}^2\,\text{моль}^{-1}\).

Стандартную функцию Ландау можно определить интегралом

\[ L(\lambda)= \frac1\pi\int_0^\infty \exp[-u\ln u-\lambda u]\sin(\pi u)\,du. \tag{5.22}\]

При больших положительных \(\lambda\)

\[ L(\lambda)\sim\frac1{\lambda^2}. \tag{5.23}\]

Наиболее вероятное значение конечно, а среднее и дисперсия идеального распределения Ландау не определены. Подробное описание флуктуаций ионизационных потерь в тонких кремниевых детекторах дано в [1].

viewof sampleLandau = Inputs.range([1000, 50000], {
  value: 15000,
  step: 1000,
  label: "Число событий N"
})

viewof landauScale = Inputs.range([0.3, 2.5], {
  value: 1.0,
  step: 0.05,
  label: "Масштаб"
})

viewof landauShift = Inputs.range([-3.0, 3.0], {
  value: 0.0,
  step: 0.1,
  label: "Сдвиг"
})

viewof regenerateLandau = Inputs.button("Новая статистика", {
  value: 0,
  reduce: value => value + 1
})

randnLandau = () => {
  let u = 0;
  let v = 0;
  while (u === 0) u = Math.random();
  while (v === 0) v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

histogramLandau = (values, bins, min, max) => {
  const width = (max - min) / bins;
  const counts = Array(bins).fill(0);
  for (const value of values) {
    if (value < min || value >= max) continue;
    counts[Math.floor((value - min) / width)] += 1;
  }
  return counts.map((count, i) => ({
    x0: min + i * width,
    x1: min + (i + 1) * width,
    density: count / (values.length * width)
  }));
}

// Landau-подобная случайная величина.
// Это alpha-stable распределение с alpha = 1 и максимальной правой асимметрией.
// Оно имеет тот же физический смысл: длинный правый хвост и отсутствие обычной дисперсии.
landauStableRandom = (random) => {
  const beta = 1.0;
  const V = Math.PI * (random() - 0.5);
  let W = 0;
  while (W === 0) W = -Math.log(random());

  return (2 / Math.PI) * (
    (Math.PI / 2 + beta * V) * Math.tan(V)
    - beta * Math.log(
      (Math.PI / 2 * W * Math.cos(V)) /
      (Math.PI / 2 + beta * V)
    )
  );
}

landauValues = {
  const random = seededUniformCh5(
    91009
    + 1019 * regenerateLandau
    + sampleLandau
    + Math.round(100 * landauScale)
    + Math.round(100 * landauShift)
  );
  const values = new Array(sampleLandau);
  for (let i = 0; i < sampleLandau; ++i) {
    values[i] = landauShift + landauScale * landauStableRandom(random);
  }
  return values;
}

landauSorted = [...landauValues].sort((a, b) => a - b)

quantileLandau = (arr, q) => {
  const pos = q * (arr.length - 1);
  const lo = Math.floor(pos);
  const hi = Math.ceil(pos);
  const h = pos - lo;
  return arr[lo] * (1 - h) + arr[hi] * h;
}

landauMedian = quantileLandau(landauSorted, 0.50)
landauQ16 = quantileLandau(landauSorted, 0.16)
landauQ84 = quantileLandau(landauSorted, 0.84)
landauRobustSigma = 0.5 * (landauQ84 - landauQ16)

histLandau = histogramLandau(landauValues, 120, -8, 25)

normalForComparison = Array.from({length: 501}, (_, i) => {
  const x = -8 + 33 * i / 500;
  const z = (x - landauMedian) / landauRobustSigma;
  return {
    x,
    y: Math.exp(-0.5 * z * z) / (landauRobustSigma * Math.sqrt(2 * Math.PI))
  };
})

landauPlot = Plot.plot({
  width: 760,
  height: 390,
  marginLeft: 55,
  x: {label: "Потеря энергии Δ", domain: [-8, 25]},
  y: {label: "Плотность", domain: [0, 0.42]},
  grid: true,
  marks: [
    Plot.rectY(histLandau, {
      x1: "x0",
      x2: "x1",
      y: "density",
      fill: "#4ea1ff",
      fillOpacity: 0.70,
      stroke: "white"
    }),
    Plot.line(normalForComparison, {
      x: "x",
      y: "y",
      stroke: "#ffb347",
      strokeWidth: 2.8
    })
  ]
})

html`
<div class="book-applet-panel">
  ${landauPlot}
  <div class="book-applet-summary">
    <span>Синий: ландау-подобное распределение;</span>
    <span>оранжевый: гаусс с той же медианой и центральной шириной;</span>
    <span>медиана = ${landauMedian.toFixed(2)};</span>
    <span>центральная ширина = ${landauRobustSigma.toFixed(2)}.</span>
  </div>
</div>
`

В идеализированной модели верхнюю границу передачи энергии удаляют. Тогда хвост \(p(\omega)\propto1/\omega^2\) даёт логарифмически расходящийся первый момент и линейно расходящийся второй:

\[ \int^\infty \omega\frac{d\omega}{\omega^2} =\int^\infty\frac{d\omega}{\omega}, \qquad \int^\infty \omega^2\frac{d\omega}{\omega^2} =\int^\infty d\omega. \tag{5.24}\]

Обычная ЦПТ к такой модели неприменима. В реальном столкновении \(W_{\max}\) конечно. Распределение Ландау, столь любимое экспериментаторами для описания потерь энергии в тонких детекторах - это приближение. Если хотите точнее, используйте распределение Вавилова. При увеличении толщины слоя последнее переходит к гауссову распределению.

Теперь понятна и конструкция из начала главы. Усечённое среднее подавляет редкие большие \(Q_i\), созданные \(\delta\)-электронами. Цена этой операции — потеря части измерений и необходимость отдельно определить распределение и смещение оценки \(\overline Q_{\mathrm{tr}}\).

Биография

Лев Давидович Ландау

1908–1968

Лев Ландау - выдающийся советский физик-теоретик. Он разработал квантовую теорию диамагнетизма, теорию фазовых переходов второго рода, теорию сверхтекучести жидкого гелия и теорию ферми-жидкости. Вместе с Виталием Гинзбургом он построил феноменологическую теорию сверхпроводимости. За работы по теории конденсированного состояния, прежде всего жидкого гелия, Ландау получил Нобелевскую премию по физике 1962 года.

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

Лев Ландау и Евгений Лифшиц создали многотомный «Курс теоретической физики», воспитавший множество поколений учёных в нашей стране и во всём мире.

Источник: Nobel Prize; фото архив Нобелевского фонда, общественное достояние

5.5 Коррелированные слагаемые

Ковариационная сумма в уравнении 5.9 исчезает только для некоррелированных слагаемых. Рассмотрим простой случай:

\[ \operatorname{Var}(X_i)=\sigma^2, \qquad \operatorname{corr}(X_i,X_j)=\rho, \quad i\ne j. \tag{5.25}\]

Для среднего

\[ \overline X_n=\frac1n\sum_{i=1}^{n}X_i \tag{5.26}\]

получаем

\[ \begin{aligned} \operatorname{Var}(\overline X_n) &=\frac1{n^2}\left[ n\sigma^2+2\frac{n(n-1)}2\rho\sigma^2 \right]\\ &=\frac{\sigma^2}{n}\left[1+(n-1)\rho\right]. \end{aligned} \tag{5.27}\]

При \(\rho=0\) остаётся \(\sigma^2/n\). При фиксированном \(\rho>0\)

\[ \operatorname{Var}(\overline X_n) \xrightarrow[n\to\infty]{}\rho\sigma^2. \tag{5.28}\]

Матрица с одинаковой корреляцией всех пар допустима при \(-1/(n-1)\le\rho\le1\). Поэтому фиксированное отрицательное \(\rho\) нельзя сохранять при неограниченном росте \(n\). Положительная корреляция оставляет конечную ширину среднего.

5.5.1 Общий сдвиг как систематическая ошибка

Пусть каждое измерение содержит одну общую и одну независимую составляющую:

\[ X_i=\mu+C+\varepsilon_i. \tag{5.29}\]

В ансамбле повторных калибровок и измерений зададим

\[ C\sim\mathcal N(0,\sigma_C^2), \qquad \varepsilon_i\sim\mathcal N(0,\sigma_\varepsilon^2), \tag{5.30}\]

причём \(C\) и все \(\varepsilon_i\) независимы. Тогда

\[ \overline X_n =\mu+C+\frac1n\sum_{i=1}^{n}\varepsilon_i \tag{5.31}\]

и

\[ \operatorname{Var}(\overline X_n) =\sigma_C^2+\frac{\sigma_\varepsilon^2}{n}. \tag{5.32}\]

Статистический член уменьшается с числом измерений. Общий сдвиг повторяется во всех \(X_i\) и остаётся целиком. Так ведут себя, например, неопределённость энергетической шкалы, нормировки потока или светимости, коэффициента усиления и общий временной сдвиг.

viewof nCorr = Inputs.range([1, 500], {value: 10, step: 1, label: "Число измерений n"})
viewof sigmaCommon = Inputs.range([0, 2.0], {value: 0.6, step: 0.05, label: "Общая ошибка σ_C"})
viewof sigmaStat = Inputs.range([0.1, 2.0], {value: 1.0, step: 0.05, label: "Статистическая ошибка σ_ε"})
viewof sampleCorr = Inputs.range([1000, 30000], {value: 12000, step: 1000, label: "Число опытов N"})
viewof regenerateCorr = Inputs.button("Новая статистика", {
  value: 0,
  reduce: value => value + 1
})

corrMeans = {
  const random = seededUniformCh5(
    111031
    + 1021 * regenerateCorr
    + 23 * nCorr
    + sampleCorr
    + Math.round(100 * sigmaCommon)
    + Math.round(100 * sigmaStat)
  );
  const values = new Array(sampleCorr);
  for (let i = 0; i < sampleCorr; ++i) {
    const c = sigmaCommon * randnCauchy(random);
    let sum = 0;
    for (let j = 0; j < nCorr; ++j) {
      sum += c + sigmaStat * randnCauchy(random);
    }
    values[i] = sum / nCorr;
  }
  return values;
}

corrTheorySigma = Math.sqrt(sigmaCommon * sigmaCommon + sigmaStat * sigmaStat / nCorr)

histCorr = histogramCauchy(corrMeans, 90, -4, 4)
normalCorr = Array.from({length: 401}, (_, i) => {
  const z = -4 + 8 * i / 400;
  return {z, density: Math.exp(-0.5 * z * z / (corrTheorySigma * corrTheorySigma)) / (corrTheorySigma * Math.sqrt(2 * Math.PI))};
})

corrPlot = Plot.plot({
  width: 760,
  height: 390,
  marginLeft: 55,
  x: {label: "Среднее", domain: [-4, 4]},
  y: {label: "Плотность", domain: [0, 1.4]},
  grid: true,
  marks: [
    Plot.rectY(histCorr, {x1: "x0", x2: "x1", y: "density", fill: "#4ea1ff", fillOpacity: 0.70, stroke: "white"}),
    Plot.line(normalCorr, {x: "z", y: "density", stroke: "#ffb347", strokeWidth: 2.8})
  ]
})

html`
<div class="book-applet-panel">
  ${corrPlot}
  <div class="book-applet-summary">
    <span>Var(среднего) = σ_C² + σ_ε²/n;</span>
    <span>σ(среднего) = ${corrTheorySigma.toFixed(3)};</span>
    <span>n = ${nCorr}.</span>
  </div>
</div>
`

В этой модели \(\overline X_n\) имеет гауссово распределение при любом \(n\). Однако его ширина стремится к \(\sigma_C\) и не обращается в нуль. Гауссова форма гистограммы сама по себе не подтверждает закон \(1/\sqrt n\).

5.6 Смеси распределений

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

\[ N_1\sim\operatorname{Pois}(1000), \qquad N_2\sim\operatorname{Pois}(100) \tag{5.33}\]

независимы, то

\[ N_1+N_2\sim\operatorname{Pois}(1100). \tag{5.34}\]

В смеси каждое событие принадлежит одному из нескольких классов. Введём скрытую переменную \(K\in\{1,2\}\) и зададим

\[ N\mid K=1\sim\operatorname{Pois}(1000), \qquad N\mid K=2\sim\operatorname{Pois}(100). \tag{5.35}\]

Если класс события не зарегистрирован, в данных остаётся только \(N\). Даже узкие распределения внутри классов образуют широкую или многопиковую смесь.

5.6.1 Среднее и дисперсия смеси

Пусть \(X\) при фиксированном классе \(K\) имеет среднее \(\mu_K\) и дисперсию \(\sigma_K^2\). Полное среднее равно

\[ \mathbb E[X]=\mathbb E_K[\mu_K]. \tag{5.36}\]

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

\[ \operatorname{Var}(X) =\mathbb E_K\!\left[\sigma_K^2\right] +\operatorname{Var}_K(\mu_K). \tag{5.37}\]

Первое слагаемое усредняет дисперсии внутри классов. Второе измеряет разброс между их средними.

Для двух равновероятных пуассоновских классов из уравнения 5.35

\[ \mathbb E[N]=\frac{1000+100}{2}=550. \tag{5.38}\]

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

\[ \mathbb E_K[\operatorname{Var}(N\mid K)] =\frac{1000+100}{2}=550, \tag{5.39}\]

а дисперсия условных средних составляет

\[ \operatorname{Var}_K(\mathbb E[N\mid K]) =\frac{(1000-550)^2+(100-550)^2}{2} =202\,500. \tag{5.40}\]

Полная дисперсия

\[ \operatorname{Var}(N)=550+202\,500=203\,050 \tag{5.41}\]

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

viewof mixWeight = Inputs.range([0.0, 1.0], {value: 0.5, step: 0.02, label: "Доля первого класса w"})
viewof mixSeparation = Inputs.range([0.0, 6.0], {value: 2.5, step: 0.1, label: "Расстояние между средними"})
viewof mixSigma2 = Inputs.range([0.2, 2.5], {value: 1.0, step: 0.05, label: "Ширина второго класса σ₂"})
viewof sampleMix = Inputs.range([1000, 40000], {value: 15000, step: 1000, label: "Число событий N"})
viewof regenerateMix = Inputs.button("Новая статистика", {
  value: 0,
  reduce: value => value + 1
})

mixValues = {
  const random = seededUniformCh5(
    131041
    + 1031 * regenerateMix
    + sampleMix
    + Math.round(100 * mixWeight)
    + Math.round(100 * mixSeparation)
    + Math.round(100 * mixSigma2)
  );
  const values = new Array(sampleMix);
  const mu1 = -mixSeparation / 2;
  const mu2 = mixSeparation / 2;
  for (let i = 0; i < sampleMix; ++i) {
    if (random() < mixWeight) values[i] = mu1 + randnCauchy(random);
    else values[i] = mu2 + mixSigma2 * randnCauchy(random);
  }
  return values;
}

mixMean = mixValues.reduce((s, x) => s + x, 0) / mixValues.length
mixVar = mixValues.reduce((s, x) => s + (x - mixMean) ** 2, 0) / mixValues.length
mixSigma = Math.sqrt(mixVar)

histMix = histogramCauchy(mixValues, 100, -8, 8)
normalMix = Array.from({length: 401}, (_, i) => {
  const x = -8 + 16 * i / 400;
  return {x, y: Math.exp(-0.5 * (x - mixMean) ** 2 / mixVar) / (mixSigma * Math.sqrt(2 * Math.PI))};
})

mixPlot = Plot.plot({
  width: 760,
  height: 390,
  marginLeft: 55,
  x: {label: "x", domain: [-8, 8]},
  y: {label: "Плотность", domain: [0, 0.55]},
  grid: true,
  marks: [
    Plot.rectY(histMix, {x1: "x0", x2: "x1", y: "density", fill: "#4ea1ff", fillOpacity: 0.70, stroke: "white"}),
    Plot.line(normalMix, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 2.8})
  ]
})

html`
<div class="book-applet-panel">
  ${mixPlot}
  <div class="book-applet-summary">
    <span>Синий: смесь;</span>
    <span>оранжевый: гаусс с теми же средним и дисперсией;</span>
    <span>среднее = ${mixMean.toFixed(2)};</span>
    <span>σ = ${mixSigma.toFixed(2)}.</span>
  </div>
</div>
`

5.6.2 Пример: жидкий сцинтиллятор

Сигнал жидкосцинтилляционного детектора можно схематично записать как

\[ Q=Q(E,K)+Q_{\mathrm{stat}}+Q_{\mathrm{el}}, \tag{5.42}\]

где \(Q_{\mathrm{stat}}\) описывает фотонную статистику, \(Q_{\mathrm{el}}\) — шум электроники, а класс \(K\) включает тип частицы, топологию события, положение в детекторе, долю черенковского света и другие детали отклика.

Для \(\gamma\)-кванта энергия передаётся нескольким вторичным электронам:

\[ \gamma\longrightarrow e_1+e_2+\cdots+e_m, \qquad E_\gamma=\sum_{a=1}^{m}E_a. \tag{5.43}\]

Если световой выход электрона \(L(E)\) нелинеен, средний заряд имеет вид

\[ Q_\gamma\simeq\sum_{a=1}^{m}L(E_a). \tag{5.44}\]

При одном и том же \(E_\gamma\) разные наборы \(E_a\) дают разные значения \(Q_\gamma\). Разброс топологий добавляет к фотонной статистике второй член из формулы полной дисперсии 5.37. Описание сигнала одним пуассоновским или гауссовым распределением этот вклад потеряет.

5.7 Отношение случайных величин

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

Пусть \(X\) и \(Y\) независимы и

\[ X,Y\sim\mathcal N(0,1). \tag{5.45}\]

Для отношения \(R=X/Y\) совместная плотность после замены \(x=ry\) даёт

\[ \begin{aligned} f_R(r) &=\int_{-\infty}^{+\infty}|y|f_X(ry)f_Y(y)\,dy\\ &=\frac1{2\pi}\int_{-\infty}^{+\infty} |y|\exp\!\left[-\frac{(1+r^2)y^2}{2}\right]dy\\ &=\frac1{\pi(1+r^2)}. \end{aligned} \tag{5.46}\]

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

\[ A=\frac{N_1-N_2}{N_1+N_2}. \tag{5.47}\]

Линейное распространение ошибок для отношения применимо, когда знаменатель с высокой вероятностью далёк от нуля и его относительная неопределённость мала. Если знаменатель может оказаться около нуля, редкие значения отношения уходят далеко в хвост. Симметричная запись \(r\pm\sigma_r\) перестаёт описывать такое распределение.

viewof muDenom = Inputs.range([0.0, 5.0], {value: 0.0, step: 0.1, label: "Среднее знаменателя μ_Y"})
viewof sigmaDenom = Inputs.range([0.1, 2.0], {value: 1.0, step: 0.05, label: "Ширина знаменателя σ_Y"})
viewof sampleRatio = Inputs.range([1000, 40000], {value: 15000, step: 1000, label: "Число опытов N"})
viewof regenerateRatio = Inputs.button("Новая статистика", {
  value: 0,
  reduce: value => value + 1
})

ratioValues = {
  const random = seededUniformCh5(
    151051
    + 1033 * regenerateRatio
    + sampleRatio
    + Math.round(100 * muDenom)
    + Math.round(100 * sigmaDenom)
  );
  const values = [];
  for (let i = 0; i < sampleRatio; ++i) {
    const x = randnCauchy(random);
    const y = muDenom + sigmaDenom * randnCauchy(random);
    const r = x / y;
    if (Number.isFinite(r)) values.push(r);
  }
  return values;
}

histRatio = histogramCauchy(ratioValues, 120, -20, 20)
cauchyRatio = Array.from({length: 501}, (_, i) => {
  const x = -20 + 40 * i / 500;
  return {x, y: 1 / (Math.PI * (1 + x * x))};
})

ratioPlot = Plot.plot({
  width: 760,
  height: 390,
  marginLeft: 55,
  x: {label: "R = X/Y", domain: [-20, 20]},
  y: {label: "Плотность", domain: [0, 0.35]},
  grid: true,
  marks: [
    Plot.rectY(histRatio, {x1: "x0", x2: "x1", y: "density", fill: "#4ea1ff", fillOpacity: 0.70, stroke: "white"}),
    Plot.line(cauchyRatio, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 2.6})
  ]
})

html`
<div class="book-applet-panel">
  ${ratioPlot}
  <div class="book-applet-summary">
    <span>X ∼ N(0, 1);</span>
    <span>Y ∼ N(${muDenom.toFixed(1)}, ${sigmaDenom.toFixed(2)}²);</span>
    <span>оранжевый: распределение Коши для Y ∼ N(0, 1).</span>
  </div>
</div>
`

5.8 Причины отклонения от гауссова приближения

В таблице 5.1 собраны рассмотренные причины.

Таблица 5.1: Причины отклонения от простого гауссова приближения

Механизм	Какое условие требует проверки	Что видно в данных
Бесконечная дисперсия	Существует ли \(\operatorname{Var}(X_i)\)	Среднее может не сужаться
Тяжёлый хвост	Насколько быстро наступает предельный режим	Редкие большие значения сохраняются при конечном \(n\)
Доминирующий вклад	Мал ли каждый вклад по сравнению со всей суммой	Форму хвоста задаёт один механизм
Корреляции	Независимы ли слагаемые	Ширина не обязана убывать как \(1/\sqrt n\)
Смесь классов	Однородна ли выборка	Возникают дополнительная ширина, асимметрия или несколько пиков
Нелинейное преобразование	Является ли величина суммой	Отношение или асимметрия наследует границы и хвосты знаменателя

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

5.9 Итоги главы

• Конечная дисперсия и отсутствие доминирующего слагаемого входят в условия гауссова предела. У среднего величин Коши ширина не уменьшается с ростом выборки.
• Длинный правый хвост ионизационных потерь создают редкие \(\delta\)-электроны. В тонком слое используется ландау-подобное описание; кинематическая граница передачи энергии приводит к более точному распределению Вавилова.
• Положительные корреляции оставляют конечный вклад в дисперсию среднего. Общая систематическая ошибка может сохранить гауссову форму и одновременно нарушить закон \(1/\sqrt n\).
• Смесь классов получает дополнительную дисперсию от различия условных средних. Этот вклад описывается формулой полной дисперсии.
• Отношение двух независимых стандартных гауссовых величин имеет распределение Коши. Для нелинейной величины распределение нужно выводить отдельно.

5.10 Задачи

Задача 1. Среднее распределения Коши

Для независимых \(X_i\sim\mathrm{Cauchy}(0,1)\) покажите, что дисперсия не существует, а \[ \overline X_n=\frac{1}{n}\sum_iX_i \] имеет то же распределение Коши. Почему это не противоречит ЦПТ?

Задача 2. Тяжёлый хвост с конечной дисперсией

Пусть \[ p(x)\sim\frac{C}{|x|^{1+\alpha}}. \]

1. При каких \(\alpha\) существуют среднее и дисперсия?
2. Почему при \(\alpha\) немного больше 2 сходимость к гауссу может быть медленной?

Задача 3. Один доминирующий вклад

Пусть \[ S=X_0+\epsilon_1+\dots+\epsilon_n, \] где \(\epsilon_i\sim\mathrm{N}(0,\sigma^2)\), а \(X_0\) негауссов и \[ \mathrm{Var}(X_0)\gg n\sigma^2. \] Объясните форму распределения \(S\).

Задача 4. Общая систематическая ошибка

Пусть \[ X_i=\mu+C+\epsilon_i, \] где \(C\sim\mathrm{N}(0,\sigma_C^2)\) и \(\epsilon_i\sim\mathrm{N}(0,\sigma_\epsilon^2)\). Найдите \(\mathbb{E}[\overline X_n]\), \(\mathrm{Var}(\overline X_n)\) и предел дисперсии при \(n\to\infty\).

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

Пусть \(n\) величин имеют одинаковую дисперсию \(\sigma^2\) и одинаковый коэффициент корреляции \(\rho\) для любой пары.

1. Получите дисперсию их среднего.
2. Укажите допустимый диапазон \(\rho\) при фиксированном \(n\).
3. Рассмотрите случаи \(\rho=0\) и \(\rho>0\) при увеличении \(n\).

Задача 6. Смесь двух гауссов

Пусть \[ X|A\sim\mathrm{N}(\mu_A,\sigma_A^2),\qquad X|B\sim\mathrm{N}(\mu_B,\sigma_B^2), \] а вероятность класса \(A\) равна \(w\). Найдите полные среднее и дисперсию. Когда смесь плохо похожа на один гаусс?

Задача 7. Скрытая топология в жидком сцинтилляторе

Моноэнергетический \(\gamma\)-квант с энергией \[ E_\gamma=1\ \text{МэВ} \] может передать энергию вторичным электронам двумя способами: \[ \theta=A:\quad T_1=1\ \text{МэВ}, \] \[ \theta=B:\quad T_1=T_2=0.5\ \text{МэВ}. \]

Из-за нелинейности светового выхода средние сигналы различаются: \[ Q|\theta=A\sim\mathrm{Pois}(1000), \qquad Q|\theta=B\sim\mathrm{Pois}(980). \] Считайте классы равновероятными.

1. Найдите полное среднее \(\mathbb E[Q]\).
2. Вычислите среднюю внутриклассовую дисперсию \[ \mathbb E_\theta[\mathrm{Var}(Q|\theta)]. \]
3. Вычислите топологический вклад \[ \mathrm{Var}_\theta(\mathbb E[Q|\theta]). \]
4. Найдите полную дисперсию и стандартное отклонение \(Q\).
5. Сравните результат с пуассоновской моделью \[ Q\sim\mathrm{Pois}(\mathbb E[Q]). \] Какой вклад не описывается этой моделью?

Задача 8. Отношение гауссовых величин

Для независимых \(X,Y\sim\mathrm{N}(0,1)\) выведите плотность отношения \(R=X/Y\) и покажите, что это распределение Коши.

Литература

[1] Bichsel, H. Straggling in Thin Silicon Detectors. Reviews of Modern Physics, 60, 663–699. 1988. doi:10.1103/RevModPhys.60.663.

[2] Feller, W. An Introduction to Probability Theory and Its Applications. John Wiley & Sons. 1971.