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

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

R позволяет также гененировать случайные выборки из различных распределений заданного объема. Для этого используются функции типа rdist(), где вместо dist указывается необходимое распределение. Первый аргумент этой функции всегда объем выборки, а затем – параметры распределения. Доступные распределения можно найти здесь.

Бинарное распределение

Например, сгененируем выборку объема 8 из бинарного распределения с параметром \(p=\frac{1}{4}\). Для этого нужно попросить R сгенерировать выборку из биномиального распределения с параметрами \(n=1\) и \(p=\frac{1}{4}\) (так как бинарное распределение – это частный случай биномиального распределения).

rbinom(8, 1, 1/4)
## [1] 0 0 0 0 0 0 1 0

Каждый раз при реализации функции, генерирующей случайную выборку, у Вас будет появляться разная последовательность значений. Так и должно быть, ведь каждое из наблюдений выборки – случайная величина. В среднем, мы ожидаем, что четвертая часть выборки, то есть 2 наблюдения, должна содержать значения 1. Но при малельком размере выборки возможными также будут выборки, где 1 встречается чаще или реже.

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

Мы можем примерно оценить, как часто нам выпадет выборка, где содержится 1 единица. Это вероятность одного успеха в серии из 8 испытаний Бернулли с тем же самым параметром вероятности успеха:

dbinom(1, 8, 1/4)
## [1] 0.2669678
dbinom(0, 8, 1/4)
## [1] 0.1001129
dbinom(8, 8, 1/4)
## [1] 1.525879e-05

Вы видим, что примерно в 10 случаях из 100 выборка не будет содержать ни одной 1, а вероятность сплошных единиц очень мала – следовательно, такая выборка (одни единицы) неправдоподобна.

Примечание. Для того, чтобы воспроизвести сгенерированную случайную выборку при повторной реализации функции, нужно прописать перед ней функцию set.seed(any_number). Например:

set.seed(111)

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

Задание 3 из семинара 10

to be updated

Гистограммы в R

Гистограмма – график, который является аналогом графика функции плотности для распределний. Все наблюдения в выборке упорядочиваются по возрастанию. Затем упорядоченная выборка начиная от некоторого стартового значения (иногда это минимальное значение в выборке, а иногда – некоторое значение, меньшее, чем минимум) разбивается на \(n\) равных интервалов. После этого подсчитывается количество наблюдений, которые попали в каждый интервал. Это количество откладывется на оси Y. По оси X откладываются значения наблюдений в выборке.

Вызовем встроенную в R базу данных, где содержатся измерения температуры тела двух бобров в течение дня. После реализации команды data() в памяти R (окошко Environment) появятся две базы данных – по одной на бобра. Описание базы данных можно найти здесь.

data(beavers)

Посмотрим на базу данных. Выведем первые 6 строчек базы данных beaver1 – это делается с помощью функции head().

head(beaver1)

Нам нужна колонка (по-другому колонки в базе данных называются переменными) temp. Чтобы обратиться к колонке внутри базы данных, нужно сначала прописать название базы данных, затем поставить знак $, а затем прописать название переменной.

beaver1$temp
##   [1] 36.33 36.34 36.35 36.42 36.55 36.69 36.71 36.75 36.81 36.88 36.89
##  [12] 36.91 36.85 36.89 36.89 36.67 36.50 36.74 36.77 36.76 36.78 36.82
##  [23] 36.89 36.99 36.92 36.99 36.89 36.94 36.92 36.97 36.91 36.79 36.77
##  [34] 36.69 36.62 36.54 36.55 36.67 36.69 36.62 36.64 36.59 36.65 36.75
##  [45] 36.80 36.81 36.87 36.87 36.89 36.94 36.98 36.95 37.00 37.07 37.05
##  [56] 37.00 36.95 37.00 36.94 36.88 36.93 36.98 36.97 36.85 36.92 36.99
##  [67] 37.01 37.10 37.09 37.02 36.96 36.84 36.87 36.85 36.85 36.87 36.89
##  [78] 36.86 36.91 37.53 37.23 37.20 37.25 37.20 37.21 37.24 37.10 37.20
##  [89] 37.18 36.93 36.83 36.93 36.83 36.80 36.75 36.71 36.73 36.75 36.72
## [100] 36.76 36.70 36.82 36.88 36.94 36.79 36.78 36.80 36.82 36.84 36.86
## [111] 36.88 36.93 36.97 37.15

Чтобы построить гистограмму, нужна команда hist(). Аргументом является вектор: набор чисел, сохраненный в отдельный объект, или колонка из базы данных.

