Данный ноутбук является домашкой по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2017-2018). Автор ноутбука вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
Приветствую вас внутри первой домашки. У каждой домашки должны быть свои герои (но это неточно). Моими героями являетесь вы! Вы - герои не потому, что я хочу кого-то простебать или обидеть, а потому что я всех вас очень люблю и хочу увековечить. Кто-то называет именами любимых людей вновь открытые звёзды и острова, кто-то вновь созданный вордовский файл для нира, а кто-то делает их героями своих историй. Я не Джордж Мартин, в моих историях никто не погибнет.
Краткий брифинг:
Ближе к делу. С чего начинается любой скрипт? Правильно! С подгрузки пакетов :)
library("ggplot2") # Пакет для красивых графиков
library("grid") # Пакет для субплотов
# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
options(repr.plot.width=4, repr.plot.height=3)
# Внимание! Если вы делаете дз в Rstudio, то вам не нужны пакеты grid, repr и т.п.
# Вам нужен только пакет ggplot2!
rep(2,5) # повторяет число 2 целых 5 раз
seq(0,1,0.2) # выдаёт все числа от 0 до 1 с шагом 0.2
x <- c(1,2,3)
y <- c(2,-1,4)
z <- c(3,5,2)
pmin(x,y,z) # находит поэлементные минимумы во всех векторах
which.max(x) # выдаёт позицию максимума
В R куча встроенных наборов данных. Например, в наборе данных chickwts
лежат веса куриц и тип корма, который используется для их выращивания.
head(chickwts,5) # команда head позволяет посмотреть на
# первые пять строк таблички
weight | feed |
---|---|
179 | horsebean |
160 | horsebean |
136 | horsebean |
227 | horsebean |
217 | horsebean |
Постройте гистограмму с распределением веса куриц с 20 столбиками.
# легче лёгкого!
qplot(chickwts$weight, bins=20)
Постройте эмпирическую функцию распределения для веса куриц.
# скопипастим код из лекции
ggplot(chickwts, aes(x = weight))+ # на какой таблице строи график
stat_ecdf( ) # какой именно график строим
# Можно сделать график чуть красивше! Розовым, пунктирным и большим!
ggplot(chickwts, aes(x = weight))+
stat_ecdf(linetype = 'dashed', color = "magenta", size = 5)
Найдите базовые характеристики распределения:
# сохраним вес курицы в отдельный вектор, чтобы было проще
x <- chickwts$weight
mean(x) # среднее
median(x) # медиана
var(x) # дисперсия
sd(x) # стандартное отклонение (корень из дисперсии)
q = c(0.05, 0.5, 0.95)
quantile(x,q)
# конечно правда! Медиана это и есть 50% квантиль
median(x) == quantile(x, 0.5, names = FALSE) # names отвечает за подписывание имён
Сделали упражнение? Ответили на вопросы? Вы офигенны! Вы уже можете больше, чем среднестатистический аудитор из большой четвёрки ;) Двигайтесь дальше в таком же темпе.
Случайные величины $X_1, \ldots, X_5$ имеют равномерное распределение на отрезке $[0;1]$ и независимы. С помощью симуляций оцените:
а) $P(\min\{X_1, \ldots, X_5\} > 0.2)$
б) $P(\min\{X_1, \ldots, X_5\} > 0.2 \mid X_1 + X_2 < 0.5)$
в) $E(\min\{X_1, \ldots, X_5\})$
г) $E(\min\{X_1, \ldots, X_5\} \mid X_1 + X_2 < 0.5)$
Насколько логично то, что вы получили? Попробуйте порассуждать на эту тему. Постройте гистограмму для распределения минимума из этих случайных величин гистограмму для распределения максимума из этих случайных величин. Вспомните как выглядит теоретическая плотность распределения максимума и теоретическая плотность распределения минимума. Мы ещё с ними столкнёмся ;)
# немного осмысленной копипасты из лекции
n_obs <- 10^6 # число наблюдений
# пяток случайных величин
x_1 <- runif(n_obs, min = 0, max = 1)
x_2 <- runif(n_obs, min = 0, max = 1)
x_3 <- runif(n_obs, min = 0, max = 1)
x_4 <- runif(n_obs, min = 0, max = 1)
x_5 <- runif(n_obs, min = 0, max = 1)
# находим минимальное значение в каждой пятёрке чисел
mn <- pmin(x_1, x_2, x_3, x_4, x_5)
qplot(mn) # посмотрим на гистограмму, интересно же!
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
success <- mn > 0.2 # TRUE будет стоять в тех местах, где это правда
sum(success) / n_obs # а вот и итоговая оценка вероятности! :)
# также легко решаем пункт b)
uslovie = x_1 + x_2 < 0.5 # находим позиции, где выполнено условие
success <- mn[uslovie] > 0.2 # делаем срез по этому условию
sum(success) / n_obs # наконец получаем оценку вероятности
# Вполне логично, что она уменьшилась по сравнению с предыдущей, мы же наложили ограничения!
# решаем пункт c)
mean(mn)
# и также пункт d)
# Вполне логично, что математическое ожидание меньше, у нас же ограничения!
mean(mn[uslovie])
Иван Фёдорович Крузенштерн (ШТО?!) случайным образом с возможностью повторов выбирает $10$ натуральных чисел от $1$ до $100$. Пусть $X$ — минимум из этих чисел, а $Y$ — максимум. С помощью симуляций оцените
а) $ P(Y > 3X)$
б) $E(X \cdot Y)$
в) $P(Y > 3X \mid Y < X^2)$
г) $E(X \cdot Y \mid Y < X^2)$
д) $E \left( \frac{X}{X + 2Y} \right)$
е) $Corr(X,Y)$
# Ивану Фёдоровичу грех не помочь!
# Основная сложность этой задачи в генерации выборки, давайте попрубуем
# разобраться пошагово и для начала сгенерируем только одну выборку
natural <- sample(1:100, size = 5, replace = TRUE)
natural
min(natural) # случайная величина X
max(natural) # случайная величина Y
# нам таких выборок понадобится очень много, придётся написать цикл для их генерации...
n_obs = 10^6
x = rep(0, n_obs) # создаём два вектора из нулей, так будет быстрее работать
y = rep(0, n_obs) # чем ежели создавать пустой вектор (почему объясню на паре)
# в вектор x будем записывать минимумы, в y максимумы
for(i in 1:n_obs){
natural = sample(1:100, size = 10, replace = TRUE) # сгенерировали выборку
x[i] = min(natural) # нашли минимум и максимум для неё
y[i] = max(natural)
}
# Вот и вся генерация, осталось только ответить на вопросы задачи
# a)
success = y > 3*x
sum(success) / n_obs
# б)
mean(x*y)
# в)
usl = y < x^2
success <- y[usl] > 3*x[usl]
sum(success) / n_obs
# г)
mean(x[usl]*y[usl])
# д)
z = x/(x + 2*y)
mean(z)
# e)
cor(x,y)
Миша поступил на первый курс экономического факультета, потому что у него в крови всегда была предпренимательская жилка. На первом курсе она довела его, и он решил учёбу забросить и стать комерсом. На паблике "Бизнес-молодость" Миша подсмотрел шикарную бизнес-модель и решил, что будет торговать еженедельными газетами. На первом курсе Миша прослушал курс микроэкономики и понимает, что в условиях совершенной конкуренции, цену определяет рынок. Миша никак не может влиять на неё, однако он может влиять на то, какое количество газет следует закупить. На втором курсе Миша не учился. Из-за этого он не знает теорию вероятностей и вынужден для решения своих повседневных проблем нанять аналитика.
Покупать газету Миша решил по $15$ рублей, а продавать за $30$. Количество потенциальных покупателей - случайная величина с распределением Пуассона. Опытным путём было установлено, что среднее значение этой величины равно $50$. Нераспроданные газеты ничего не стоят. Пусть $n$- количество газет, максимизирующее ожидаемую прибыль Мишы.
a) С помощью компьютера найдите оптимальное значение $n$ и ожидаемую прибыль.
b) Правильно ли поступил Миша, что бросил учёбу и ушёл в коммерцию?
# Ох уж этот Миша!
# Найдём оптимальное количество газет для закупки!
n_obs = 10^6 # число наблюдений
cost = 15 # издержки на покупку газеты
price = 30 # цена, по которой будем продавать газету
x = rpois(n_obs, lambda = 50) # среднее значение совпадает с lambda
# посмотрим на прибыль, еслм мы будем продавать 30 газет!
# покупатели не могут купить больше 30 газет
x[x > 30] = 30 # заменим всё, что больше 30 на 30
# найдём все прибыли
profit = x*price - 30*cost
# найдём среднюю прибыль
mean(profit)
# Отлично! Дело осталось за малым, перебрать все возможноые варианты покупок и выбрать оптимальный
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)
# b) На этот пункт у меня для вас ответа нет.
Глеб раздобыл кубик и собирается как следует изучить его свойства. Для этого он подбросил его $n$ раз. Величина $X_1$ — число выпадений единицы, а $X_6$ — число выпадений $6$. Глебу интересно какова корреляцие между этими двумя случайными величинами. Что вам на этот счёт подсказывает интуиция? Какой знак будет у этой корреляии: положительным или отрицательным?
Оцените с помощью симуляций $Corr(X_1,X_6)$. Пусть $n=10$. Обратите внимание, что $n$ в данном случае это не количество симуляций. Десять подбрасываний это всего лишь одно испытание. После симуляций попытайтесь найти эту корреляцию руками. на бумажке.
# подбросим кубик 10 раз и вычислим значения X1 и X6
smpl = sample(size = 10, 1:6, prob = rep(1/6,6), replace = TRUE)
sum(smpl == 1) # это X1
sum(smpl == 6) # это X6
n_obs = 10^6 # сделаем 10 подбрасываний n раз
x1 = rep(0, n_obs) # сюда будем записывать количество 1 в каждой серии испытаний
x6 = rep(0, n_obs) # сюда будем записывать число 6 в каждой серии испытаний
for(i in 1:n_obs){
smpl = sample(size = 10, 1:6, prob = rep(1/6,6), replace = TRUE) # первая серия испытаний!
x1[i] = sum(smpl == 1)
x6[i] = sum(smpl == 6)
}
cor(x1, x6)
Случайная величина $X$ распределена равномерно на отрезке $[0;1]$. Женя подкидывает эту случайную величину. Она принимает значение $x$. После подбрасывания Женя изготавливает монетку, которая выпадает «орлом» с вероятностью $x$ и передаёт её Дане. Даня, не зная $x$, подкидывает монетку один раз. Она выпала «орлом».
a) Какова вероятность, что она снова выпадет «орлом»?
# смодулируем одну итерацию всего этого процесса
x = runif(1, min = 0, max = 1) # Женя делает монетку
cat('Вероятность орла на монетке:',x) # Команда cat это более пафосный print
# Даня два раза подкидывает монетку
y = sample(c('О','Р'), prob = c(x, 1-x), size = 2, replace = TRUE)
cat('\nУ Дани выпало:', y)
Вероятность орла на монетке: 0.7324506 У Дани выпало: О О
# Отлично! Теперь мы должны сделать как можно больше разных монет и подкинуть их как можно больше раз.
# Из-за того, что наш эксперимент происходит последовательно, придётся подкидывать каждую монетку Жени много раз.
n_zenia = 100 # сколько монеток сделал Женя
x = runif(n_obs, min = 0, max = 1) # Женя делает свой ход n_zenia раз
x[1:5] # первые пять монеток, созданные Женей :)
# Знаете, иногда так бывает, что тебе лень что-то делать и ты откладываешь это до лучших времён.
# Вот и я отложил. Сейчас, когда я пишу это, до пары остаётся два часа. Можно было бы взять и
# написать тут два цикла для генерации вероятностей вместо это текста, но я не хочу их писать.
# Они будут большими и страшными. И вообще я бы эту задачку лучше бы убрал из домашки, так как
# она для генераций не очень удачная.
b) Как выглядит ответ, если Дане известно, что монетка при $n$ подбрасываниях $k$ раз выпадала орлом?
# Ня! Решение пропущено. Или отложено до лучших времён. Не знаю.
А сможете найти эти вероятности вручную?
Пусть $X \sim Exp(2)$. Аня помнит, что такое распределение называется экспоненциальным и что в такой ситуации функция распределения случайной величины $X$ выглядит как:
$$ F(x) = 1 - e^{-2 x}$$Аня хочет понять какое распределение будет у случайной величины $Y$, если
$$ Y = 1 - e^{-2 X}$$Попробуйте ссимулировать такую случайную величину, посмотреть на её гистограмму и дать ответ на вопрос Ани. Подумайте почему ответ получился именно таким.
n_obs = 10^6
x = rexp(n_obs, rate = 2)
y = 1 - exp(-2*x)
qplot(y) # равномерная на отрезке [0; 1]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Маша и Виолетта живут в общежитии в одной комнате. Маша постоянно что-то говорит про нейронные сети. Виолетта не очень понимает о чём идёт речь, но хотела бы разобраться, чтобы быть с подругой на одной волне. Для этого она залезла в мудрую книгу и стала читать. В одной из первых глав она узнала, что в качестве функции активации в нейронах часто используют сигмоида, которая представляет из себя функцию распределения для логистической случайной величины. Функция выглядит как-то так:
$$ F(x) = \frac{1}{1 + e^{-x}} $$Виолетта пока не поняла что такое нейрон, сигмоида и функция активации, зато она очень хорошо знает что такое распределение случайной величины. Ей жутко интересно как оно выглядит.
а) Помогите Виолетте сгенирировать случайную величину $X$, имеющую логистическое распределение. Постройте для неё гистограмму, плотность и функцию распределения прямо как на паре.
б) Также Виолетте хотелось бы узнать какое распределение будет у случайной величины $Y = F(X)$. Попробуйте посмотреть на него с помощью симуляций. Как следует подумайте почему ответ получился именно таким. Функция rlogis
с параметрами location = 0
и scale = 1
вам в помощь!
# a)
n_obs = 10^3 # не будем генерировать особо много, а то функция не построится и всё зависнет
x = rlogis(n_obs, location=0, scale=1) #сгенерировали
# построим функцию распределения
df = data.frame(sample = x)
ggplot(df, aes(x = sample)) + stat_ecdf( )
# построи гистограмму и плотность
ggplot(df, aes(x = sample))+
# Наносим гистограмму
geom_histogram(aes(y = ..density..))+ # опция ..density.. нормирует площадь гистограмы к 1
# попробуйте стереть её и посмотреть что будет
# Наносим плотность распределения
stat_function(fun = dlogis, args = c(0,1), color = "darkred", size=1)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# b) сгенерируем выборку в нормальном объёме
n_obs = 10^6
x = rlogis(n_obs, location=0, scale=1) #сгенерировали
y = 1/(1 + exp(-x))
qplot(y) # равномерная на отрезке [0; 1]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Воу, воу, воу, теория вероятностей, палехче! В обеих задачах получился одинаковый ответ. Видимо, в случайной величине $Y = F(X)$, где $F(x)$ - функция распределения случайной величины $X$, есть какая-то магия. Подумайте о том какая и постарайтесь доказать наличие этой магии на бумаге. На семинаре мы попробуем сделать это совместно.
Было бы нечестно, если после задач $7.1$ и $7.2$ следовала бы задачка с номером $8$. Поэтому тут задачка с номером $9$!
Есть две монеты, у которых вероятности выпадения орла при однократном подбрасывании равны $0.6$ и $0.4$. Камила подбрасывает первую монетку до появления орла. Вика делает то же самое со второй монеткой. Найдите с помощью симуляций вероятность того, что
a) Вика сделает больше подбрасываний, чем Камила;
b) Камила сделает больше подбрасываний, чем Вика;
c) Камила и Вика сделают одинаковое число подбрасываний.
d) Верно ли, что в сумме эти три вероятности дают единицу? Логично ли это?
e) Найдите математическое ожидание суммарного количества бросков.
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] # в векторе стоит число подбрасываний до первого орла
# a) # логично, Вике больше повезло с монткой!
sum(x_vika > x_kamila)/n_obs
# b)
sum(x_vika < x_kamila)/n_obs
# c)
sum(x_vika == x_kamila)/n_obs
# d)
# да, эти вероятности дадут в сумме единицу, потому что других исходов не бывает
sum(x_vika > x_kamila)/n_obs + sum(x_vika < x_kamila)/n_obs + sum(x_vika == x_kamila)/n_obs
# e)
mean(x_vika + x_kamila)
Продемонстрируйте с помощью симуляций, что у геометрического распределения нет памяти. Определение отсутствия памяти ищи тут. Что это означает для Вики и Камилы?
Читаем определение! У распределения нет памяти, если количество прошлых неудач не влияет на количество будущих неудач
$$ P( Y > m + n \mid Y \ge m) = P(Y > n).$$Ровно это и закодим.
n_obs <- 10^7
p <- 0.2
y <- rgeom(n_obs, prob = p)
sum(y > 3)/n_obs # Правая часть равенства
# Левая часть равенства
y = y[y >= 2]
sum(y > 2 + 3)/length(y)
Эту задачку тоже можно решить руками на бумаге. Слабо?
Света очень долго сидела и думала. В конечном счёте у неё возникло очень много разнообразных вопросов. Помогите Свете ответить на них.
a) Пусть $X \sim N(5,9)$, а $Y \sim N(3, 16)$. Верно ли, что сумма нормальных случайных величин будет нормальной случайной величиной? Какими будут математическое ожидние и дисперсия у итоговой случайной величины. Проверьте это с помощью симуляций.
___ВНИМАНИЕ:___ код с симуляциями ни в коем случае не является доказательством верных фактов, их строгие доказательства можно найти в лекциях. Но при этом он является опровержением неверных фактов. Помните о том, что в математике да это всегда да, а нет это хотя бы один раз нет.
n_obs = 10^4
x <- rnorm(n_obs, mean=5, sd=3) # обратите внимание, что мы указываем sd а не дисперсию!
y <- rnorm(n_obs, mean=3, sd=4)
mean(x+y) # 8 = 5 + 3 всё как в лекциях!
var(x+y) # 25 = 9 + 16 снова всё как в лекциях!
qplot(x + y) # посмотрим на гистограмку :3
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
b) Пусть $X \sim N(\mu, \sigma^2)$. Верно ли, что $\frac{X-\mu}{\sigma} \sim N(0,1)$? Проверьте это с помощью симуляций.
# конечно верно! такая процедура называется стандартизацией случайной величины
# хорошо бы попробовать код для разных примеров нормальных случайных величин
x = rnorm(n_obs, mean=5, sd=3)
y = (x - 5)/3
mean(y)
var(y)
qplot(y)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
с) Пусть $X \sim U[0;1]$ и $Y \sim U[0;1]$. Верно ли, что $X + Y$ также будет иметь равномерное распределение? Проверьте это с помощью симуляций?
# конечно неверно, это чушь! сумма равномерных величин иметт распределение треугольника
x <- runif(n_obs, min = 0, max = 1)
y <- runif(n_obs, min = 0, max = 1)
qplot(x+y)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
d) Пусть $X \sim U[0;1]$. Какое распределение будет у $Y = 5 + 3\cdot X$? Узнайте это с помощью симуляций.
# она будет равномерной на отрезке от 5 до 8 :)
# мы сдвинули распределение вправо на 5 и растянули в три раза
x <- runif(n_obs, min = 0, max = 1)
y = 5 + 3*x
qplot(y)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Когда сумма случайных величин из одного семейства даёт случайную величину из этого же семейства, то такое распределение называется устойчивым относительно суммирования. Получилось ли, что нормальное и равномерное распределения устойчивы относительно суммирования? Как думает распределение Пуассона, экспоненциальное и Хи-квадрат устойчивы относительно суммирования?
Найдите в лекциях по теории вероятностей весь материал, который касается свойств нормального распределния и повторите его.