Байесовская статистика и вероятностное программирование

В сегодняшнем занятия будут использоваться новая библиотека — pymc3.

Следует отметить, что сейчас (февраль 2019 года) активно идёт разработка новой, четвёртой, версии библиотеки PyMC, API которой претерпел значительные изменения из-за перехода с Theano на TensorFlow, которые теперь используется для всех вычислений. По этой причине данный туториал, вероятно, потеряет актуальность через некоторое время.

Введение в байесовский статистический вывод

Если данные говорят с вами, значит вы — байесовец.

Формально, теорема Байеса – это просто математическая формула. Однако её значение гораздо глубже. Теорема Байеса подводит нас к тому, что необходимо иначе взглянуть на процесс выдвижения и проверки идей.

С XVIII в. существует две основных философских позиции относительно того, как применять вероятностные модели. Первая определяет вероятность как нечто, заданное внешним миром. Вторая утверждает, что вероятность существует в головах людей. В русле первого подхода возникли вначале классическая и затем развивающая её частотная концепции вероятности. Второй подход нашёл выражение в концепции байесовской вероятности.

Сторонники классического подхода исходят из того, что истинные параметры модели не случайны, а аппроксимирующие их оценки случайны, поскольку они являются функциями наблюдений, содержащих случайный элемент. Параметры модели считаются не случайными из-за того, что классическое определение вероятности исходит из предположения равновозможности как объективного свойства изучаемых явлений, основанного на их реальной симметрии. На такое представление о вероятности повлияло то, что в начале своего развития теория вероятности применялась прежде всего для анализа азартных игр. Суждение вида «вероятность выпадения шестёрки при бросании игрального кубика равняется 1/6» основывается на том, что любая из шести граней при подбрасывании на удачу не имеет реальных преимуществ перед другими, и это не подлежит формальному определению. Таким образом, вероятностью случайного события $A$ в её классическом понимании будет называться отношение числа несовместимых (не могущих произойти одновременно) и равновозможных элементарных событий $m$ к числу всех возможных элементарных событий $n$: $$ P(A)=\frac{m}{n} $$

Однако такое определение наталкивается на некоторые непреодолимые препятствия, связанные с тем, что не все явления подчиняются принципу симметрии. Например, из соображений симметрии невозможно определить вероятность наступления дождливой погоды. Для преодоления подобных трудностей был предложен статистический или частотный способ приближенной оценки неизвестной вероятности случайного события, основанный на длительном наблюдении над проявлением или не проявлением события $A$ при большом числе независимых испытаний и поиске устойчивых закономерностей числа проявлений этого события. Если в результате достаточно многочисленных наблюдений замечено, что частота события $A$ колеблется около некоторой постоянной, то мы скажем, что это событие имеет вероятность. Данный тип вероятности был выражен Р. Мизесом в следующей математической формуле: $$ p=\lim_{x\to\infty}\frac{\mu}{n}, $$ где $\mu$ — количество успешных испытаний, $n$ — количество всех испытаний. Вероятность здесь понимается как частота успешных исходов я является чисто объективной мерой, поскольку зависит лишь от точного подсчёта отношения количества успешных и неуспешных событий.

Основываясь на этом подходе, статистика занималась созданием вероятностных моделей, которые включали в себя параметры, которые, как предполагалось, связаны с характеристиками исследуемой выборки. Параметры никогда не могут быть известны с абсолютной точностью до тех пор, пока мы не исследуем всю генеральную совокупность. До тех пор всегда существует вероятность отклонить гипотезу, когда она на самом деле верна, т. е. совершить ошибку первого рода. Для обозначения вероятности такой ошибки частотники используют понятие уровня значимости $\alpha$. Именно вероятность ошибки первого рода частотники ставят во главу анализа, определяя вероятность события. После каждого своего утверждения они обычно добавляют «... на доверительном уровне в 95%», подразумевая, что исследователь допускает вероятность ошибки в пяти процентах случаев (при $\alpha = 0,05$).

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

Частотный подход доминировал в XX веке, придя на смену другому пониманию вероятности, связанном с именем английского математика Томаса Байеса. Сущность байесовского подхода составляют три элемента: априорная вероятность, исходные статистические данных данные, постаприорная вероятность.