hist(beaver1$temp)

Настроим гистограмму.

hist(beaver1$temp,
     xlab = "Темература тела бобра", # изменить подпись оси OX
     ylab = "Количество наблюдений", # изменить подпись оси OY
     main = "Моя первая сливовая гистограмма про бобров", # изменить название гистограммы
     col = "plum" # поменять цвет гистограммы
     )

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

Мы можем наложить на гистограмму кривую нормального распределения с параметрами, которые соответствуют выборочному среднему и выборочному стандартному отклонению. Иными словами, мы можем сравнить, похожа ли гистограмма на нормальное распределение, которое мы ожидали бы увидеть с заданными параметрами?

Прежде всего поменяем шкалу оси OY, указав аргумент freq = F. Тогда по оси Y будут отложены значения плотности: так, чтобы площадь под гистограммой была равна 1.

hist(beaver1$temp,
     xlab = "Темература тела бобра", # изменить подпись оси OX
     ylab = "Плотность", # изменить подпись оси OY
     main = "Моя первая сливовая гистограмма про бобров", # изменить название гистограммы
     col = "plum", # поменять цвет гистограммы
     freq = F # плотность
     )

После этого рассчитаем выборочные средние и стандартное отклонение через команды mean() и sd() соответственно. Сохраним в новые объекты.

mean_beaver <- mean(beaver1$temp)
sd_beaver <- sd(beaver1$temp)

Теперь нарисуем на гистограмме график плотности распределения через функцию curve(dnorm()).

hist(beaver1$temp,
     xlab = "Темература тела бобра", # изменить подпись оси OX
     ylab = "Плотность", # изменить подпись оси OY
     main = "Моя первая сливовая гистограмма про бобров \nс кривой нормального распределения", # изменить название гистограммы
     col = "plum", # поменять цвет гистограммы
     freq = F # плотность
     )

curve(dnorm(x, # рассчитываем значения функции плотности нормального распределения
            mean = mean_beaver, # с выборочным средним
            sd = sd_beaver), # и выборочным стандартным отклонением
      add = T # дорисовываем кривую на предыдущем графике 
      )

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

Иллюстративно

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

Сгенерируем игрушечную “генеральную совокупность”: выборку из равномерного распределения, определенного на отрезке \([0,~10]\).

set.seed(545)
population <- runif(1000, min = 0, max = 10)
hist(population, col = "paleturquoise1", main = "Выборка из равномерного распределения")

Рассчитаем “истинное среднее”.

mean(population)
## [1] 5.004849

Возьмем 1000 выборок объемом 20, для каждой из них посчитаем выборочное среднее, и для этих 1000 средних нарисуем гистограмму.

set.seed(545)
sample20_means <- c()
for (i in 1:1000) {
  sample20_mean <- mean(sample(population, 20))
  sample20_means <- c(sample20_means, sample20_mean)
}

hist(sample20_means, col = "paleturquoise1", 
     main = "Среднее для 1000 подвыборок размером 20",
     xlab = "Выборочное среднее",
     xlim = c(2,8), # границы на оси OX
     ylim = c(0, 0.6),
     freq = F)

curve(dnorm(x, mean = mean(sample20_means), sd = sd(sample20_means)), add = T)

Несмотря на то, что игрушечная “генеральная совокупность” не распределена нормально, распределение выборочных средних из этой “генеральной совокупности” подчиняется нормальному распределению! Этот вывод дает нам центральная предельная теорема.

Посчитаем среднее и стандартное отклонение для этих 1000 средних:

mean(sample20_means)
## [1] 5.000925
sd(sample20_means)
## [1] 0.6730099

Закон больших чисел

Теперь возьмем 1000 выборок объемом 100, для каждой из них посчитаем выборочное среднее, и для этих 1000 средних нарисуем гистограмму. Границы графика оставим те же, что и для предыдущей гистограммы.

set.seed(545)
sample100_means <- c()
for (i in 1:1000) {
  sample100_mean <- mean(sample(population, 100))
  sample100_means <- c(sample100_means, sample100_mean)
}

hist(sample100_means, col = "paleturquoise1", 
     main = "Среднее для 1000 подвыборок размером 100",
     xlab = "Выборочное среднее",
     xlim = c(2,8) # границы на оси OX
     )

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

mean(sample100_means)
## [1] 4.994113
sd(sample100_means)
## [1] 0.2803005

Среднее стало (немного) ближе к истинному среднему, а стандартное отлоклоение – меньше. Это значит, что с увеличением объема выборки оценка среднего становится точнее! Этот вывод дает нам закон больших чисел.