Центральная предельная теорема и нормальное распределение

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

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

План лекции

• Псевдослучайные числа
• Суммы случайных величин
• Гаусс, квантили и интервалы
• Совместимость измерений

Псевдослучайные числа

Генератор псевдослучайных чисел

• В конце первой лекции мы подбирали \(\sigma\) для лучшего описания псевдоданных.
• Но как мы получили псевдоданные?
• Для генерации псевдослучайных чисел используются алгоритмы, которые выдают последовательность чисел, статистически похожих на случайные.
• Такие алгоритмы называются генераторами псевдослучайных чисел.
• Достаточно иметь генератор равномерного распределения \(X\sim U(-1,1)\), чтобы построить генераторы для других распределений.

Генерация и визуализация псевдослучайных чисел

• Код
• Данные
• График

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(7)
x = rng.uniform(-1, 1, size=20_000)
preview = pd.DataFrame({"i": np.arange(1, 101), "x": x[:100]})

• default_rng(7) создаёт воспроизводимый генератор.
• uniform(-1,1,size=20_000) даёт выборку из равномерного распределения.
• preview сохраняет первые 100 значений для просмотра.

Первые 100 значений выборки:

i	x
1	0.250
2	0.794
3	0.551
4	-0.550
5	-0.400
6	0.747
7	-0.989
8	0.642
9	0.594
10	-0.064
11	-0.394
12	-0.443
13	-0.490
14	-0.110
15	0.009
16	0.107
17	0.991
18	0.585
19	0.244
20	0.978
21	-0.569
22	-0.680
23	0.225
24	-0.912
25	-0.929
26	0.030
27	-0.068
28	0.834
29	0.258
30	0.028
31	-0.006
32	-0.505
33	-0.976
34	-0.615
35	0.384
36	-0.599
37	-0.261
38	-0.993
39	0.660
40	-0.691
41	-0.465
42	0.761
43	0.020
44	0.694
45	0.279
46	0.484
47	-0.817
48	0.082
49	0.016
50	0.743
51	-0.277
52	0.196
53	-0.881
54	-0.225
55	-0.354
56	-0.700
57	0.633
58	-0.241
59	0.957
60	0.180
61	0.210
62	0.276
63	0.353
64	-0.698
65	-0.119
66	-0.521
67	-0.195
68	-0.807
69	0.936
70	-0.570
71	0.344
72	-0.399
73	0.748
74	0.324
75	-0.737
76	0.690
77	0.890
78	0.808
79	0.139
80	-0.709
81	-0.615
82	0.856
83	0.105
84	-0.639
85	0.768
86	0.283
87	0.139
88	-0.247
89	-0.178
90	-0.521
91	-0.924
92	0.752
93	-0.065
94	0.095
95	-0.356
96	0.503
97	-0.950
98	-0.256
99	-0.939
100	-0.754

Суммы случайных величин

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

• Рассмотрим случайную величину \[ S_m=X_1+\dots+X_m. \]

• Вопрос: каким становится распределение \(S_m\), если \(m\to \infty\)?

• Для примера возьмём независимые \(X_i\sim U(-1,1)\).

• У каждой величины: \[ \begin{aligned} \mathbb{E}[X_i]&=0,\\ \mathrm{Var}(X_i)&=1/3. \end{aligned} \]

• Короткий вывод: \[ \mathbb{E}[X_i^2]=\int_{-1}^{1} x^2\cdot \frac12\,dx=\frac13, \qquad \mathrm{Var}(X_i)=\mathbb{E}[X_i^2]-\mathbb{E}[X_i]^2=\frac13. \]

Среднее суммы

• Определение среднего: \[ \mathbb{E}[S_m] =\mathbb{E}\!\left[\sum_{i=1}^m X_i\right]=\sum_{i=1}^m\mathbb{E}[X_i]. \]
• Если \(X_i\sim U(-1,1)\), то \[ \boxed{\mathbb{E}[S_m]=0.} \]

Дисперсия суммы

По определению, \[ \mathrm{Var}(S_m)=\mathbb{E}\!\left[(S_m-\mathbb{E}[S_m])^2\right]. \]

\[ \mathrm{Var}(S_m) =\mathbb{E}\!\left[\left(\sum_{i=1}^m (X_i-\mathbb{E}[X_i])\right)^2\right]. \]

\[ \mathrm{Var}(S_m) =\sum_{i=1}^m \mathbb{E}\!\left[(X_i-\mathbb{E}[X_i])^2\right] +2\sum_{i<j}\mathbb{E}\!\left[(X_i-\mathbb{E}[X_i])(X_j-\mathbb{E}[X_j])\right]. \]

Смешанные члены подсказывают новое понятие: ковариацию \[ \boxed{ \mathrm{cov}(X,Y)=\mathbb{E}\!\left[(X-\mathbb{E}[X])(Y-\mathbb{E}[Y])\right] =\mathbb{E}[XY]-\mathbb{E}[X]\mathbb{E}[Y].} \] В частности, \(\mathrm{Var}(X)=\mathrm{cov}(X,X)\).

Завершение вычисления

