В этой лекции мы замахнёмся на великое — на центральную предельную теорему в её классической форме. Мы напомним её формулировку, а затем проверим её в R экспериментально. Лекция частично основана на этой статье.

Немного теории

Итак, пусть \(x_1, x_2, \ldots, x_n\) — выборка, полученная в результате реализации некоторой случайной величины \(X\). Иными словами, \(X\) — это такой «генератор случайных чисел», с помощью которого мы сгенерировали \(n\) чисел. Можно думать, например, что есть большой ящик с бумажками, на каждой бумажке написано какое-то число, какие-то числа встречаются чаще, какие-то реже, какие-то вообще не встречаются. Когда задан такой ящик, это и означает, что задана случайная величина \(X\). Вытаскивая каждый раз наугад из этого ящика бумажку (и затем возвращая бумажку в ящик), мы будем получать последовательность чисел, которая и обозначается через \(x_1, x_2, \ldots, x_n\). Мы тщательно перемешиваем бумажки в ящике после того, как достаём очередную бумажку и результат каждого следующего вытаскивания не зависит от результата предыдущего. Сама по себе выборка тоже случайна (мы можем записать на бумажке \(n\) чисел, а потом запустить процедуру их генерирования ещё раз, и получить какие-то другие \(n\) чисел).

Предположим, что у случайной величины \(X\) есть математическое ожидание \(EX=m\) и стандартное отклонение (корень из дисперсии) \(\sigma\).

Рассмотрим выборочное среднее \(\bar x=\frac{1}{n}(x_1+x_2+\ldots+x_n)\). Это тоже случайная величина, она каждый раз меняется, когда мы создаём новую выборку, в отличие от \(m\), которая не меняется.

Закон больших чисел учит нас, что если \(n\) достаточно большое, то \(\bar x\) становится близко к \(m\). Но, конечно, для каждой конкретной выборки \(\bar x\) может отклоняться от \(m\). Как именно сильно оно отклоняется? Можно сказать, что чем больше \(n\), тем меньше отклонение, но чтобы ответить на этот вопрос более точно, нужно составить величину, обозначаемую обычно буквой \(Z\):

\[Z=\frac{\bar x-m}{\sigma/\sqrt{n}}=\sqrt{n}\frac{\bar x - m}{\sigma}\]

Как видно, она меряет, в какой мере \(\bar x\) отклоняется от \(m\), компенсируя тот факт, что это отклонение уменьшается при увеличении \(n\) (умножение на \(\sqrt{n}\)). Центральная предельная теорема гласит, что \(Z\) распределена «примерно» по стандартному нормальному закону, причём чем больше \(n\), тем лучше работает это «примерно».

Проверяем в R

Наш план состоит в том, чтобы сгенерировать много-много выборок из какого-нибудь распределения, посчитать для каждой из них значение \(\bar x\) и \(Z\), записать их в списки (каждую в свой) и построить гистограммы того, что получилось. Для этого нам придётся немножко познакомиться с программированием в R.

Цикл for

Цикл for в R устроен примерно так же, как в Python — он служит для перебора элементов списка.

lst=c(2,6,1)
for (x in lst) {
  print(x)
}
## [1] 2
## [1] 6
## [1] 1

В отличие от Python, тело цикла должно быть заключено в фигурные скобки, но в остальном они работают очень похоже. Самое простое применение циклов — запустить фрагмент кода заданное число раз, делается так:

for (i in 1:5) {
  print("Hello, R")
}
## [1] "Hello, R"
## [1] "Hello, R"
## [1] "Hello, R"
## [1] "Hello, R"
## [1] "Hello, R"

Если бы мы хотели создать список квадратов натуральных чисел от 1 до n с помощью циклов, то это можно было сделать так.

n <- 5
squares <- rep(0,n)

Здесь мы инициализировали список squares нулями, чтобы было, куда записывать результат. (Хотя можно было бы вместо этого создать пустой список и добавять к нему элементы с помощью append.) Функция rep повторяет заданный в первом аргументе элемент столько раз, сколько сказано во втором аргументе.

for(i in 1:n) {
  squares[i] <- i**2
}
squares
## [1]  1  4  9 16 25

Конечно, это не самый эффективный способ создать список из квадратов в R, гораздо проще и эффективнее было бы написать (1:n)**2, но зато мы показали, как работают циклы, и теперь готовы к решению нашей задачи.

