Введение

Данный короткий обзор имеет целью представить способы использования языка R для решения некоторых несложных задач, а именно:

  1. Моделирование заданного распределения и его гисторамма. Примеры непрерывного распределения, дискретного распределения и гистограмма после стандартизации дискретного распределения.
  2. Как нарисовать на одном графике плотность распределения и гистограмму по выборке.
  3. Как нарисовать на одном графике функцию распределения и график накопленных частот.

Моделирование различных распределений

В данном разделе рассмотрим примеры непрерывных и дискретных распределений, промоделируем некоторые из них, и построим соответствующие гистограммы.

Непрерывные распределения

Нормальное (Гауссовское) распределение: моделирование, гистограмма, функция плотности

Предположим, что перед нами встала задача моделирования одномерного нормального распределения \(\mathrm N(0,1)\). В данном случае можно пользоваться семейством следующих функций:

  1. rnorm() — для моделирования нормальных случайных величин;

  2. dnorm() — для вычисления функции плотности нормального распределения в точке (в частности, для построения графика плотности);

  3. pnorm() — для вычисления функции распределения \(\mathrm N(a,\sigma^2)\) в точке (и для построения её графика);

Напомним, что для уточнения работы той или иной функции можно использовать справку прямо в консоли R; например:

? rnorm

Итак, промоделируем \(100\) псевдослучайных величин из нормального распределения со стандартными параметрами с помощью функции rnorm().

set.seed(1) #установка начала псевдослучайной последовательности, нужно делать это в самом начале программы для воспроизводимости вычислений
size_of_sample <-  100
mean <-  0
sd <- 1
sample <- rnorm(size_of_sample, mean, sd)

Теперь, с помощью функции head() посмотрим на несколько первых элементов выборки.

head(sample)
## [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684

Несложно построить гистограмму по значениям из sample.

hist(sample, main = "Гистограмма выборки из N(0,1)", ylab = "Частота", col = "light blue")

После построения гистограммы логично на том же графике построить график плотности. Для этого будем использовать функцию lines() для построения кривой на уже прорисованном графике с гистограммой. Однако, чтобы на оси \(y\) гистограмма и плотность отображались в корректном масштабе, в функции hist() теперь следует задать параметр freq = FALSE.

hist(sample, main = "Гистограмма выборки из N(0,1) и плотность", ylab = "Плотность", col = "light blue", freq = FALSE)
x <- seq(-3, 3, length.out = 100) #строим последовательность из 100 чисел, заполняющих промежуток от -3 до 3 как значения аргумента для подстановки в функцию плотности
y <- dnorm(x, mean, sd) #подставляем x в функцию плотности и получаем вектор
lines(x, y, lwd = 2, col = "orange") #строим кривую поверх гистограммы

Построим теперь на одном и том же графике функцию распределения и график накопленных частот (эмпирической функции распределения) для sample. Как и с гистограммой, сначала вызываем функцию, которая строится на нашей выборке sample — plot.ecdf(), а затем уже поверх полученного строим теоретическую функцию распределения.

plot.ecdf(sample, col = "orange", cex = 0.7) #параметр cex отвечает за кратное увеличение/уменьшение размера точки на оранжевом графике
x <- seq(-3, 3, length.out = 100)
y <- pnorm(x, mean, sd)
lines(x, y, lwd = 3, col = "light blue")

Разумно сейчас посмотреть на полученную гистограмму, а также на эмпирическую функцию распределения в сравнении с плотностью и теоретической функцией распределения, и усомниться в том, что всё хорошо: визуально они не совсем похожи. Давайте увеличим размер выборки до \(10000\) и посмотрим, как будут выглядеть те же графики. Сначала гистограмма с плотностью:

set.seed(1) #пусть выборка начинается с тех же чисел, что и было раньше, для этого установим сид
size_of_sample <- 10000 #изменили размер выборки, дальше всё то же самое
mean <-  0
sd <-  1
sample <- rnorm(size_of_sample, mean, sd) 

hist(sample, main = "Гистограмма выборки из N(0,1) и плотность", ylab = "Плотность", col = "light blue", freq = FALSE)
x <- seq(-3, 3, length.out = 100) #строим последовательность из 100 чисел, заполняющих промежуток от -3 до 3 как значения аргумента для подстановки в функцию плотности
y <- dnorm(x, mean, sd) #подставляем x в функцию плотности и получаем вектор
lines(x, y, lwd = 2, col = "orange") #строим кривую поверх гистограммы

Выглядит уже весьма правдоподобно, если обратить внимание на симметричность и убывание хвоста.

Эмпирическая функция распределения вместе с теоретической:

plot.ecdf(sample, col = "orange", cex = 0.7) 
x <- seq(-3, 3, length.out = 100)
y <- pnorm(x, mean, sd)
lines(x, y, lwd = 3, col = "light blue")

Всё отлично, кажется тут можно обойтись без комментариев.

Другие непрерывные распределения

Несложно найти функции, действующие аналогично rnorm(), dnorm(), pnorm(), и для многих других непрерывных распределений в базовом пакете stats. Приведем несколько примеров:

  • Экспоненциальное распределение (rexp(), dexp(), pexp());
  • Распределение Коши (rcauchy(), dcauchy(), pcauchy());
  • Гамма распределение (rgamma(), dgamma(), pgamma()).

Обычный поиск в RStudio скорее всего приведёт к нужным функциям. Если нет, то в CRAN обычно можно найти пакеты, которые содержат такие же функции для специфических распределений; установка этих пакетов в RStudio тоже проста.

Дискретные распределения

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

Пример с распределением Пуассона

Для начала возьмём распределение Пуассона и просто построим гистограмму на основании выборки размера \(10000\). Так как выборка берётся из дискретного распределения, самостоятельно выбираем расположение столбцов гистограммы с помощью аргумента breaks, а также right = FALSE, так как нужно, чтобы первый столбец отвечал за \(0\), второй за \(1\) и т.д.

set.seed(1)
size_of_sample <- 10000
lambda <- 1
sam <- rpois(size_of_sample, lambda) #генерируем выборку
hist_pois <- hist(sam, breaks = 0:(max(sam)+1), freq = F, right = FALSE, col = "light blue", xlim = c(0,(max(sam)))) #рисуем гистограмму и сохраняем объект гистограмма в переменную hist_pois

hist_pois$density
## [1] 0.3700 0.3599 0.1862 0.0653 0.0161 0.0022 0.0002 0.0001

Существует много различных способов изобразить две гистограммы на одном рисунке. Мы воспользуемся методом, который не требует знания пакета ggplot, а использует базовую функцию barplot(), которая строит столбчатые диаграммы. Всё это делается в несколько этапов:

  1. Сначала строим столбчатую диаграмму на основании теоретического распределения, используя значения функции вероятности dpois(). У распределения Пуассона вероятности убывают очень быстро, поэтому не будем рисовать столбцы большие максимального значения в нашей выборке. Из нового: параметры density и angle задают соответственно плотность и угол поворота прямых, которые позволяют отличать столбцы при их совместном отображении. Параметр names.arg задаёт имена каждого столбца и позволяет отметить числа на оси \(x\)l;

  2. Далее строим эмпирическую гистограмму поверх старой, задав аргумент add = TRUE. В качестве основы берём объект гистограммы hist_pois, который мы создали выше (“забираем” у неё вычисленные плотности посредством обращения hist_pois$density). Меняем цвет, плотность и поворот прямых у столбцов;

  3. Указываем название посредством функции title();

  4. Добавляем легенду для графика с помощью функции legend().

barplot(dpois(0:max(sam), lambda), col = "blue", border = "blue", density = 10, angle = 45, xlab = "Значение", ylab = "Частота", names.arg = c(0:max(sam)))
barplot(hist_pois$density, col = "orange", border = "orange", density = 10, angle = -45, add = TRUE)
title(main = "Сравнение теоретических вероятностей и эмпирических частот")
legend("topright", legend = c("Теоретические вероятности", "Эмпирические частоты"), angle = c(45, -45), density = c(20,20), border = c("blue", "orange"))

Для того, чтобы сейчас построить график накопленных частот, главное знать о существовании функции cumsum(), которая, получая на вход вектор из вероятностей, на выходе выдаёт вектор накопленных вероятностей. Покажем на примере, чтобы точно было понятно:

dpois(0:max(sam), lambda)
## [1] 3.678794e-01 3.678794e-01 1.839397e-01 6.131324e-02 1.532831e-02
## [6] 3.065662e-03 5.109437e-04 7.299195e-05
cumsum(dpois(0:max(sam), lambda))
## [1] 0.3678794 0.7357589 0.9196986 0.9810118 0.9963402 0.9994058 0.9999168
## [8] 0.9999898

Поэтому построение графиков накопленных частот практически ничем не отличается от ситуации выше. Вместо barplot() используем функции plot() и lines, меняем теоретические вероятности на накопленные посредством cumsum(), то же делаем и с эмпирическими частотами. Получаем следующий код и график:

plot(0:max(sam), cumsum(dpois(0:max(sam), lambda)), col = "light blue", xlab = "Значение", ylab = "Накопленная частота", type = "s", lwd = 4, lty = 2)
lines(0:max(sam), cumsum(hist_pois$density), col = "orange", lwd = 2, type = "s")
title(main = "Сравнение накопленных вероятностей и частот")
legend("bottomright", legend = c("Накопленные вероятности", "Накопленные частоты"), col = c("light blue", "orange"), lwd = c(4, 2), lty = c(2, 1))

Пример с геометрическим распределением

Проделаем похожие манипуляции с геометрическим распределением и построим гистограммы и графики сравнения частот/вероятностей и накопленных частот/вероятностей. В этот раз используем функции rgeom() и dgeom().

set.seed(1)
size_of_sample <- 10000
prob <- 0.3
sam <- rgeom(size_of_sample, prob) #генерируем выборку
hist_geom <- hist(sam, breaks = 0:(max(sam)+1), freq = F, right = FALSE, col = "light blue", xlim = c(0,(max(sam)))) #рисуем гистограмму и сохраняем объект гистограмма в переменную hist_geom

Опять же, сравнение теоретических вероятностей и эмпирических частот:

barplot(dgeom(0:max(sam), prob), col = "blue", border = "blue", density = 15, angle = 45, xlab = "Значение", ylab = "Частота", names.arg = c(0:max(sam)))
barplot(hist_geom$density, col = "orange", border = "orange", density = 15, angle = -45, add = TRUE)
title(main = "Сравнение теоретических вероятностей и эмпирических частот")
legend("topright", legend = c("Теоретические вероятности", "Эмпирические частоты"), angle = c(45, -45), density = c(20,20), border = c("blue", "orange"))

А также сравенение накопленных теоретических вероятностей и эмпирических частот:

plot(0:max(sam), cumsum(dgeom(0:max(sam), prob)), col = "light blue", xlab = "Значение", ylab = "Накопленная частота", type = "s", lwd = 3, lty = 2)
lines(0:max(sam), cumsum(hist_geom$density), col = "orange", lwd = 2, type = "s")
title(main = "Сравнение накопленных вероятностей и частот")
legend("bottomright", legend = c("Накопленные вероятности", "Накопленные частоты"), col = c("light blue", "orange"), lwd = c(3, 2), lty = c(2, 1))

Стандартизация

Также иногда требуется стандартизировать дискретное распределение (например, для того, чтобы более осмысленно сравнивать разные по природе выборки или, например, “похожесть” дискретного распределения на нормальное). Это делается несложным преобразованием: вычитанием среднего и делением на стандартное отклонение \[\frac{X_i - \bar{x}}{s},\] где \(X_i\) — элемент из выборки \(\{X_i\}_{i = 1}^{n}\), \(\bar{x}\) — среднее, а \(s\) — стандартное отклонение.

Промоделируем выборку из Пуассоновского распределения с достаточно большим параметром \(\lambda\), например пусть \(\lambda = 100\).

set.seed(10000)
lambda <- 100
size_of_sample <- 10000
sample <- rpois(size_of_sample, lambda)
hist_pois2 <- hist(sample, freq = FALSE, right = FALSE, col = "light blue")

Допустим, мы хотим сравнить это дискретное распределение со стандартным нормальным распределением, хотя бы визуально. Для этого можно стандартизировать нашу выборку и построить гистограмму вместе с плотностью нормального распределения. (!) Заметим, что вместе со стандартизацией выборки мы также применим то же преобразование и к вектору hist_pois2$breaks, чтобы сохранить ширину и позицию столбцов гистограммы при стандартизации. Если этого не сделать, то столбцы гистограммы после стандартизации могут существенно изменить свою форму, чего происходить не должно, так как кроме масштабирования и сдвига мы ничего не сделали.

stand_sample <- (sample - mean(sample))/sd(sample)
stand_breaks <- (hist_pois2$breaks-mean(sample))/sd(sample)

hist(stand_sample, freq = FALSE, col = "light blue", xlim = c(-4, 4), main = "Гистограмма стандартизированной выборки и плотность N(0,1)", ylab = "Плотность", breaks = stand_breaks) #xlim задаёт границы для оси x
x <- seq(-4, 4, length.out = 100) #строим последовательность из 100 чисел, заполняющих промежуток от -4 до 4 как значения аргумента для подстановки в функцию плотности
y <- dnorm(x, mean = 0, sd = 1) #подставляем x в функцию плотности и получаем вектор
lines(x, y, lwd = 2, col = "orange") #строим кривую поверх гистограммы

В данном случае можно сделать интересное наблюдение, которое соответствует действительности: при больших \(\lambda\) Пуассоновское распределение неплохо аппроксимируется нормальным распределением со средним и дисперсией, равными \(\lambda\).