Теперь формула принимает вид \[ \mathrm{Var}(S_m) =\sum_{i=1}^m\mathrm{Var}(X_i) +2\sum_{i<j}\mathrm{cov}(X_i,X_j). \]

Для независимых величин ковариации равны нулю: \[ \mathrm{Var}(S_m)=\sum_{i=1}^m\mathrm{Var}(X_i). \]

Если \(X_i\sim U(-1,1)\), то \[ \boxed{\mathrm{Var}(S_m)=m/3.} \]

Формулировка ЦПТ

Теорема

Если \(X_1,\dots,X_m\) независимы и имеют конечные дисперсии, то нормированная сумма \[ Z_m=\frac{S_m-\mathbb{E}[S_m]}{\sqrt{\mathrm{Var}(S_m)}} \] при большом \(m\) стремится к нормальному распределению \(\mathrm{N}(0,1)\).

ЦПТ для равномерных величин

Для \(X_i\sim U(-1,1)\): \[ Z_m=\frac{X_1+\dots+X_m}{\sqrt{m/3}}. \]

• При \(m=1\) распределение ещё равномерное.
• При росте \(m\) форма быстро становится похожей на гауссову.

ЦПТ: \(m=1\)

ЦПТ: \(m=2\)

ЦПТ: \(m=5\)

ЦПТ: \(m=10\)

Интерактив: ЦПТ

viewof m = Inputs.range([1, 50], {
  value: 5,
  step: 1,
  label: "m"
})

viewof sampleSize = Inputs.range([100, 30000], {
  value: 12000,
  step: 1000,
  label: "N"
})

viewof regenerate = Inputs.button("Новая статистика")

zValues = {
  regenerate;
  const scale = Math.sqrt(m / 3);
  const values = new Array(sampleSize);

  for (let i = 0; i < sampleSize; ++i) {
    let sum = 0;
    for (let j = 0; j < m; ++j) {
      sum += 2 * Math.random() - 1;
    }
    values[i] = sum / scale;
  }

  return values;
}

histogram = (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)
  }));
}

hist = histogram(zValues, 64, -4, 4)

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

zMean = zValues.reduce((sum, value) => sum + value, 0) / zValues.length

zVariance = zValues.reduce(
  (sum, value) => sum + (value - zMean) ** 2,
  0
) / zValues.length

cltPlot = Plot.plot({
  width: 960,
  height: 580,
  marginLeft: 55,
  x: {label: "z", domain: [-4, 4]},
  y: {label: "density", domain: [0, 0.55]},
  grid: true,
  marks: [
    Plot.rectY(hist, {
      x1: "x0",
      x2: "x1",
      y: "density",
      fill: "#4ea1ff",
      fillOpacity: 0.70,
      stroke: "white"
    }),
    Plot.line(normalCurve, {
      x: "z",
      y: "density",
      stroke: "#ffb347",
      strokeWidth: 2.8
    })
  ]
})

html`
<div class="clt-panel">
  ${cltPlot}
  <div class="clt-summary">
    <span>z = (X1 + ... + Xm)/sqrt(m/3)</span>
    <span>m = ${m}</span>
    <span>mean = ${zMean.toFixed(3)}</span>
    <span>variance = ${zVariance.toFixed(3)}</span>
  </div>
</div>
`

Гаусс, квантили и интервалы

Почему гаусс так важен

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

В этой лекции: одномерный гаусс

• Как устроена плотность \(\mathrm{N}(\mu,\sigma)\).
• Как читать вероятности, квантили и интервалы в единицах \(\sigma\).
• Как сравнивать и усреднять независимые измерения.
• Как переносить ошибку через функцию.
• Многомерный гаусс и корреляции подробно разберём в другой лекции.

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

\[ f(x;\mu,\sigma)= \frac{1}{\sigma\sqrt{2\pi}} \exp\!\left[ -\frac{(x-\mu)^2}{2\sigma^2} \right] \]

• \(\mu\) задаёт наиболее вероятное значение распределения \(x\)
• \(\sigma\) задаёт характерную ширину распределения
• Дисперсия равна \(\sigma^2\)
• Нормировка: \(\int_{-\infty}^{+\infty} f(x;\mu,\sigma)\,dx = 1\)

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

Если \[ X\sim \mathrm{N}(\mu,\sigma), \] то после замены \[ Z=\frac{X-\mu}{\sigma} \] получаем стандартное нормальное распределение \[ \boxed{ Z\sim \mathrm{N}(0,1), \qquad f_Z(z)=\frac{1}{\sqrt{2\pi}}e^{-z^2/2}.} \]

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

• В лекции о случайных величинах и распределениях мы ввели кумулятивную функцию распределения: \[ \begin{aligned} \Phi(z)=P(Z\le z)&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{z}e^{-t^2/2}\,dt =\\ &=\frac12\left[1+\operatorname{erf}\!\left(\frac{z}{\sqrt2}\right)\right]. \end{aligned} \]

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

• Квантиль уровня \(\beta\) для \(Z\sim \mathrm{N}(0,1)\) — это число \(z_\beta\), такое что
  \[ \Phi(z_\beta)=\beta. \]

