7 Коэффициент Байеса
7.1 Формула Байеса опять
\[P(\theta|Data) = \frac{P(Data|\theta) \times P(\theta) }{P(Data)}\]
Рассмотрим какой-то простой случай, который мы уже видели много раз.
Немного упрощая данные из статьи (Rosenbach 2003: 394), можно сказать что носители британского английского предпочитают s-генитив (90%) of-генитиву (10%). Проведите байесовский апдейт, если Вы наблюдаете в интервью британского актера из 120 контекстов 92 s-генитивов. Априорное распределение берите соразмерное данным.
Если мы не будем следовать простой дорожкой, которую мы обсуждали несколько разделов назад, а будем все делать согласно формуле Байеса, то получатся следующие компоненты:
- априорное распределение
tibble(x = seq(0, 1, 0.001),
prior = dbeta(x = x, shape1 = 120*0.9, shape2 = 120*0.1)) %>%
ggplot(aes(x, prior))+
geom_line(color = "red")
- функция правдоподобия
tibble(x = seq(0, 1, 0.001),
likelihood = dbinom(x = 92, size = 120, prob = x)) %>%
ggplot(aes(x, likelihood))+
geom_line()
- их произведение (пропорционально апостериорному распределению)
tibble(x = seq(0, 1, 0.001),
prior = dbeta(x = x, shape1 = 120*0.9, shape2 = 120*0.1),
likelihood = dbinom(x = 92, size = 120, prob = x),
product = prior*likelihood) %>%
ggplot(aes(x, product))+
geom_line()
- предельное правдоподобие, которое позволяет сделать получившееся распределение распределением вероятностей
<- integrate(function(p){
marginal_likelihood dbinom(92, 120, p) * dbeta(p, 120*0.9, 120*0.1)},
lower = 0,
upper = 1)
marginal_likelihood
0.0009531395 with absolute error < 0.000044
- … и в результате получается апостериорное распределение!
tibble(x = seq(0, 1, 0.001),
prior = dbeta(x = x, shape1 = 120*0.9, shape2 = 120*0.1),
likelihood = dbinom(x = 92, size = 120, prob = x),
product = prior*likelihood,
posterior = product/marginal_likelihood[[1]]) %>%
ggplot(aes(x, posterior))+
geom_line(color = "darkgreen")+
geom_line(aes(y = prior), color = "red")
… которое мы умеем доставать и быстрее:
tibble(x = seq(0, 1, 0.001),
prior = dbeta(x = x, shape1 = 120*0.9, shape2 = 120*0.1),
likelihood = dbinom(x = 92, size = 120, prob = x),
product = prior*likelihood,
posterior = product/marginal_likelihood[[1]],
posterior_2 = dbeta(x = x, shape1 = 120*0.9+92, shape2 = 120*0.1+120-92)) %>%
ggplot(aes(x, posterior))+
geom_line(color = "darkgreen", size = 2)+
geom_line(aes(y = prior), color = "red")+
geom_line(aes(y = posterior_2), linetype = 2, color = "yellow")
Представим себе, что у нас есть \(k\) гипотез \(M\). Тогда формула Байеса может выглядеть вот так:
\[P(M_k|Data) = \frac{P(Data|M_k) \times P(M_k) }{P(Data)}\] В данном занятии мы рассмотрим только случай двух модели, но можно рассматривать и случаи, когда моделей много. Посмотрим на соотношение апостериорных распределений двух моделей:
\[\underbrace{\frac{P(M_1 \mid Data)}{P(M_2 \mid Data)}}_{\text{posterior odds}} = \frac{\frac{P(Data|M_1) \times P(M_1) }{P(Data)}}{\frac{P(Data|M_2) \times P(M_2) }{P(Data)}}=\underbrace{\frac{P(Data \mid M_1)}{P(Data \mid M_2)}}_{\text{Bayes factor}}\times\underbrace{\frac{P(M_1)}{P(M_2)}}_{\text{prior odds}}\] Таким образом байесовский коэффициент это соотношение апосториорных распределений деленное на соотношение априорных распределений.
\[BF_{12}= \frac{P(M_1 \mid Data)/P(M_2 \mid Data)}{P(M_1)/P(M_2)}=\frac{P(M_1 \mid Data)\times P(M_2)}{P(M_2 \mid Data)\times P(M_1)}\]
В результате получается, что коэффициент Байеса – это соотношение предельных правдоподобий (знаменатель теоремы Байеса):
\[BF_{12}= \frac{P(Data|\theta, M_1))}{P(Data|\theta, M_2))}=\frac{\int P(Data|\theta, M_1)\times P(\theta|M_1)}{\int P(Data|\theta, M_2)\times P(\theta|M_2)}\]
Важно заметить, что если вероятности априорных моделей равны, то байесовский коэффициент равен просто соотношению функций правдоподобия.
Надо отметить, что не все тепло относятся к сравнению моделей байесовским коэффициентом (см. Gelman, Rubin 1994).
7.2 Категориальные данные
Для примера обратимся снова к датасету, который содержит спамерские и обычные смс-сообщения, выложенному UCI Machine Learning на kaggle, и при помощи пакета udpipe
токенизируем и определим часть речи:
<- read_csv("https://raw.githubusercontent.com/agricolamz/2022_da4l/master/data/spam_sms_pos.csv")
sms_pos glimpse(sms_pos)
Rows: 34
Columns: 3
$ type <chr> "ham", "ham", "ham", "ham", "ham", "ham", "ham", "ham", "ham", "h…
$ upos <chr> "ADJ", "ADP", "ADV", "AUX", "CCONJ", "DET", "INTJ", "NOUN", "NUM"…
$ n <dbl> 4329, 5004, 5832, 5707, 1607, 3493, 1676, 12842, 1293, 2424, 1144…
%>%
sms_pos group_by(type) %>%
mutate(ratio = n/sum(n),
upos = fct_reorder(upos, n, mean, .desc = TRUE)) %>%
ggplot(aes(type, ratio))+
geom_col()+
geom_label(aes(label = round(ratio, 3)), position = position_stack(vjust = 0.5))+
facet_wrap(~upos, scales = "free_y")
Давайте полученные доли считать нашей моделью: сумма всех чисел внутри каждого типа (ham
/spam
) дает в сумме 1. Мы получили новое сообщение:
Call FREEPHONE 0800 542 0825 now!
Модель udpipe
разобрала его следующим образом:
VERB NUM NUM NUM NUM ADV PUNCT
Если мы считаем наши модели равновероятными:
<- tibble(model = c("ham", "spam"),
first_update prior = 0.5,
likelihood = c(0.135, 0.096),
product = prior*likelihood,
marginal_likelihood = sum(product),
posterior = product/marginal_likelihood)
first_update
Если же мы примем во внимание, что наши классы не равноправны, то сможем посчитать это нашим априорным распределением для моделей.
%>%
sms_pos uncount(n) %>%
count(type) %>%
mutate(ratio = n/sum(n)) ->
class_ratio class_ratio
<- tibble(model = c("ham", "spam"),
second_update prior = class_ratio$ratio,
likelihood = c(0.135, 0.096),
product = prior*likelihood,
marginal_likelihood = sum(product),
posterior = product/marginal_likelihood)
second_update
# Bayes factor
$marginal_likelihood[1]/first_update$marginal_likelihood[1] second_update
[1] 1.098469
7.4 Биномиальные данные
Рассмотрим простенькую задачу, которую мы видели раньше:
Немного упрощая данные из статьи (Rosenbach 2003: 394), можно сказать что носители британского английского предпочитают s-генитив (90%) of-генитиву (10%), а носители американского английского предпочитают s-генитив (85%) of-генитиву (15%). Мы наблюдаем актера, который в интервью из 120 контекстов использует в 92 случаях s-генитивы. Сравните модели при помощи байесовского коэффициента.
tibble(x = seq(0, 1, by = 0.001),
y = dbeta(x, 120*0.9, 120*0.1),
z = dbeta(x, 120*0.85, 120*0.15)) %>%
ggplot(aes(x, y))+
geom_line(color = "red")+
geom_line(aes(y = z), color = "lightblue")+
geom_vline(xintercept = 92/120, linetype = 2)
<- function(p) dbinom(92, 120, p) * dbeta(p, 120*0.9, 120*0.1)
m1 <- function(p) dbinom(92, 120, p) * dbeta(p, 120*0.85, 120*0.15)
m2
integrate(m1, 0, 1)[[1]]/integrate(m2, 0, 1)[[1]]
[1] 0.0672068
В работе (Coretta 2016) собраны данные
длительности исландских гласных (столбец vowel.dur
).
Отфильтруйте данные, произнесенные носителем tt01
(переменная speaker
), посчитайте байесовский коэффициент
(\(B_{12}\)) для двух априорных
моделей:
- нормального распределения со средним 87 и стандартным отклонением 25. (\(m_1\))
- нормального распределения со средним 85 и стандартным отклонением 30. (\(m_2\))
Ответ округлите до трёх или менее знаков после запятой.