Байесовская статистика начинает построение своей модели при помощи понятия априорной вероятности, с помощью которой описывается текущее состояние наших знаний, относительно параметров распределения. Априорная вероятность, таким образом, — это степень нашей уверенности в том, что исследуемый параметр примет то или иное значение ещё до начала сбора исходных статистических данных. На этом основании байесовское понимание вероятности относят к группе субъективистских трактовок вероятности. Чаще всего предполагается, что для оценки степени уверенности необходимо привлечь экспертов, чьё субъективное свидетельство позволит избежать действительной многократной реализации интересующего нас эксперимента.

Следующий элемент — это исходные статистические данные. По мере их поступления статистик пересчитывает распределение вероятностей анализируемого параметра, переходя от априорного распределения к апостериорному, используя для этого формулу Байеса: $$ P(\theta|y)=\frac{P(y|\theta)P(\theta)}{P(y)}, $$ где

• $P(\theta)$ —— априорная вероятность гипотезы $\theta$, выражает наши знания о значениях параметров, предшествующие построению модели.
• $P(y|\theta)$ —— апостериорная вероятность, представляет собой распределение вероятностей для параметра в нашей модели, рассчитанное с учётом априорных знаний и правдоподобия новых.
• $P(\theta|y)$ —— правдоподобие, вероятность наблюдать данные $y$ при истинности гипотезы $\theta$
• $P(y)$ —— свидетельство (evidence, marginal likelihood), полная вероятность гипотезы $\theta$. Другими словами — вероятность наблюдения данных, усредненая по всем возможным значениям, которые могут принимать параметры. Можно рассматривать её как плотность распределения всех объектов, но, по сути, она нам не очень нужна: выбор максимальной вероятности принадлежности к классам не зависит от этого значения. Можно расценивать этот член как некоторый нормирующий коэффициент. Без него формула станет даже проще:

$$p(\theta|y)∝ p(y|\theta)p(\theta)$$