• Смысл: \[ P(-\infty < Z \le z_\beta) = \beta \]

• Для центральной области вероятности \(\gamma\): \[ P\!\left(|Z|\le z_{(1+\gamma)/2}\right)=\gamma. \] Если задан уровень значимости \(\alpha=1-\gamma\), то та же граница равна \(z_{1-\alpha/2}\).

Уровень \(\beta\)	\(z_\beta\)
0.90	1.282
0.95	1.645
0.975	1.960
0.995	2.576

Центральные интервалы

• Код
• График

from statistics import NormalDist
import numpy as np

dist = NormalDist()
z = np.linspace(-4, 4, 800)
f = np.exp(-z**2 / 2) / np.sqrt(2 * np.pi)
p1 = dist.cdf(1) - dist.cdf(-1)
p2 = dist.cdf(2) - dist.cdf(-2)
p3 = dist.cdf(3) - dist.cdf(-3)

• p1, p2, p3 дают вероятности попадания в интервалы [-1,1], [-2,2], [-3,3].
• Для стандартной нормали это и есть центральные интервалы 1σ, 2σ, 3σ.
• На графике они сопоставлены с шириной гауссовой кривой для 1σ, 2σ, 3σ.

Совместимость измерений

Совместимость двух измерений

• Пусть \(x_1\pm \sigma_1\), \(x_2\pm \sigma_2\) независимо измеряют одну и ту же величину \(\mu\).

• Как проверить, что эти измерения совместимы, то есть не противоречат друг другу?

• Например, какие из этих пар измерений совместимы, а какие нет?

  • первая пара (оценки сильно отличаются, но ошибки велики): \[ 1.0\pm0.5, \qquad 2.0\pm0.5. \]
  • вторая пара (оценки близки, но ошибки малы): \[ 1.010\pm0.001, \qquad 1.020\pm0.001. \]

Распределение разности

• Проверяем не сами значения по отдельности, а разность: \[ D=X_1-X_2. \]

• Если оба измерения оценивают одну величину \(\mu\), то при гипотезе совместимости \[ \mathbb{E}[D]=\mathbb{E}[X_1-X_2]=\mu-\mu=0. \]

• Дисперсия разности: \[ \mathrm{Var}(D)=\mathrm{Var}(X_1-X_2)=\sigma_1^2+\sigma_2^2-2\,\mathrm{cov}(X_1,X_2). \]

Распределение разности

• Для независимых измерений \[ \mathrm{cov}(X_1,X_2)=0, \qquad \sigma_D^2=\sigma_1^2+\sigma_2^2. \]

• Если ошибки гауссовы, то \[ D\sim \mathrm{N}(0,\sigma_D), \qquad \sigma_D=\sqrt{\sigma_1^2+\sigma_2^2}. \]

• Поэтому \(Z=D/\sigma_D\) имеет стандартное нормальное распределение: \[ Z=\frac{D}{\sigma_D} =\frac{D}{\sqrt{\sigma_1^2+\sigma_2^2}} \sim \mathrm{N}(0,1). \]

• Для наблюдённых значений \(x_1\), \(x_2\) получаем \[ \boxed{ z_{\rm obs} = \frac{x_1-x_2}{\sqrt{\sigma_1^2+\sigma_2^2}}. } \]

Алгоритм проверки совместимости двух измерений

• Мы свели задачу к стандартной нормальной величине \(Z\sim\mathrm{N}(0,1)\).

• Теперь спрашиваем: попало ли наблюдаемое значение \(z_{\rm obs}\) в центральную область стандартного гаусса?

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

• Если попало в один из хвостов, такая разность слишком редкая для гипотезы совместимости.

• Для данных считаем \[ \boxed{ z_{\rm obs} = \frac{x_1-x_2}{\sqrt{\sigma_1^2+\sigma_2^2}} } \]

Алгоритм проверки совместимости двух измерений

• Выбираем центральную область совместимости: \[ \gamma=0.90, \qquad \alpha=1-\gamma=0.10. \]

• Критическая граница: \[ c=z_{(1+\gamma)/2}=z_{1-\alpha/2}. \]

• Для \(\gamma=0.90\): \[ c=z_{0.95}=1.64. \]

• Правило: \[ |z_{\rm obs}|\le c \] совместимость не отвергается; при \(|z_{\rm obs}|>c\) измерения несовместимы на уровне значимости \(\alpha\).

Совместимость: важна не разность, а значимость

• Пример A: оценки сильно отличаются, но ошибки велики. \[ 1.0\pm0.5, \qquad 2.0\pm0.5 \] \[ z_{\rm obs} = \frac{1.0-2.0}{\sqrt{0.5^2+0.5^2}} = -1.41 \]

• Пример B: оценки близки, но ошибки малы. \[ 1.010\pm0.001, \qquad 1.020\pm0.001 \] \[ z_{\rm obs} = \frac{1.010-1.020}{\sqrt{0.001^2+0.001^2}} = -7.07 \]

• Дальше двигаем \(z_{\rm obs}\) в апплете и смотрим, попадает ли он в центральную область.