Данный ноутбук является домашкой по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2019). Автор ноутбука - вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
Приветствую вас внутри третьей домашки. Краткий брифинг:
Как оформлять домашки:
На странице курса висит видео-инструкция по оформлению.
Ближе к делу. С чего начинается любой скрипт? Правильно! С подгрузки пакетов :)
library('maxLik') # пакет для метода макс. правдоподобия
library('dplyr') # пакет для работы с таблицами
library("ggplot2") # Пакет для красивых графиков
library("grid") # Пакет для субплотов
# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
options(repr.plot.width=4, repr.plot.height=3)
Практически каждая задачка строится по принципу от простого к сложному. В конце почти каждой задачи есть небольшая мораль с заделом на будущее. Если вы в какой-то момент осознаёте, что пошёл хард, и непонятно что делать. Либо бросайте эту задачу и переходите к началу следующей, либо, если вы настоящий самурай, пишите мне в личку. Никого к чёрту не пошлю, всем отвечу и помогу.
Начнём наше погружение в моделирование с небольшого разогрева. Пишите в этих задачках сразу же хороший код. Вы сможете затем его переиспользовать в более сложных задачах. Например, код для нормального распределения понадобится вам около $3-4$ раз.
Оцените методом максимального правдоподобия какие сериалы Филипп посмотрел за последний год. Для каждого сериала дайте обоснование, почему вы решили, что я его смотрел. Каждый угаданный сериал принесёт вам дополнительный балл. Каждый неправильно предположенный сериал отнимет у вас один балл. Правильно угаданные сериалы без обоснования тоже отберут у вас один балл.
# кажется, что тут не нужно писать код
Пусть $X \sim Exp(0.05)$. Сгенерируйте из этого распределения выборку размера $100$, а потом методом максимального правдоподобия оцените параметр $\alpha$.
# ваш код
Постройте график для функции правдоподобия. По оси $x$ отложите значения $\alpha$, по оси $y$ логарифм правдоподобия. Отметьте на картинке точку оптимума.
# легко, но надо посидеть и разобраться в том как строятся картинки с пары
Пусть $X \sim N(5, 4)$. Сгенерируйте из этого выборку размера $100$, затем методом максимального правдоподобия оцените параметры $\mu$ и $\sigma^2$. Постройте для параметра $\mu$ асимптотический $80\%$ доверительный интервал. Нарисуйте такую же картинку, как в предыдущей задачке. Нанесите на неё доверительный интервал.
# легко, надо только над картинкой посидеть
Проверьте тестом отношения правдоподобий гипотезу о том, что $\mu = 4$ на уровне значимости $5\%$.
# изи, переиспользуем код из 2 задачи
В условиях той же задачки тестом отношения правдоподобий проверьте гипотезу, что $\mu = 4$ и $\sigma = 3$. Обратите внимание на число ограничений и количество степеней свободы.
# изи
Давайте немного повысим ставки. Будет всё ещё просто, но пунктов больше.
В табличке лежит классический датасет с пассажирами, которые плыли на Титанике. Кто-то из них умер, кто-то остался в живых. В будущем, вы будете решать на этом датасете свою первую в жизни задачку классификации. А пока что давайте попробуем поанализировать, как часто на титанике умирают люди.
df = read.csv('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/titanic.csv', sep=',')
head(df)
X | pclass | survived | name | sex | age | sibsp | parch | ticket | fare | cabin | embarked | boat | body | home.dest |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 1st | 1 | Allen, Miss. Elisabeth Walton | female | 29.0000 | 0 | 0 | 24160 | 211.3375 | B5 | Southampton | 2 | NA | St Louis, MO |
2 | 1st | 1 | Allison, Master. Hudson Trevor | male | 0.9167 | 1 | 2 | 113781 | 151.5500 | C22 C26 | Southampton | 11 | NA | Montreal, PQ / Chesterville, ON |
3 | 1st | 0 | Allison, Miss. Helen Loraine | female | 2.0000 | 1 | 2 | 113781 | 151.5500 | C22 C26 | Southampton | NA | Montreal, PQ / Chesterville, ON | |
4 | 1st | 0 | Allison, Mr. Hudson Joshua Crei | male | 30.0000 | 1 | 2 | 113781 | 151.5500 | C22 C26 | Southampton | 135 | Montreal, PQ / Chesterville, ON | |
5 | 1st | 0 | Allison, Mrs. Hudson J C (Bessi | female | 25.0000 | 1 | 2 | 113781 | 151.5500 | C22 C26 | Southampton | NA | Montreal, PQ / Chesterville, ON | |
6 | 1st | 1 | Anderson, Mr. Harry | male | 48.0000 | 0 | 0 | 19952 | 26.5500 | E12 | Southampton | 3 | NA | New York, NY |
В колонке survived
лежат нолики и единицы. Единичка означает, что человек выжил.
[а] Случайная величина survived
для каждого человека имеет распределение Бернулли. Выпишите логорифмическое правдоподобие. Вбейте функцию в R.
# изи
[б] Найдите оценку параметра $p$ методом максимального правдоподобия на компьютере.
# изи, пакет maxLik и семинар в помощь
[в] Постройте для параметра $p$ асимптотический $80\%$ доверительный интервал.
# по матрице гессе, как на паре
[г] Нарисуйте картинку, где по оси $x$ будет отложен парметр $p$, а по оси $y$ значения функции правдоподобия. Отметьте на этой картинке точку оптимума, а также нарисуйте по оси $p$ доверительный интервал. Насколько сильно в рамках этого интервала изменяется правдоподобие?
# берём картинку из 1 задачи и немного меняем её
[д] Слушайте, а что если $p = \frac{1}{2}$? Постройте тест отношения правдоподобий для проверки этой гипотезы. Проверьте её на уровне значимости $20\%$. Добавьте точку $p=\frac{1}{2}$ на картинку из предыдущего пункта. Лежит ли эта точка в доверительном интервале?
# исользуем результаты из пункта а)
[е] Усложняем ситуацию. В колонке sex
находится информация о том, мужчина перед нами или женщина. Пусть $p_w$ - вероятность того, что умрёт женщина. При этом $p_m$ - вероятность того, что умрёт мужчина. Выпишите функцию правдоподобия для такой модели. Найдите оценки максимального правдоподобия.
# надо не запутаться в выписывании функции
Если хотите проверить себя, можно посчитать среднюю выживаемость мужчин и женщин. Логично, что это и есть оценки максимального правдоподобия. Что мы видим, то и есть правда.
[ж] Насколько сильно отличаются $p_m$ и $p_w$? Правда ли, что функцию правдоподобия можно разбить на две и решать две отдельные задачи максимизации? Означает ли это, что оценки параметров независимы?
# код тут не нужен
[з] Последний этап. Давайте с помощью теста отношения правдоподобий проверим гипотезу о том, что $p_m = p_w$. Какое распределение будет у статистики?
# у нас уже есть всё необходимое в пунктах а),б) и е), надо только выписать статистику
Мораль: только что мы с вами попробовали оценить вероятность того, что человек умрёт на Титанике методом максимального правдоподобия. По ходу решения задачи мы поняли, что для мужчин и женщин эта вероятность очень разная. То есть вероятность смерти должна быть какой-то функцией от пола. Кроме пола, в таблице с данными есть огромное количество других признаков на любой вкус. Скорее всего, вероятность смерти — это функция ещё и от них.
Можно добавить их в нашу модель, описывающую вероятность смерти и построить более сложную штуку, которая называется логистической регрессией. В будущем вы обязательно займётесь этим. Оценивать эту модель вы будете ровно этим же самым методом максимального правдоподобия. В прошлом году я даже пытался рассказать про это, сохранилась тетрадка, но пара что-то не зашла. Видимо, я бежал впереди паровоза. Поэтому в этом году она на самостоятельном изучении для всех тех, кто тоже любит бегать от паровозов.
Каждый из группы второкурсников подкинул монетку два раза и никому не рассказывал, что у него выпало. В зависимости от того, что выпало на монетках, а также от того пробовал ли человек наркотики или нет, он сказал одну из двух фраз: очевидно или автомат у Козко.
? | пробовал | не пробовал |
---|---|---|
выпал ОО | автомат у Козко | очевидно |
другое | очевидно | автомат у Козко |
Ответы распределились следующим образом: Очевидно сказали $6$ человек, автомат у Козко сказали $8$ человек.
[a] Оцените методом максимального правдоподобия долю второкурсников, пробовавших наркотики. Сделайте все необходимые выводы на бумажке, после вбейте в R и найдите оценку.
# подсказок давать не буду, надо не запутаться :)
[б] Найдите $\hat{Var}(\hat p)$ и постройте $90\%$ доверительный интервал для доли пробовавших наркотики.
# ваш код
Мораль: иногда на практике собрать выборку довольно сложно. Мало того, что она может оказаться смещённой и с ошибками, так люди ещё и склонны врать. А ещё у людей куча стереотипов. Например, за ними наблюдается эффект якорения. Иногда, чтобы выборка на выходе оказалась хорошей, приходится придумывать разные сложные техники её сбора и разметки. Эксперимент, описанный выше помогает избежать вранья при ответе на неудобные вопросы. Всегда можно отмазаться, сказав, что монетка выпала не так. При этом в выборке будет какой-то сигнал, который можно восстановить по формуле полной вероятности.
Разметка данных очень сложный процесс, с которым сталкивается практически каждый аналитик. В интернете есть большое количество сервисов, которые позволяют за денюжку эту разметку провести. Один из таких сервисов, яндекс.толока. Например, на хабре есть описание разметки фоток Киркорова в толоке.
В табличке лежит информация о стоимости квартир в Москве и об основных параметрах этих квартир. На паре по доверительным интервалам мы уже работали с этой таблицей. Там мы предположили, что цены нормально распределены, и построили для них точные доверительные интервалы. Это было не очень корректно. Мы не проверили предпосылку о нормальности в данных, хотя обычно это нужно делать. На следующих парах мы посмотрим на то, как это делать.
df = read.csv('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/flat.csv', sep='\t')
nrow(df) # сколько всего строчек в табличке
head(df)
n | price | totsp | livesp | kitsp | dist | metrdist | walk | brick | floor | code |
---|---|---|---|---|---|---|---|---|---|---|
1 | 81 | 58 | 40 | 6 | 12.5 | 7 | 1 | 1 | 1 | 3 |
2 | 75 | 44 | 28 | 6 | 13.5 | 7 | 1 | 0 | 1 | 6 |
3 | 128 | 70 | 42 | 6 | 14.5 | 3 | 1 | 1 | 1 | 3 |
4 | 95 | 61 | 37 | 6 | 13.5 | 7 | 1 | 0 | 1 | 1 |
5 | 330 | 104 | 60 | 11 | 10.5 | 7 | 0 | 1 | 1 | 3 |
6 | 137 | 76 | 50 | 9 | 11.0 | 7 | 1 | 1 | 1 | 8 |
[а] Сразу прыгаем в пекло. Предположим, что цена на квартиру имеет нормальное распределение, $price \sim N(\mu, \sigma^2)$. Выпишите для этой задачки функцию правдоподобия, найдите оценки для параметров, а также оценку их ковариационной матрицы. Сделайте всё это руками на листочке.
# на листочке ручками
[б] Найдите оценки для $\mu$ и $\sigma^2$ на компьютере. Правда ли, что они асимптотически независимы? Правда ли, что они просто независимы (тут надо вспомнить вторую домашку, доказывать это опять не надо, хватит ответа да или нет). В чём разница между независимостью и асимптотической независимостью?
# можно переиспользовать код из 2 задачи
[в] Постройте для средней цены $90\%$ доверительный интервал.
# изи
[г] С чего всё началось? Сундук выплюнул на нас выборкe а мы выдвинули предпосылку, что она из нормального распределения. Постройте для цен гистограмму. Как думаете, верна ли эта предпосылка?
# изи
[д] Попробуем исправиться. Давайте возьмём от цен логарифм. Постройте для $\ln(price)$ гистограмму. Стало ли распределение более похожим на нормальное?
# совсем изи
[е] Иногда исследователи предполагают, что цены имеют логнормальное распределение. Если вы никогда не слышали про него, бегом на википедию! Обратите внимание на раздел Occurrence and applications, там есть примеры, в которых используют эту модель.
Есть случайная величина $price$. Для неё у нас есть $E(price) = \mu$ и $Var(price) = \sigma^2$. Есть случайная величина $\ln(price) \sim N(a, \tau^2)$. Параметры этих двух случайных величин взаимосвязаны между собой. И вы даже знаете как, если сходили на википедию. Методом максимального правдоподобия найдите оценку для параметров $a$ и $\tau^2$ (используйте функцию, вбитую в пункте б).
# примените старую функцию для ln(price)
[ж] Вспомните свойство инвариантности оценок максимального правдоподобия, состоящее в том, что если $\hat \theta$ это оценка для $\theta$, тогда $g(\hat \theta)$ это оценка для $g(\theta)$. Примените это свойство и получите оценки для $\mu$ и $\sigma^2$.
# вбить формулки
[з] Продолжаем вспоминать свойства оценок максимального правдоподобия. Оценки, полученные этим методом имеют асимптотически нормальное распределение. Интересующая нас оценка
$$\hat \mu = g(\hat a, \hat \tau^2)$$это функция от оценок максимального правдоподобия. Мы хотели бы построить для $\hat \mu$ доверительный интервал. Беда в том, что мы не знаем для него дисперсии.
Вспоминаем другую замечательную штуку, дельта-метод. Он говорил нам, что если
$$X \sim N(\mu; \sigma^2),$$тогда
$$g(X) \sim N(g(\mu); \sigma^2 \cdot (g'(\mu))^2).$$Как мы это получили? Просто разложили функцию $g$ в ряд Тэйлора. Все распределения были асимптотическими.
В нашей ситуации беда в том, что функция $g$ зависит сразу от двух асимптотически нормальных случайных величин. Нам нужно обобщить дельта-метод на случай функции от нескольких аргументов и раскладывать в ряд Тэйлора по двум аргументам. На паре мы этим занимались, но я не добавил разложение в конспект. Поэтому давайте ещё разок посмотрим на то, как это работает на конкретном примере.
Пущай у нас есть две оценки $\hat a$ и $\hat \tau^2$. И нас интересует $g(a, \tau^2) = \frac{\tau^2}{a}$. Если вы правильно сделали всё для поиска ковариационной матрицы, у вас должно было получится, что
$$ [-H]^{-1} = \begin{pmatrix} \frac{\tau^2}{n} & 0 \\ 0 & \frac{\tau^2}{2n} \end{pmatrix} $$Если мы подставим сюда оценку $\hat \tau^2$, то мы получим оценку ковариационной матрицы для вектора параметров.
Когда мы говорили про одномерный случай, мы просто умножали дисперсию на квадрат производной. Сейчас у нашей функции два параметра. Значит производных тоже две. Найдём её градиент.
$$ \nabla g = \begin{pmatrix} \frac{\partial g}{\partial a} \\ \frac{\partial g}{\partial \tau^2} \end{pmatrix} = \begin{pmatrix} \frac{-\tau^2}{a^2} \\ \frac{1}{a} \end{pmatrix} $$На паре мы немного поупражнялись в работе с квадратичными формами и доказали, что выражение $\nabla g^T \cdot [-H]^{-1} \cdot \nabla g$ это ровно то же самое, что умножение на квадрат производной в одномерном случае. Добиваем задачку:
Чтобы получить оценку для дисперсии, нужно подставить в формулу $\hat a$ и $\hat \tau^2$. Дальше можно стоить для $g$ доверительный интервал.
Подытожим. Если вектор $\hat \theta$ имеет асимптотическое многомерное нормальное распределение
$$ \hat \theta \sim N(\theta; Var(\theta)), $$и у нас есть дифференцируемая функция $g(\theta)$, тогда случайная величина
$$ g(\hat \theta) \sim N \left( g(\theta); \nabla g^T \cdot Var(\theta) \cdot \nabla g \right) $$Что я хочу от вас: Найдите распределение средней цены и постройте для неё $90\%$ доверительный интервал. Насколько координально он отличается от интервала, поулученного в пункте [в]?
# тут вот уже посложнее
Мораль: мы попробовали сходу замоделировать цены, отталкиваясь от предпосылки, что они распределены нормально. Мы напоролись на проблемы, связанные с тем, что наши предпосылки не выполнялись. Настоящий анализ данных - это борьба за предпосылки. Часто приходится довольно долго подбирать такую спецификацию модели, которая лучше ляжет на данные. Обычно перед тем, как ринуться в бой, проверяют серию из гипотез. Мы посмотрим на то, как это делается в следующие разы, когда резко погрузимся в финансовый мир.
Логарифмирование переменной — это довольно частый приём, с помощью которого у распределения пытаются обрубить длинный хвост и сгладить его. Благодаря этому сглаживанию модели, основанные на нормальности (например тот же самый МНК), начинают работать лучше. Есть целый набор нормализующих преобразований, который называется преобразованием Бокса-Кокса. На эконометрике вы будете его обсуждать. Наверное.
В этой задачке нет реальных данных, но можно ещё разок на своей шкуре прочувствовать, как работают методы моментов и максимального правдоподобия. Пусть $X_1, \ldots, X_n$ независимые случайные величины с плотностью распределения
$$ f_X(x) = \begin{cases} (a+1) \cdot x^a, \quad x \in [0;1] \\ 0, \text{ иначе} \end{cases} $$[а] Найдите руками для параметра $a$ оценку методом моментов. Найдите оценку методом максимального правдоподобия.
[б] Найдите асимптотические распределения этих оценок. Для максимального правдоподобия через гессиан, для метода моментов через дельта-метод.
# можно пользоваться компухтером
[в] Постройте на основе этих распределений $95\%$ доверительные интервалы. Для какой оценки доверительный интервал оказался короче?
# нужно пользоваться компухтером
[г] Какая из оценок является асимптотически более эффективной? Для ответа на этот вопрос найдите $$\lim_{n \to \infty} \frac{\hat Var(\hat a_{ML})}{\hat Var(\hat a_{MM})}.$$
Правда ли, что для более эффективной оценки доверительный интервал уже? Почему?
В датасете faithful
лежит информация о времени, которое проходит между извержениями гейзера "Старый служака".
data("faithful")
head(faithful)
eruptions | waiting |
---|---|
3.600 | 79 |
1.800 | 54 |
3.333 | 74 |
2.283 | 62 |
4.533 | 85 |
2.883 | 55 |
y = faithful$waiting
Давайте немного поисследуем это время. Филипп как-то раз на паре говорил, что время между разными событиями иногда моделируется с помощью экспоненциального распределения. Он даже показывал пример с данными по игре в Киллера!
Можно поверить ему, и немного повозиться этим, но мы будем умнее. К чёрту этого Филиппа с его чушью, давайте посмотрим на данные и сами примем решение!
[а] Постройти для времени между извержениями гистограмму. Что вы на ней видите?
# изи
Ну дела! Распределение времени для гейзера оказывается бимодальным. Википедия подсказывает, что так происходит из-за того, что вода входит в гейзер на нескольких различных уровнях. Получается, что время, которые мы получаем на выходе это результат смеси двух распределений (вспомните как вы генерировали смесь в первой домашке, тут за вас это сделала природа).
Задача разделения смесей — это классическая задача обучения без учителя из машинного обучения. И уже сейчас мы с вами можем попытаться её решить. Без учителя она, потому что мы не знаем какое извержение было из какого распределения и лишь можем догадываться об этом, исходя из полученных данных.
Давайте предположим, что с вероятностью $p$ данные приходят из нормального распределения со средним $\mu_1$ и дисперсией $\sigma_1^2$ и с вероятностью $1-p$ из нормального распределения со средним $\mu_2$ и дисперсией $\sigma_2^2$. В таком случае плотность распределения случайной величины будет выглядеть как:
$$ f(x) = p \cdot \frac{1}{\sqrt{2 \pi} \sigma_1} e^{- \frac{(x - \mu_1)^2}{2 \sigma_1^2}} + (1 - p) \cdot \frac{1}{\sqrt{2 \pi} \sigma_2} e^{- \frac{(x-\mu_2)^2}{2 \sigma_2^2}} $$Нам необходимо оценить пять параметров.
[б] Вбейте логорифмическую функцию правдоподобия.
# не запутайтесь с вбиванием вектора параметров
К сожалению, для такой функции правдоподобия не существует глобального максимума. Почему, можно попытаться понять самому, либо прочитать в книжке Лагутина на странице 117 в примере 8. Тем не менее, при большом числе наблюдений, почти наверное, есть локальный максимум, дающий хорошие оценки. Про это тоже есть на странице 117.
Пакет maxLik
почему-то начинает тупить на этом примере. Поэтому, давайте ради разнообразия воспользуемся другой функцией для оптимизации, встроенной в R. Она носит гордое имя nlm
и умеет только минимизировать. Давайте домножим наше правдоподобие на минус один.
# если надо, название поменяйте!
psevdo_lnL <- function(theta,y){
return(-1*mixt_lnL(theta, y))
}
Теперь проведём процедуру оптимизации. Параметр hessian = TRUE
отвечает за то, что на выходе посчитается гессиан. Когда вы запустите код ниже, на выходе вы получите кучу красного текста и неадекватные оценки на параметры.
result <- nlm(mixt_lnL, c(0.5, 1, 1, 1, 1), y, hessian=TRUE)
result
Это происходит из-за того, что мы плохо подобрали стартовые точки и не смогли пробиться в нужный нам локальный оптимум. Давайте попробуем подобрать стартовые точки на основе нашей интуиции более грамотно. Посмотрим на гистограмму, которую мы построили выше. Там где-то треть лежит в левом куполе. Берём эту цифру как стартовую для $p$.
[в] Для двух средних и дисперсий нам тоже нужны приближения. Придумайте как найти их и найдите.
# скорее всего, первая же пришедшая в голову идея - правильная
Теперь нужно провести оптимизацию. Подставьте в вектор c( )
свои приближения. Если вы их подобрали удачно, процедура оптимизации сойдётся куда нужно.
result <- nlm(psevdo_lnL, c( ), y, hessian=TRUE)
# Почти угадали оценки на глазок! А всё потому что метод максимального
# правдоподобия говорит: че видишь, то и есть на самом деле!
result$estimate
[г] Ура! Мы оценили параметры. Теперь давайте найдём их ковариационную матрицу. Наши оценки будут иметь асимптотически нормальное многомерное распределение. Вы же ещё не успели забыть, как искать её через матрицу Гессе?
# изи
[д] Интересно, насколько сильно различается среднее время извержения гейзера в разных состояниях. Постройте доверительный интервал для разности $\mu_2 - \mu_1$.
# изи
[е] Постройте на одно картинке гистограмму и оценённую теоретическую плотность распределения времени извержений.
# надо помучиться и что-нибудь получится
[ж] Гейзер может работать в двух режимах. Судя по полученным оценкам, оба режима по своему разбросу, в плане времени, оаботают одинаково. Проверьте с помощью теста отношения правдоподобий, гипотезу, что $\sigma_1 = \sigma_2$ против двухсторонней альтернативы.
# вы можете сделать это!
Мораль: «Жизнь как коробка шоколадных конфет: никогда не знаешь, какая начинка тебе попадётся». Когда оцениваете модели, всегда держите в голове то, что это делается численными методами. Для хорошей сходимости, для алгоритма не помешает хорошая инициализация. Когда вы дальше будете смотреть на машинное обучение и нейросетки, инициализация весов перед обучением будет занимать там одно из важнейших мест.
Кстати говоря, довольно часто задачи, где выскакивают очень сложные функции правдоподобия, решают с помощью EM— алгоритма. Но мы про него не будем говорить в нашем коротеньком курсе, хотя и могли бы. Но мы можем поговорить о нём на каком-нибудь ещё факультативе :)
Подробнее про задачу разделения смесей, EM и кластеризацию в R можно посмотреть тут. Написано не очень понятно, но впечатляет.
В группе мемы про машинное обучение для взрослых мужиков юных леди постят мемы про машинное обучение. Юные леди смотрят мемы про машинное обучение, лайкают их, комментируют и репостят. В итоге рождается табличка со статистикой.
df = read.csv('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data//memes_likes.csv', sep=';')
nrow(df)
head(df)
X | Date | likes | comments | views | reposts |
---|---|---|---|---|---|
0 | 2018-05-04 20:23:54 | 6 | 4 | 185 | 0 |
1 | 2018-05-02 01:20:28 | 13 | 9 | 309 | 0 |
2 | 2018-04-29 16:48:14 | 0 | 0 | 147 | 0 |
3 | 2018-04-22 17:03:31 | 3 | 0 | 326 | 0 |
4 | 2018-04-16 14:47:33 | 4 | 0 | 417 | 0 |
5 | 2018-04-12 17:37:35 | 9 | 0 | 504 | 1 |
Нам предстоит очень серьёзно поработать с этой табличкой. Важно: постарайтесь в плане кода быть в этой задачке полаконичнее.
[a] Постройте для числа лайков, числа комментов, числа репостов и числа просмотров гистограммы. Как вы думаете, какие именно распределения выплюнули из сундука на нас эти выборки?
# изи
Гипотезы о распределениях мы проверять пока не умеем, поэтому сразу же идём дальше.
[б] Окей. У нас есть $4$ разных случайных величины. Давайте попробуем сконструировать фабрику по производству моделей. Мы попробуем предположить для каждой случайной величины по очереди три вещи:
Давайте для всех трёх ситуаций вобьём логарифмические функции правдоподобия.
# первая функция правдоподобия
# вторая функция правдоподобия
# третья функция правдоподобия
[в] Предположим, что мы исследуем лайки. Мы пытаемся сначала описать их нормальным распределением, потом распределением Пуассона, потом геометрическим. Таким образом, мы получаем для лайков три модели. Возникает вопрос: какая из моделей лучше всего описывает выборку? Как выбрать конкретную?
Для ответа на этот вопрос нужно вспомнить что такое правдоподобие. Правдоподобие — это вероятность того, что мы получим именно нашу выборку при фиксированных параметрах. Значит, мы должны для каждой из трёх моделей посчитать значение правдоподобия в точке оптимума. Где оно больше, там больше вероятность рождения именно нашей выборки. Ту модель и выбираем.
Давайте возьмём лайки и оценим для них все три модели. Пакет maxLik
выводит в протоколе значение логарифма правдоподобия в точке оптимума. Вытащите это значение из всех трёх моделей и возьмите от него экспоненту. Это вероятность получить нашу выборку. Для какой модели правдоподобие оказалось самым большим?
# ТАКААААААЯ МАЛЕНЬКАЯ?!
[г] У схемы, которую мы описали в предыдущем пункте, есть проблема. Она заключается в том, что в нормальном распределении я оцениваю два параметра, в распределении Пуассона и геометрическом, один параметр. Чем больше параметров я оцениваю, тем больше у меня свободы и тем выше правдоподобие. Чтобы корректно сравнивать между собой эти две модели, нужно ввести какой-то штраф на число переменных.
Например, можно сделать это вот так ($k$ - число параметров):
$$ AIC = 2 \cdot k - 2 \cdot \ln L $$Тогда чем меньше величина $AIC$, тем лучше модель. Такой способ введения штрафа называют информационным критерием Акаике. Можно ввести штраф более извращённым способом. Например, вот так:
$$ SC = k \cdot \ln n - 2 \cdot \ln L. $$Такой способ штрафа более жёсткий. Он включает в себя не только число оцениваемых параметров $k$, но и число наблюдений $n$. Его называют критерием Шварца. Иногда говорят, что это байесовский информационный критерий, потому что его можно вывести с помощью формулы Байеса.
На эконометрике вы будете более пристально смотреть на эти критерии. Почему-то люди не понимают, что это просто оштрафованное правдоподобие. Поэтому я решил познакомить вас с ними уже сейчас.
Возьмите все три модели для лайков и найдите для них значение критерия Акаике. Сделать это можно командой AIC(result)
. Какая модель оказывается лучшей? Согласуется ли это с тем, что вы решили в пункте [а], посмотрев на гистограмму?
# не пугайтесь этого пункта!
[д] Оцените все три модели для репостов. Какая оказывается лучшей? Скорее всего, это оказалось геометрическое распределение. Проинтерпретируйте параметр $p$ для него.
# вот тут стоит быть лаконичным
[е] Оцените все три модели для просмотров. Какая оказалась лучшей? Согласуется ли это с вашими ощущениями из пункта [а]?
# снова лаконичность
[ж] Проделайте эту процедуру в последний раз, для комментариев.
# и ещё разок лаконичность
Повышаем планочку. Давайте посмотрим на гистограмму для числа комментариев. В нуле у них довольно большой пик, который выбивается из распределения Пуасона.
qplot(df$comments)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Такой пик в нуле выскакивает довольно часто. Даже специальную модель для борьбы с ним придумали. Она называется zero inflated model. Давайте попробуем её на вкус. Правда сначала её надо построить.
Делай раз: мы хотим, чтобы распределение Пуассона для нас работало, начиная с $X = 1$. В таком случае нам его надо сдвинуть вправо так, чтобы сумма вероятностей по-прежнему оставалась равна единице. Мы знаем, что
$$ \sum_{k=0}^{\infty} P(X = k) = \sum_{k=0}^{\infty} \frac{\lambda^k e^{-\lambda}}{k!} = e^{-\lambda} + \sum_{k=1}^{\infty} \frac{\lambda^k e^{-\lambda}}{k!} = 1. $$Если мы решим оставить только сумму, начиная с единицы, получится, что
$$ \sum_{k=1}^{\infty} \frac{\lambda^k e^{-\lambda}}{k!} = 1 - e^{-\lambda}. $$Чтобы перед нами было полноценное распределение и все вероятности в сумме давали $1$, нам надо поделить сумму слева на $1 - e^{-\lambda}$. Получается, что для распределения Пуассона, обрезанного со стороны нуля, формула для поиска вероятности выглядить как
$$ P(X = k) = \frac{1}{1 - e^{-\lambda}} \cdot \frac{\lambda^k e^{-\lambda}}{k!}. $$Можно получить эту формулу исходя не из интуции, а из формулы условной вероятности:
$$ P(X = k \mid X > 0) = \frac{P(X = k \cap X > 0)}{P(X > 0)} = \frac{\frac{\lambda^k e^{-\lambda}}{k!}}{1 - e^{-\lambda}} $$Делай два: Теперь давайте построим смесь из двух распределений. Случайная величина $X$ будет принимать с вероятностью $p$ значение $0$, и с вероятностью $1 - p$ будет распределена по Пуассону со сдвигом:
$$ \begin{aligned} & P(X = 0) = p \\ & P(X = k) = (1 - p) \cdot \frac{1}{1 - e^{-\lambda}} \cdot \frac{\lambda^k e^{-\lambda}}{k!}. \end{aligned} $$Построенная модель — это ещё не совсем zero inflated model. У неё есть кое-какой минус.
[з] Выведите для такой модели функцию правдоподобия. Искать её максимум придётся по двум параметрам: $\lambda$ и $p$. Вбейте её и найдите оценки.
# тут посложнее, но терпимо
[и] Посчитайте для оценённой модели критерий Акаике и сравните её с лучшей из предыдущих, построенных для комментариев, трёх моделей.
# одна строчка
У такой формулировки модели есть минус. Невозможно проверить гипотезу о том, что в пике нет никаких особенностей. Если $p = 0$, то у нас просто-напросто не бывает нулевых значений. Хочется, чтобы у нас была возможность протестировать такую гипотезу. Для этого ноль выносится в отдельную категорию не в результате обрезания распределения Пуассона, а немного иначе. Давайте домножим $P(X = k)$ на $(1-p)$, а потом просто вынесем $(1 - p) \cdot P(X = 0)$ в отдельное слагаемое. И тогда получится моделька:
$$ \begin{aligned} & P(X = 0) = p + (1 - p) \cdot e^{-\lambda} \\ & P(X = k) = (1 - p) \cdot \frac{\lambda^k e^{-\lambda}}{k!}. \end{aligned} $$Что мы видим здесь? Если $p=0$, то у нас получается распределение Пуассона. У нас возникает возможность проверить гипотезу $H_0$, состояшую в том, что в нуле действительно есть особенности.
[к] В пункте [з] параметр $p$ можно было проинтерпретировать, как вероятность возникновения нуля. Как он будет интерпретироваться в такой модели?
Ответ:
[л] Вбейте функцию правдоподобия для zero inflated model. Оцените значения параметров. Правда ли, что по AIC эта модель никак не отличается от обрезанного Пуассона?
# если вы справились с тем, что было выше, то и это не проблема
[м] Последний пункт! Нам осталось проверить гипотезу о том, что $p = 0$ в нашей zero inflated model. Сделаем это с помощью теста отношения правдоподобий на уровне значимости $5\%$. Если гипотеза о том, что $p = 0$ не отвергается, значит усложнять модель нет никакого смысла, так как поведение случайной величины в нуле ничем не отличается от Пуассоновского.
# да и это тоже не проблема
Мораль: Во-первых, можно брать какие-то уже знакомые вам модели и дорабатывать их до новых. Метод максимального правдоподобия позволяет проверить гипотезу о том, что доработки имели смысл. Не забывайте о том, что чем больше параметров вы оцениваете, тем более сложной получается модель, и тем лучше она описывает данные. Ситуация, когда мы берём слишком подробную модель для маленького объёма данных, и вылизываем их, называется переобучением. Это очень большая головная боль, с которой вы в будущем будете бороться, как с самой большой чумой 21 века.
Во-вторых, можно оценивать много моделей, а после сравнивать их между собой. Для этого можно использовать не только правдоподобие и информационные критерии, но и разные другие метрики. Но об этом тоже позже.
Впервые европейцы появились на берегах Австралии в $1606$ году. Голландцы под руководством Виллема Янсзона сошли на её прекрасные берега. Приплыли они в пустыню. Поэтому сразу же взошли на борт и свалили. В $1788$ году к Австралии пришвартовались несколько Британских кораблей. Пришвартовались они с правильной стороны материка. Это позволило основать в Австралии колонию для осуждённых, которая впоследствии превратилась в Сидней. С тех пор Австралия использовалась англичанами как тюрьма.
Со временем на материке появились и другие, никак не связанные с преступниками колонии. В $1901$ году Австралийские колонии создали свою собственную федерацию. На $2014$ год население Австралии составляет $23$ млн. человек (население Москвы - около $15$ млн.).
В $1858$ году колонисты привезли в Австралию $24$ кролика. Год спустя этих длинноухих зверьков можно было увидеть за $100$ километров от места высадки — и на север, и на запад. А через $3-4$ года число кроликов возросло до $750$ миллионов. И начался хаос.
Кролики на пастбищах начали соперничать с овцами и коровами. Десять кроликов съедали столько травы, сколько съедала одна овца. Только вот от овцы можно получить мяса в три раза больше. Фермеры уже подсчитывали свои убытки и объявили кроликам настоящую войну. Более того, кролики начали вытеснять из привычных зон обитания местные уникальные виды. К борьбе с кроликами подключились и ученые.
Сама борьба с кроликами принесла немало бед для Австралийской флоры и фауны. Первоначально решили завезти естественных врагов кроликов – лисиц, хорьков, кошек, горностаев, ласок. Но попытка не увенчалась успехом. Привезенные виды также стали активно размножаться. Вместо того, чтобы охотиться на быстрых кроликов, они стали поедать местных сумчатых животных и птиц, которые были не так быстры и не могли сопротивляться новым хищникам.
Тогда обратились к традиционным методам – ядохимикаты, отстрел, взрывание нор. Это было не эффективно, учитывая огромную численность животных. В штате Западная Австралия в период с $1901$ по $1907$ гг. построили огромный проволочный забор. Он так и называется – «Забор от кроликов №1». Забор постоянно патрулируется на машинах, кроличьи подкопы засыпаются, кролики отстреливаются.
Сначала забор патрулировали на верблюдах, привезённых из Африки (ничему жизнь не учит австралийских колонистов). После появления автомашин, верблюдов за ненадобностью выпустили на волю, они расплодились, стали уничтожать пастбища, и в Австралии появилась новая проблема.
В середине $50$-х гг. $20$ века для борьбы с кроликами стали использовать достижения медицины. В Австралию привезли кроличьих блох и комаров, зараженных вирусом миксоматоза. Это заболевание вызывает опухоли и смерть кроликов. Таким образом было уничтожено около $90\%$ заболевших животных. Но оставшиеся кролики выработали иммунитет, с течением времени стали реже заболевать и ещё реже умирать. Так что на данный момент проблема кроликов в Австралии до сих пор не решена.
В этой задачке мы попробуем немного поработать с результатами исследований $1956$ года. В нём Феннер и его коллеги изучали разные вирусы, из-за которых умирают кролики. Табличка с собранными данными есть внутри R, в пакете emdbook
. Установим пакет и подгрузим данные.
#install.packages('emdbook')
library('emdbook')
# подгрузили данные
data(MyxoTiter_sum)
df = MyxoTiter_sum #переименовал для удобства
# глянули на них
head(df)
grade | day | titer |
---|---|---|
1 | 2 | 5.207 |
1 | 2 | 5.734 |
1 | 2 | 6.613 |
1 | 3 | 5.997 |
1 | 3 | 6.612 |
1 | 3 | 6.810 |
В колонке grade
находится штамм вируса titr. Всего их в исследовании $5$ разновидностей. Пока что мы сосредоточимся на первом. В колонке titer
находится показатель концентрации вируса в коже кролика. Измерение делалось в разные дни после заражения. Мы будем игнорировать этот факт и притворимся, что все измерения пришли к нам из одного распределения, которое не зависит от времени. То есть, мы предположим, что выборка независима и одинакого распределена (на самом деле это неправда и чуть позже мы откажемся от такого подхода).
Сконцентрируем внимание на первом штамме вируса.
df = df[df['grade'] == 1, ]
Будем препдполагать, что данные о концентрации вируса пришли к нам из гамма-распределения. Из того самого распределения из ваших контрольных! Почему именно оно? Концентрация принимает положительные значения в рамках какого-то фиксированного промежутка, а дальше начинается хвост распределения. Надо только понять где этот промежуток и насколько резко он обрывается. Параметры гамма-распределения помогают сделать это.
[а] Посмотрите на плотность гамма-распределения на Википедии. Попытайтесь проинтерпретировать за что отвечают его параметры.
Экспоненциальное распределение, кстати говоря частный случай гамма-распределения. Давайте в этом убедимся.Выпишите плотность гамма-распределения.
Ответ:
[б] Методом максимального правдоподобия оцените параметры гамма-распределения. Если есть проблемы с вбиванием функции правдоподобия, воспользуйтесь приёмом с dgamma
из лекции (мы так хитрили при авариях на шахтах с Пуассоном).
# КАКОЙ ТАКОЙ ПРИЁМ?!
[в] На одном графике постройте гистограмму и оценённую в ходе метода максимального правдоподобия плотность распределения.
# снова мучения с картинками
Пришло время отказаться от предпосылки, что наши наблюдения по кроликам не зависят от времени. Это гнусная ложь. Усложняем модель! Концентрация вируса изменяется в зависимости от количества дней, прошедших после заражения.
[г] О, мой любитель графиков! Построй картинку, на которой будет видна зависимость концентрации вируса от времени (qplot
поможет). Какой характер носит эта зависимость?
# одна строка
Попробуем побаловаться с моделями времени жизни. Так можно замоделировать динамику различных популяций, объём дневных осадков. А ещё экономический рост. Простейшим примером модели времени жизни является модель Рикера. Давайте попробуем использовать её, чтобы промодлировать зависимость концентрации вируса в кроликах от времени.
\begin{equation*} \begin{aligned} & m = a \cdot t \cdot e^{-b \cdot t} \\ & X \sim Gamma(shape = s, scale = m) \\ \end{aligned} \end{equation*}Будем снова использовать гамма-распределение, чтобы описывать концентрацию вируса $X$ в момент $t$. Параметр scale
отвечает за интенсивность вируса в кролике. Мы хотим, чтобы он зависил от времени. Пусть вирусв крови кролика копится экспоненциально, то есть ведёт себя как функция $f(t) = a \cdot t \cdot e^{-b \cdot t}$.
Чтобы нормально справится с этим пунктом, загляните в семинар, в часть, где мы оценивали похожую модель для игры в Киллера.
[д] Методом максимального правдоподобия оцените параметры $a, b$ и $s$.
# вы же поняли, что делать? очень похоже на задачку про киллера с семинара
Мораль: Зайку жалко :(