Данный ноутбук является конспектом по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2019). Автор ноутбука - вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
Посмотрим на то, что уже заботали под новым углом
Немного углубим знания, выработаем интуицию
Будем смотреть презы, писать код, болтать, решать задачки на доске
Посмотрим на примеры задач, которые приходится решать аналитикам
Конечно же в R
R очень лёгок в освоении
Одна из самыз красивых визуализаций данных
Очень хорош в работе со статистикой
Много готовых пакетов и большое комьюнити
Используется многими компаниями в работе, вот устаревший мини-список
Не делаешь сам $\Rightarrow$ никогда не научишься $\Rightarrow$ $4$ большие огромные домашки, за каждую можно набрать $100$ баллов
Домашки делаются в командах по три человека
Баллы ставятся на команду
Для зачёта и оценки $10$ надо набрать на команду $100$ баллов
В каждой домашке много рыбёшек и большой кит, можно решать только то, к чему лежит душа
Две команды, набравшие наибольшее количество баллов свожу почилить с винишком, пиццей и настолками на крышу Яндекса
Команда может претендовать на крышу, если она набрала больше $100$ баллов
Мы изучали вероятности
Мы изучали события
Мы изучали случайные величины и их распределения
ЗАЧЕМ МЫ ЭТО ДЕЛАЛИ?
Лаплас: детерминизм, мы могли бы идеально прогнозировать вселенную, если бы измерили точное положение каждого атома. Издержки этого огромны.
Между совершенством природы и несовершенством человеческого познания огромный разрыв.
Неопределённость — результат этого разрыва. Случайность — это результат нашего невежества, а вероятность способ это невежество измерить.
Наука не может рассматривать вероятность как субъективную меру невежества.
Можно оценивать вероятности только тех событий, которые происходят более одного раза.
Вопрос «Какова вероятность, что кандидат N победит на выборах?» не имеет ответа, так как событие уникально и не обладает «частотой».
Сундук — различные процесс порождения данных. Теория вероятностей изучает этот сундук. В реальности мы его не видим.
В реальном мире мы видим выборки, которые сундук выплёвывает на нас. Математическая статистика изучает испражнения сундука и по ним пытается восстановить его внутренности.
Все манипуляции по восстановлению внутренностей сундука по его испражнениям позволяет делать ряд теорем.
Некоторые из них вы уже выучили!
Что это за теоремы?
ЗБЧ утверждает, что среднее арифметическое большого числа похожих случайных величин «стабилизируется» с рочтом их числа
Как бы сильно случайные величины не отклонялись от своего среднего значения, эти отклонения взаимно гасятся и среднее арифметическое приближается к постоянной величине
В обеих принимают роды, выяснилось, что в одной из больниц оценка вероятности появления мальчика составила $0.7$
В какой из больниц это произошло и почему?
Есть две больницы: маленькая и большая
Пусть $X_1, \ldots, X_n$ попарно независимые и одинаково распределённые случайные величины с конечным вторым моментом, $E(X_i^2) < \infty$, тогда имеет место сходимость:
$$ \frac{X_1 + \ldots + X_n}{n} \overset{p}{\to} E(X_1) $$Пусть $X_1, \ldots, X_n$ попарно независимые и одинаково распределённые случайные величины с конечным вторым моментом, $E(X_i^2) < \infty$, причём $\sum_{k=1}^{\infty} \frac{Var(X_k)}{k^2} < \infty$, тогда имеет место сходимость
$$ \frac{X_1 + \ldots + X_n}{n} \overset{п.н.}{\to} E(X_1) $$Примерно с $1600$-х годов ЗБЧ позволяет зарабатывать деньги на страховании
Упражнение: для $25$-летней девушки вероятность прожить ещё год составляет $0.9$. Страховка в год стоить $1000$ рублей при взносе в $110$ рублей. Какой будет средняя прибыль компании с одной страховки?
$X_i$ — прибыль с одной страховки
$X_i$ | $110$ | $-890$ |
---|---|---|
$P(\ldots)$ | $0.9$ | $0.1$ |
Средняя прибыль компании составит $\frac{1}{n} \sum X_i$. По закону больших чисел:
$$ \frac{1}{n} \sum_{i=1}^n X_i \overset{p}{\to} E(X_1) = 0.9 \cdot 110 - 0.1 \cdot 890 = 10 $$Другой ништяк, который нам разрешает ЗБЧ — генерация случайных величин для оценки разных математических ожиданий и т.п.
Не можешь посчитать? Сгенерируй!
Обычно такие генерации называют методом Монте-Карло
library("ggplot2") # Пакет для красивых графиков
# Если вы работаете в R-studio, вы можете избежать подгрузки пакетов ниже
# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
library("grid") # Пакет для субплотов
options(repr.plot.width=4, repr.plot.height=3)
Чтобы сварить в R любую случайную величину, нужно знать четыре буквы: $r$, $d$, $p$ и $q$.
rnorm
эта команда сгенерирует выборку из нормального распределенияdnorm
эта команда вычислит значение плотности в указанной точкеpnorm
эта команда находит вероятностьqnorm
эта команда находит квантилиХочу сгенерировать нормальную случайную величину $X \sim N(\mu, \sigma^2)$:
$$ f(x) = \frac{1}{\sigma \sqrt{2 \pi}} \cdot e^{-\frac{(x - \mu)^2}{2 \sigma^2}} $$x <- rnorm(4, mean=5, sd=3) # если дисперсия 9,
x # то стандартное отклонение 3
Хочу узнать $f(3)$
dnorm(3, mean=5, sd=3)
Хочу узнать $F(3) = P(X < 3) = \int_{-\infty}^3 f(x)dx$
pnorm(3, mean=5, sd=3)
Хочу узнать $P(4 < X < 9)$ - ?
Хочу узнать
$$P(4 < X < 9) = \int_4^9 f(x) dx = F(9) - F(4)$$pnorm(9, mean=5, sd=3) - pnorm(4, mean=5, sd=3)
Квантиль уровня $\gamma$ это такое число $q$, что:
$$ P(X < q) = \gamma $$qnorm(1 - 0.05/2, mean=0, sd=1)
x <- c(0.95, 0.975, 0.995)
qnorm(x, mean = 0, sd = 1)
# квантили для другого распределения!
x <- c(0.95, 0.975, 0.995)
qt(x, df = 10)
Что сгенерирует код?
rnorm(10)
x <- rnorm(1000, mean=5, sd=3)
mean(x) # среднее
var(x) # выборочная дисперсия
sd(x) # выборочное стандартное отклонение
median(x) # выборочная медиана
library("ggplot2") # пакет для красивых картинок
x <- rnorm(1000, mean=5, sd=3)
qplot(x) # гистограма
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
library("ggplot2") # пакет для красивых картинок
x <- seq(-5, 15, by=0.01)
y <- dnorm(x, mean=5, sd=3)
qplot(x, y, geom="line")
library("ggplot2") # пакет для красивых картинок
x <- seq(-5, 15, by=0.01)
y <- pnorm(x, mean=5, sd=3)
qplot(x, y, geom="line")
Чтобы решать реальные проблемы! Например, пусть $X \sim N(5,3)$. Как найти $E \left(\frac{1}{X} \right)$?
Спосбо второй (ЗБЧ разрешает нам так делать):
n_obs <- 10^6 # число наблюдений
x <- rnorm(n_obs, mean = 5, sd = 3)
mean(1/x)
$X_1, X_2, X_3 \sim U[0;2]$, независимые
Хотим знать $P(X_1 + X_2 + X_3^2 > 5)$
n_obs <- 10^6
x_1 <- runif(n_obs, min = 0, max = 2)
x_2 <- runif(n_obs, min = 0, max = 2)
x_3 <- runif(n_obs, min = 0, max = 2)
success <- x_1 + x_2 + x_3^2 > 5
success[1:5]
sum(success) / n_obs
$X_1, X_2, X_3 \sim U[0;2]$, независимые
Хотим знать $$P(X_1 + X_2 > 0.8 \mid X_3 < 0.1)$$
n_obs <- 10^6
x_1 <- runif(n_obs, min = 0, max = 2)
x_2 <- runif(n_obs, min = 0, max = 2)
x_3 <- runif(n_obs, min = 0, max = 2)
uslovie <- x_3 < 0.1
x_1[1:3]
uslovie[1:3]
x_1[uslovie][1:3]
n_obs <- 10^6
x_1 <- runif(n_obs, min = 0, max = 2)
x_2 <- runif(n_obs, min = 0, max = 2)
x_3 <- runif(n_obs, min = 0, max = 2)
uslovie <- x_3 < 0.1
success <- x_1[uslovie] + x_2[uslovie] > 0.8
sum(success)/n_obs
sample(1:10, size = 8) # выборка без повторений
sample(1:10, size = 8, replace = TRUE) # с повторениями
# неправильная монетка
sample(c("Орёл", "Решка"), size = 5, replace = TRUE, prob = c(0.3, 0.7))
x <- sample(c("Орёл", "Решка"), size = 10^6, replace = TRUE, prob = c(0.3, 0.7))
sum(x == 'Орёл')/length(x)
При определённых условиях сумма достаточно большого числа случайных величин имеет распределение близкое к нормальному
Главное: чтобы случайные величины были похожи и не было такого, что одна резко выделяется на фоне остальных
Есть много разных ЦПТ с разными условиями
Пусть $X_1, \ldots, X_n, \ldots$ — последовательность независимых одинаковых случайных велчин с конечным вторым моментом $E(X_i^2) < \infty$. Тогда
n = 10^4
X1 <- runif(n, min = -1, max = 1)
X2 <- runif(n, min = -1, max = 1)
X3 <- runif(n, min = -1, max = 1)
X4 <- runif(n, min = -1, max = 1)
qplot(X1, bins = 70)
qplot(X1 + X2, bins = 70)
qplot(X1 + X2 + X3, bins = 70)
qplot(X1 + X2 + X3 + X4, bins = 70)
В ЗБЧ:
$$ \frac{X_1 + \ldots + X_n - n E(X_1)}{n} \overset{p}{\to} 0 $$В ЦПТ та же дробь домножается на $\sqrt{n}$, это замедляет сходимость и мы приходим к более интересному результату:
$$ \sqrt{n} \cdot \frac{X_1 + \ldots + X_n - n E(X_1)}{n} \overset{d}{\to} N(0, Var(X_1)) $$ЦПТ и ЗБЧ работают в среднестане
А что, если какая-то одна случайная величина выбивается?
Тогда мы перемещаемся из среднестана в крайнестан и сталкиваемся с проблемой тяжёлых хвостов
О тяжёлых хвостах мы ещё поговорим, они часто выскакивают в мире финансов
Две монеты, вероятности орла равны $0.6$ и $0.4$. Камила подбрасывает первую монетку до появления орла. Вика вторую.
а) Найдите с помощью симуляций вероятность того, что Вика сделает больше подбрасываний, чем Камила
а) Найдите с помощью симуляций вероятность того, что Вика сделает больше подбрасываний, чем Камила
n_obs = 10^6
p1 = 0.6
p2 = 0.4
# симулируем подбрасывания Вики и Камилы
x_kamila = rgeom(n_obs, prob = p1)
x_vika = rgeom(n_obs, prob = p2)
x_kamila[1:10] # в векторе стоит число подбрасываний до первого орла
sum(x_vika > x_kamila)/n_obs
б) Камила сделает больше подбрасываний, чем Вика;
n_obs = 10^6
p1 = 0.6
p2 = 0.4
# симулируем подбрасывания Вики и Камилы
x_kamila = rgeom(n_obs, prob = p1)
x_vika = rgeom(n_obs, prob = p2)
x_kamila[1:10] # в векторе стоит число подбрасываний до первого орла
sum(x_vika < x_kamila)/n_obs
в) Камила и Вика сделают одинаковое число подбрасываний
sum(x_vika == x_kamila)/n_obs
г) Верно ли, что в сумме эти три вероятности дают единицу? Логично ли это?
sum(x_vika > x_kamila)/n_obs + sum(x_vika < x_kamila)/n_obs + sum(x_vika == x_kamila)/n_obs
Иван Фёдорович Крузенштерн (ШТО?!) случайным образом с возможностью повторов выбирает $10$ натуральных чисел от $1$ до $100$. Пусть $X$ — минимум из этих чисел, а $Y$ — максимум. С помощью симуляций оцените:
а) $P(Y > 3X)$
natural <- sample(1:100, size = 5, replace = TRUE)
natural
min(natural) # случайная величина X
max(natural) # случайная величина Y
а) $P(Y > 3X)$
n_obs = 10^6
x = rep(0, n_obs)
y = rep(0, n_obs)
for(i in 1:n_obs){
natural = sample(1:100, size = 10, replace = TRUE) # сгенерировали выборку
x[i] = min(natural) # нашли минимум и максимум для неё
y[i] = max(natural)
}
x[1:5]
success = y > 3*x
success[1:5]
sum(success) / n_obs
б) $E(X \cdot Y)$
mean(x*y)
в) $P(Y > 3X \mid Y < X^2)$
z <- c(1,2,3,4,5,6,7)
z > 4
z[z > 4]
usl = y < x^2
usl[1:5]
success <- y[usl] > 3*x[usl]
sum(success) / sum(usl)
г) $E(X \cdot Y \mid Y < X^2)$
mean(x[usl]*y[usl])
д) $E \left( \frac{X}{X + 2Y} \right)$
z = x/(x + 2*y)
mean(z)
е) $Corr(X,Y)$
cor(x,y)
Известно, что в пятизначном номере телефона все цифры разные. Какова вероятность того, что при этом условии среди них только одна цифра чётная (номер может начинаться с нуля)?
# генерируем номер телефона
sample(size = 5, 0:9, replace = FALSE)
# проверяем все цифры на чётность
sample(size = 5, 0:9, replace = FALSE) %% 2
# смотрим сколько в сумме нечётных цифр
sum(sample(size = 5, 0:9, replace = FALSE) %% 2)
# и всё это в цикле!
n_obs = 10^6
m = 0
for( i in 1:n_obs){
# если 4 нечётные, значит одна чётная
if(sum(sample(size = 5, 0:9, replace = FALSE) %% 2) == 4){
m = m + 1
}
}
m/n_obs
Для $20$ участников соревнований, среди которых $8$ российских, в гостинице забронировано $20$ номеров. Из них $12$ с видом на море. Участники соревнований наугад получают ключи от номернов. Найдите вероятность того, что номера с видом на море достанутся всем российским спортсменам.
# вектор номеров
room = c(rep('море',12),rep('пустошь',8))
room[1:5]
sample(room)[1:5] # перемешаем номера
# Будем считать что первые 8 позиций вектора относятся к россиянам
sample(room)[1:8] == 'море'
n_obs = 10^6
# комнаты
room = c(rep('море',12),rep('пустошь',8))
m = 0
for(i in 1:n_obs){
if(sum(sample(room)[1:8] == 'море') == 8){
m = m + 1
}
}
m/n_obs
Игральная кость. Нам хочется узнать какое среднее значение на ней будет выпадать при увеличении количества бросков.
sample(1:6, size = 1, replace = TRUE)
sample(1:6, size = 2, replace = TRUE)
sample(1:6, size = 3, replace = TRUE)
sample(1:6, size = 4, replace = TRUE)
n_obs = 500
kubik = rep(0, n_obs)
for(n in 1:n_obs){ # n - количество подбрасываний кубика
# Сделали выборку
smpl = sample(1:6, size = n, replace = TRUE)
# Подсчитали среднее значение и запомнили его
kubik[n] = mean(smpl)
}
kubik[1:5]
qplot(1:n_obs, kubik, geom='line', xlab='number of observations', ylab='result')
df = data.frame(n = 1:n_obs, kubik = kubik, mean = rep(3.5,n_obs))
head(df,2)
n | kubik | mean |
---|---|---|
1 | 4 | 3.5 |
2 | 6 | 3.5 |
# https://colorscheme.ru/color-converter.html
ggplot(df, aes(n)) +
geom_line(aes(y=kubik), color="#0aa120", alpha = 0.6) +
geom_line(aes(y=mean), color="blue",size=1.2) +
xlab('number of observations') + ylab('result') +
ggtitle('Law of large numbers')
Смешная предыстория, которая не влезла на слайд.
Покупать газету Миша решил по $15$ рублей, а продавать за $30$. Количество потенциальных покупателей - случайная величина с распределением Пуассона. Опытным путём было установлено, что среднее значение этой величины равно $50$. Нераспроданные газеты ничего не стоят. Пусть $n$- количество газет, максимизирующее ожидаемую прибыль Мишы.
С помощью компьютера найдите оптимальное значение $n$ и ожидаемую прибыль.
n_obs = 10^6 # число наблюдений
cost = 15 # издержки на покупку газеты
price = 30 # цена, по которой будем продавать газету
x = rpois(n_obs, lambda = 50) # среднее значение совпадает с lambda
x[1:5]
profit = x*price - 20*cost
profit[1:5]
# найдём среднюю прибыль
mean(profit)
n_obs = 10^6 # число наблюдений
profits = rep(0,100) # будем записывать средние прибыли
for(s in 1:100){
x = rpois(n_obs, lambda = 50) # генерируем покупочки
x[x > s] = s # нельзя продать лишнее
profits[s] = mean(x*price - s*cost) # запоминем среднюю прибыль
}
qplot(1:100, profits, geom='line')
max(profits)
which.max(profits)
sqrt(var(profits))
n_obs = 10^6 # число наблюдений
profits = rep(0,100) # будем записывать средние прибыли
std_error = rep(0,100) # сюда будем записыват ошибку
for(s in 1:100){
x = rpois(n_obs, lambda = 50) # генерируем покупочки
x[x > s] = s # нельзя продать лишнее
profits[s] = mean(x*price - s*cost) # запоминем среднюю прибыль
std_error[s] = sqrt(var(x*price - s*cost)) # запоминаем дисперсию
}
df = data.frame(s = 1:100, profit = profits, std_error = std_error)
head(df,2)
s | profit | std_error |
---|---|---|
1 | 15 | 0 |
2 | 30 | 0 |
ggplot(df, aes(s)) +
geom_line(aes(y=profit), colour="blue") +
geom_ribbon(aes(ymin = profit - 2*std_error,
ymax = profit + 2*std_error), alpha=0.2)
set.seed(42)
rnorm(1, mean = 2, sd = 3)
rnorm(1, mean = 2, sd = 3)
picture <- function(x){
df = data.frame(sample = x)
# Размеры картинки
options(repr.plot.width=8, repr.plot.height=3)
binwidth = 1 # ширина бинов
p1 = ggplot(df, aes(x = sample, binwidth = binwidth))+
# Наносим гистограмму
geom_histogram(aes(y=..density..), binwidth = binwidth, colour = "white",
fill = "cornflowerblue", size = 0.1)
# цвет линий разделителей, заливка, толщина линий разделителей
p2 = ggplot(df, aes(x = sample, binwidth = binwidth))+
stat_ecdf(color = "darkred", size = 1)
# Располагаем графики рядом. Этот код нужен только для юпитерской тетрадки.
pushViewport(viewport(layout = grid.layout(1, 2)))
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
}
n = 1000
size = 20
p = 0.5
x <- rbinom(n, size = size, prob = p)
picture(x)
n = 1000
# Распределение Пуассона, Pois(lambda)
lambda =3
x <- rpois(n, lambda = lambda)
picture(x)
Распределение Пуассона часто используется для моделирования "событий-счётчиков", при этом $\lambda$ это интенсивность потока событий.
Если у нас ест несколько суммирующихся между собой потоков, их интенсивности суммируются. То есть, если
n = 1000
# Нормальное распределение, N(mu, sigma)
mu = 0; sigma = 1
x <- rnorm(n, mean = mu, sd = sigma)
picture(x)
n = 10000
# Экспоненциальное распределение, Exp(alpha)
alpha = 0.1
x <- rexp(n, rate = alpha)
picture(x)
# Равномерное распределение, U[mn;mx]
n = 10000
mn=0; mx=24
x <- runif(n, min = mn, max = mx)
picture(x)