Постановка задачи

В качестве «подопытной» случайной величины мы возьмём равномерно распределенную на отрезке \([0;1]\) – довольно-таки непохожую на нормальную.

Упражнение. Чему равно математическое ожидание и стандартное отклонение такой случайной величины?

Ответ: \(EX=m=1/2\), \(DX=\sigma^2=\frac{1}{12}\), \(\sigma=\sqrt{DX}=\frac{1}{\sqrt{12}}\).

А теперь собственно задача для R (обсуждалась на семинаре):

  1. С помощью rep создать списки xbars и zs длиной \(k\), инициализированные нулями.
  2. Построить выборку из n случайных величин, равномерно распределенных на отрезке \([0,1]\) (с помощью runif).
  3. Найти выборочное среднее для этой выборки.
  4. Найти значение \(Z\)-статистики (\(\frac{\bar x-m}{\sigma/\sqrt{n}}\)) для этой выборки.
  5. Записать полученное выборочное среднее в xbars, а \(Z\) в zs.
  6. Операции 2–5 повторять \(k\) раз.
  7. Построить гистограммы для xbars и zs с параметром freq=FALSE. На картинке с zs построить также кривую плостности стандартного нормального распределения.

Вариант решения

k <- 1000
n <- 200
xbars <- rep(0,k)
zs <- rep(0,k)
m <- 0.5
sigma=sqrt(1/12)
for (i in 1:k)
{
  s <- runif(n)
  xbar <- mean(s)
  xbars[i] <- xbar
  zs[i] <- (xbar-m)/sigma*sqrt(n)
}
hist(xbars, freq=FALSE)

hist(zs, freq=FALSE)
curve(dnorm(x), add=TRUE)

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

Контрольный вопрос: чем отличаются гистограммы для xbars и zs? Почему?

Задание на дом

В задачах статистики часто возникает ситуация, когда математическое ожидание \(m\) в формуле для \(Z\) можно считать известным, а стандартное отклонение \(\sigma\) — нет. В этом случае стандартное отклонение может быть оценено как квадратный корень из исправленной выборочной дисперсии:

\[\sigma_+^2=\frac{1}{n-1}\sum_{i=1}^n (x_i-\bar x)^2\]

Если выборка достаточно большая (обычно считается, что \(n>35\) – это большая), то в формулу для \(Z\) можно подставить \(\sigma_+\) вместо настоящего \(\sigma\) и дальше пользоваться нормальным законом как обычно. А если выборка маленькая, то нормальное распределение «перестаёт работат»: вместо кривой нормального распределения придётся использовать другую кривую (похожую), называемую кривой \(t\)-распределения или распределения Стьюдента.

Давайте посмотрим, как это выглядит на практике. Для этого нужно в приведенный выше код внести изменения:

  1. В формуле для вычисления z подставить вместо заранее найденного sigma исправленное выборочное стандартное отклонение для нашей выборки s (вычисляется с помощью функции sd).
  2. Установить число элементов выборки n в какое-нибудь небольшое целое число (например, n <- 2 или n <- 5). Построить гистограмму для zs с такими параметрами.
  3. Скорее всего в списке zs теперь встречаются очень большие или очень маленькие значения. Проверить это с помощью команды summary. Откуда берутся эти значения?
  4. Эти экстремально большие или маленькие значения приводят к тому, что на гистограмме практически ничего не видно. Это можно исправить, изменив минимальное и максимальное значение на горизонтальной оси с помощью параметров xlim и breaks у функции hist. Например, добавьте xlim = c(-10,10), breaks=100 при вызове функции hist и посмотрите, что получается. «Поиграйте» со значением breaks, чтобы получить разумно выглядящую картинку.
  5. Добавьте к гистограмме zs графики функции плотности стандартного нормального распределения (собственно, это уже сделано в исходной программе) и распределения Стьюдента (через функцию dt). Второй параметр функции dt — так называемое число степеней свободы, оно в данном случае должно быть на единицу меньше, чем \(n\) (число элементов выборки). Чтобы различать кривые, можно добавить параметр col='red' или col='green', устанавливающие соответствующий цвет. Какая из кривых лучше приближает распределение данных в zs? На каких участках? Поиграйте с различными значениями \(n\).