Суть формулы в том, что она позволяет переставить причину и следствие: по известному факту события вычислить вероятность того, что оно было вызвано данной причиной. Эту формулу также называют формулой обратной вероятности. Процесс пересмотра вероятностей, связанных с высказываниями, по мере поступления новой информации составляет существо машинного обучения и является одним из возможных способов формализации и операционализации следующего тезиса: «степень нашей разумной уверенности в некотором утверждении (касающемся, например, неизвестного численного значения интересующего нас параметра) возрастает и корректируется по мере пополнения имеющейся у нас информации относительно исследуемого явления» (Айвазян 2008). В частотном подходе данный тезис интерпретируется в свойстве состоятельности оценки неизвестного параметра: чем больше объём выборки, на основании которой мы строим свою оценку, тем большей информацией об этом параметре мы располагаем и тем ближе к истине наше заключение. Специфика байесовского подхода к интерпретации этого тезиса основана на том, что вероятность, понимаемая как количественное значение степени разумной уверенности в справедливости некоторого утверждения, пересматривается по мере изменения информации, касающейся этого утверждения. Поэтому в данном подходе вероятность всегда есть условная вероятность, при условии нынешнего состояния информации (в русле классического подхода исследователь скорее склонен рассматривать совместную вероятность.

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

Так или иначе, в 1980-х годах, стало ясно, что частотный подход к статистическому выводу не достаточно хорошо подходит для анализа нелинейных отношений в больших объёмах данных, производимых сложными системами при моделировании процессов реального мира. Для преодоления этих ограничений частотники создали нелинейные версии параметрических методов, такие как множественный нелинейный регрессионный анализ.

В то время как в частотном подходе происходили изменения, немногочисленные сторонники байесовского подхода упрямо продвигали свою точку зрения на модель статистического вывода. Как оказалось, байесовская модель лучше подходит для поиска ответов на некоторые практические вопросы, поскольку полнее учитывает прошлую информацию и располагает к предсказаниям. Например, намного важнее минимизировать вероятность ложноотрицательного диагностирования некоторой опухоли как раковой, чем вероятность её ложноположительного определения (ошибка первого рода).

Проведение тестирования на статистическую значимость оценивает лишь вероятность получения похожего результата с другим набором данных при сохранении тех же самых условий. Однако оно предоставляет ограниченную картину такой вероятности, поскольку в расчёт принимается ограниченное количество информации относительно исследуемых данных. И оно само по себе не способно вам сказать, являются ли основные положения исследования верными и будут ли подтверждены полученные результаты в различных условиях. Уровень $p$ говорит только о вероятности получения результата при (обычно) совершенно нереалистичных условиях нулевой гипотезы. А это совсем не то, что мы хотим узнать, — обычно мы хотим знать величину эффекта независимой переменной с учётом имеющихся данных. Это байесовский вопрос, а не частотный. Вместо этого значение $p$ часто интерпретируется так, будто бы оно показывало силу ассоциации.

С другой стороны, и у байесовского метода имеются несколько недостатков. Одним из них является необходимость привлекать для расчёта априорные данные, которые могут быть недоступны. А если они и доступны, то, как отмечалось выше, часто носят субъективный характер. Другой недостаток — сложность вычислений. В вышеописанном примере для вычисления байесовской вероятности нам необходимо было вычислить частотную вероятность, полную вероятность, и, наконец, собственно байесовскую вероятность. Сложность байесовских вычислений частично объясняет тот факт, что байесовские методы вновь обрели популярность с развитием вычислительной техники. Следующий недостаток байесовского метода — неинтуитивность, непонятность его результатов для обыденного сознания. Именно на этой неинтуитивности построен знаменитый парадокс Монти Холла, который легко решает с помощью формулы Байеса.

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

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

Все неопределенности, которые мы имеем, моделируются с использованием вероятностей. В других статистических парадигмах существуют разные типы неизвестных величин; в байесовских рамках все, что неизвестно, рассматривается одинаково. Если мы не знаем количество, мы назначаем ему распределение вероятностей. Затем теорема Байеса используется для преобразования предшествующего распределения вероятности $p (θ)$ (что мы знаем о данной проблеме до наблюдения данных) в апостериорное распределение $p (θ | D)$ (то, что мы знаем после наблюдения данных). Другими словами, байесовская статистика является формой обучения.

Статистический вывод

• Доверительный интервал
  • Байесовский: credible interval
  • Частотный: confidence interval
• Результат
  • Байесовский: Оценка уверенности в виде вероятности, предсказание
  • Частотный: Вердикт о правильности или неправильности гипотезы при заданном уровне доверительного интервала, уровне значимости
• Выборка
  • Частотный: Фиксирована
  • Байесовский: Имеет некоторое распределение
• Оцениваемый параметр
  • Байесовский: Случайная величина с некоторым заданным априорным распределением
  • Частотный: Фиксирован

Вероятностное программирование

Возможность автоматизации части логического вывода привела к разработке вероятностных языков программирования (PPL), которые позволяют четко разделить создание модели и логический вывод. Использую PPL пользователи могут задать вероятностную модель, написав несколько строк кода, после чего вывод делается автоматически.

Бросание монетки — частотный подход

В большинстве учебников по теории вероятностей изучение начинается с бросания монетки. То же самое и в байесовской статистики, только тут всё куда интересней. Пусть кто-то бросает монету и мы наблюдаем исходы: 1 — выпал орёл, 0 — выпала решка. Наша задача оценить вероятность выпадения орла (т.е. $\theta$ в данном случае — это количество выпадений орла). Если любите формулы и интегралы, читайте объяснение Николенко или Дьяконова, а ещё лучше целый курс.

Эта случайная величина — результаты бросания монетки — подчиняется распределению Бернулли. Мы также знаем, что последовательность независимых случайных величин $Y=X_{1}+X_{2}+\ldots +X_{n}$, имеющих одинаковое распределение Бернулли, имеет биномиальное распределение с параметрами $n$ и $p$, где $p$ — вероятность успеха в $n$ испытаниях. Это записывается в виде: $Y\sim {\mathrm {Bin}}(n,p).$

Функция вероятности данной случайной величины задаётся формулой Бернулли:

$$P_{n}^{k}=C_{n}^{k}p^{k}(1-p)^{n-k},$$ где $k$ — количество успешных исходов.

Давайте запрограммируем формулу на Python.

Формула Бернулли позволяет ответить на вопрос, какова вероятность в 10 бросках монеты получить 9 «орлов», если монека честная (вероятность «орла» составляет 50%).

Как видно, вероятность достаточно небольшая.

from math import factorial

def num_of_successes(n, k):
    return factorial(n)/(factorial(k) * factorial(n - k))

def probability_of_success(p, n, k):
    C_kn = num_of_successes(n, k)
    return C_kn * (p**k) * (1 - p)**(n - k)


def successes_not_less_k(p, n, k):
    return sum([probability_of_success(p, n, k) for k in range(k + 1)])

successes_not_less_k(p=5/6, n=10, k=4)

Вероятность успешных испытаний для бросков монетки моделируется при помощи биноминального распределения. Это дискретное распределение показывает вероятность успеха в серии из N испытаний при фиксированном значении θ, если все испытания независимы друг от друга. В терминах байесовской статистики это будет правдоподобие $p(y|\theta)$. Следующий код генерирует 9 биномиальных распределений; у каждого есть своя легенда, где показаны параметры $n$ — количество испытаний и p—вероятность успеха в каждом испытании.

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

В нашем примере в роли априорного распределения будет выступать бета-распределение. Оно выражает нашу веру в то, что $\theta$ принимает именно такое значение, какое декларируется. Мы выбрали бета-распределение, поскольку оно является сопряжённым к Биномиальному распределению, т.е. может принимать тот же вид, но при других параметрах. Также, если функцию правдоподобия для распределения Бернулли умножить на плотность бета-распределения, то мы снова получим бета-распределение. Таким образом, переходя от априорного распределения к апостериорному вид распределения не меняется, что удобно.

Здесь изображено то, что нас интересует, то есть, распределение вероятностей для значений $\theta$ после учёта имеющихся данных, принимая во внимание априорную информацию. Эти распределения удалось легко вычислить, потому что биноминальное и бета распределения сопряжённые. Предположим, однако, что априорное распределение не является сопряжённым и мы не можем решить задачу, так сказать, вручную. На практике обычно бывает именно так.

Тут на помощь приходит семплирование методом MCMC.

MCMC

MCMC генерирует сэмплы из апостериорного распределения вероятностей, создавая обратимую цепь Маркова, равновесное распределение которой – это целевое апостериорное распределение [6].

1. Задаётся начальное значение параметра
2. Перемещаемся из этой позиции в какое-то другое место (каким-либо образом меняем значение параметра).
3. Если полученное распределение с новым параметром объясняет данные лучше, чем при предыдущем параметре (т.е. вероятность данных выше), то двигаемся в этом направлении дальше.

Обычно, движение осуществляется согласно алгоритму Метрополиса — Гастингса.

PyMC3

PyMC3 позволяет выводить апостериорное распределение без непосредственного вывода формул.

$$θ ∼ Beta (α , β )$$ $$y ∼ Bin(n = 1, p = θ )$$

The function summary provides a text-format summary of the posterior. We get the mean, standard deviation, and the HPD intervals:

Since we are approximating the posterior with a finite number of samples, the first thing we need to do is to check if we have a reasonable approximation. There are several tests we can run, some are visual and some quantitative. These tests try to find problems with our samples but they cannot prove we have the correct distribution; they can only provide evidence that the sample seems reasonable. If we find problems with the sample, the solutions are:

• Increase the number of samples.
• Remove a number of samples from the beginning of the trace. This is know as burn-in.
• Re-parametrize the model, that is express the model in a different but equivalent way.
• Transform the data.

The left side shows our marginal posterior – for each parameter value on the x-axis we get a probability on the y-axis that tells us how likely that parameter value is. KDE plots should look like smooth curves. Often, as the number of data increases, the distribution of each parameter will tend to become Gaussian-like; this is due to the law of the large numbers. Of course, this is not always true.

The plot on the right should look like white noise; we are looking for good mixing. This plot shows individual sampled values at each step during the sampling. We should not see any recognizable pattern, we should not see a curve going up or down, instead we want a curve meandering around a single value.

A quantitative way to check for convergence is by using the Gelman-Rubin test. The idea of this test is to compare the variance between chains with the variance within chains, so of course we need more than one chain for this test to work. Ideally, we should expect a value of $R=1$. As an empirical rule, we will be ok with a value < 1.1; higher values are signaling a lack of convergence:

We can also visualize the $R$ for every parameter together with the mean, 50% HPD (Highest Posterior Density) and 95% HPD for each parameter distribution using the function forestplot. An HPD is the shortest interval containing a given portion of the probability density.

An ideal sample will lack autocorrelation, that is, a value at one point should be independent of the values at other points.

Assignment

Using PyMC3, change the parameters of the prior beta distribution to match those of the previous chapter and compare the results to the previous chapter. Replace the beta distribution with a uniform one in the interval [0,1]. Are the results equivalent to beta(α =1,β =1)? Is the sampling slower, faster, or the same? What about using a larger interval such as [-1, 2]?

Байесовский подход для сравнения групп

Одной из частных задач в анализе данных является сравнение групп. Пример такой задачи можно посмотреть в секции «Анализ данных для маркетологов», где в последней задаче надо было сравнить выручку, полученную с пользователей Айфонов и Айпадов. Там это было сделано при помощи двувыборочного t-test'а, который является типичным методом частотного подхода. Суть этого теста в том, чтобы проверить, значимо ли различаются средние значения в двух выборках, а его результатом является вероятность того, что средние различаются незначимо (нулевая гипотеза). Если эта вероятность меньше определённого значения (p-value), обычно $0.05$, мы можем говорить, что с вероятностью не менее 95% средние величины в двух выборках различаются значимо.

Однако, полагаться только на статистическую значимость может быть проблематично по причине того, что идея статистической значимости связана с вычислением p-значений. Существует множество исследований и эссе, показывающих, что чаще всего p-значения используются и интерпретируются неправильно даже очень компетентными учёными. В рамках байесовского вывода нам не нужно вычислять p-значения. Вместо этого мы собираемся сосредоточиться на оценке, две группы действительно различны. В конце концов, на практике мы чаще всего хотим знать величину эффекта, то есть количественную оценку силы изучаемого явления.

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

Модель для этой проблемы почти такая же, как и раньше; единственное отличие состоит в том, что теперь параметры µ и σ будут векторами вместо скалярных случайных величин. Другими словами, каждый раз, когда мы производим выборку из наших априорных значений, мы получим четыре значения для μ и четыре значения σ вместо одного, как в нашей предыдущей модели. При помощи PyMC3 это очень удобно воплотить: для задания арпиорного распределения нам нужно передать аргумент shape, а для вероятности нам нужно правильно индексировать μ и σ, и поэтому мы создаем переменную idx:

Оценим различие между группами при помощи Cohen's d, который рассчитывается как разница между средними, деленными корень из дредних квадратов отклонений и очень популярен, например, в психологии, и Probability of Superiority вероятности того, что показатель одной случайной выборки больше, чем показатель другой случайной выборки из той же генеральной совокупности.

Чтобы интерпретировать эти графики, надо сравнить значение 0 (нет разницы) с интервалом HPD. Согласно предыдущему рисунку, у нас есть только один случай, когда 95% HPD исключают 0, разницу в чаевых между четвергом и воскресеньем. Для всех остальных примеров мы не можем исключить разницу в 0 (в соответствии с критериями HPD). Но даже для этого случая средняя разница составляет ≈0,5 доллара.

Линейная регрессия в байесовском стиле

Обычно линейная регрессия выглядит так: $$Y = X\beta + \epsilon$$

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

$$Y \sim \mathcal{N}(X \beta, \sigma^2)$$

Таким образом, мы рассматриваем $Y$ как случайную переменную (или случайный вектор), элементы которой распределены нормально. Среднее значение этого нормального распределения обеспечивается нашим линейным предиктором с дисперсией $\sigma^2$.

Хотя это, по сути, одна и та же модель, есть два критических преимущества байесовского подхода:

• Наличие априорного распределения: Мы можем количественно оценить любые предварительные знания, которые у нас могут быть, выбрав соответствующие априорные распределения. Например, если мы думаем, что дисперсия, вероятно, будет небольшой, мы бы выбрали априоное распределение, где вероятность небольших значений выше.
• Количественная оценка неопределенности: мы не получаем единственную оценку 𝛽, как указано выше, но вместо этого получаем полное апостериорное распределение о том, насколько вероятны различные значения. Например, с небольшим количеством данных наша неопределенность для 𝛽 будет очень высокой, и мы получим очень широкие постериорные распределения для для этого параметра.

Сгенерируем 200 точек, которые описываются уравнением $y=2x+1+\epsilon$, и нанесём их и прямую на график.

In the GLM we thus do not only have one best fitting regression line, but many. A posterior predictive plot takes multiple samples from the posterior (intercepts and slopes) and plots a regression line for each of them. Here we are using the plot_posterior_predictive_glm() convenience function for this.

Многоуровневые регрессионные модели в байесовском стиле

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

Чтобы подчеркнуть эффект иерархической линейной регрессии, мы сначала построим неиерархическую байесовскую модель, в которой для каждого округа будет задан свой собственный свободный коэффициент. $$radon_{i, c} = \alpha_{c} + \beta_{c}*\text{floor}_{i, c} + \epsilon_c$$ Фактически, это будет просто набор отдельных моделей. Поскольку у нас нет предварительной информации о том, каким может быть перехват или регрессия, мы будем использовать нормальное распределение с центром в 0 с широким стандартным отклонением. Мы предполагаем, что измерения обычно распределены с шумом, который мы смоделируем при помощи равномерного распределения.

Вместо того, чтобы создавать отдельные модели, иерархическая модель создает параметры на уровне округов, благодаря чему округа рассматриваются не по отдельности. $$\alpha_{c} \sim \mathcal{N}(\mu_{\alpha}, \sigma_{\alpha}^2)$$ $$\beta_{c} \sim \mathcal{N}(\mu_{\beta}, \sigma_{\beta}^2)$$ Эти распределения впоследствии используются, чтобы влиять на распределение 𝛼 и 𝛽 каждого округа. Таким образом, мы предполагаем, что свободный коэффициент 𝛼 и угловой коэффициент 𝛽 происходят из нормального распределения, сосредоточенного вокруг их соответствующего среднего по группе $\mu$ с некоторым стандартным отклонением $\sigma^2$, значения (или, скорее, исходные значения), которое мы также оцениваем. Вот почему это называется многоуровневым моделированием.

Величина $mu_a$ показывает средние уровни радона в группах. $mu_b$ говорит, что отсутствие фундамента значительно снижает уровни радона (график плотности лежит ниже нуля). Мы также можем увидеть, посмотрев на график $a$, что между округами довольно много различий в уровнях радона (каждый цвет соответствует одному округу); Различная ширина связана с тем, насколько мы уверены в каждой оценке параметров - чем больше измерений в округе, тем выше будет наша достоверность.

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

Глядя на округ «CASS», мы видим, что неиерархическая оценка сильно смещена: поскольку данные этого округа содержат только домохозяйства с подвалом, регрессионная модель даёт бессмысленный результат, показывая, что даже в доме без подвала следует ожидать уровни радона как в радоновой шахте.

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

Самостоятельное задание

Изучая классические, не байесовские многоуровневые модели, мы пытались найти различия в ценностях людей из разных стран в базе данные European Social Survey. Попробуйте сделать тоже самое с опорой на байесовский подход и проведите необходимые диагностики. Описание всех переменных из ESS доступно в протоколе исследования, сами данные можно скачать по ссылке.

Дополнительные ссылки

• Введение в Байесовские методы
• Николенко. Введение в байесовский вывод
• Дьяконов. Байесовский подход
• Байесовские методы машинного обучения (курс лекций) / 2018
• GLM: Linear regression
• GLM: Hierarchical Linear Regression
• О формуле Байеса, прогнозах и доверительных интервалах
• MCMC-сэмплинг для тех, кто учился, но ничего не понял
• PyMC3 — MCMC и не только