Код для базового иерархического байесовского анализа - PullRequest
Новая
       103

Код для базового иерархического байесовского анализа

0 голосов
/ 03 марта 2019

У меня есть этот код, написанный в winBUGS: 

n <- 100

x1 <- rbinom(n,1,.7)

x2 <- rbinom(n,1,.5)

sum(x1)

sum(x2)

model{

x1 ~ dbin(p1, n) x2 ~ dbin(p2, n) p1 ~ dbeta(a1, b1) p2 ~ dbeta(a2,b2)

diff <- p1 - p2 p.value <- step(diff)

} list(n = 100, x1 = 70, x2 = 54, a1 = 1, b1 = 1, a2 = 1, b2 = 1)

У меня возникли проблемы с выполнением этого в R / JAGS.На самом деле, я даже не совсем уверен, что этот код пытается сделать (я думаю, чтобы вычислить постеры?).Я никогда раньше не использовал winBUGS, и я новичок в R. Это также мой первый класс Байеса, и я совершенно растерялся, когда был введен код.

Кроме того, как бы я вычислил среднее значение истандартное отклонение для разности пропорций?или как найти апостериорную вероятность того, что p1 больше, чем p2, и имеет ли это значение?

1 Ответ

0 голосов
/ 06 марта 2019

Если вам нужна помощь в настройке rjags , я ожидаю, что ваш TA или одноклассники могут помочь вам в этом.Этот ответ сфокусирован на объяснении того, что делает код.

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

Генеративное / прямое моделирование

Первое - это прямое моделирование, где мы уже знаем, какова вероятность перевернуть голову для двух монет. 

## how many times we will flip each coin
n <- 100

## coin 1 has 70% chance to land on heads
## simulate n flips
x1 <- rbinom(n, 1, .7)

## coin 2 has 50% chance to land on heads
## simulate n flips
x2 <- rbinom(n, 1, .5)

## count number of coin 1 heads
sum(x1) # 70

## count number of coin 2 heads
sum(x2) # 54

Задние вероятности

Теперь мы возьмем эти смоделированные данные и попробуем повернуть эксперимент вспять.То есть, если мы наблюдали 70 голов одной монеты и 54 головы другой, можем ли мы что-то сказать о вероятностях, что каждая монета должна перевернуть голову?В частности, код будет спрашивать: « При едином предположении об истинной вероятности, что каждая монета имеет, какова (задняя) вероятность того, что монета 1 с большей вероятностью перевернет голову, чем монета 2? "

Вычисление вероятностей

Способ, которым JAGS будет делать это, - MCMC, где выборки берутся из пространства всех возможных конфигураций параметров (в данном случае значения p1 и p2), взвешенный по тому, насколько вероятно, что каждая конфигурация генерирует наблюдаемые данные.Следовательно, главное, что JAGS-код должен выполнить, это определить, как генерировать значения параметров и как вычислять вероятность данных, учитывая эти значения.Это первая часть кода модели.

Вычисление интересующих значений

Вторая часть кода в разделе модели тестирует поставленный нами вопрос.Это делается путем вычисления некоторых дополнительных переменных, которые не влияют на вероятность конфигураций.Вместо этого они сообщают нам информацию о конфигурации выборки.В частности, diff будет отслеживать распределение различий между двумя вероятностями монет, а p.value будет отслеживать, как часто p1 больше, чем p2.

model{

## COMPUTE LIKELIHOODS
## compute likelihood that heads resulted from coins
## with given probabilities after `n` coin flips
x1 ~ dbin(p1, n) 
x2 ~ dbin(p2, n) 

## SAMPLE PARAMETERS
## randomly select coin probabilities from Beta distributions
## Note that since these are all 1, it's really just a Uniform[0,1]
p1 ~ dbeta(a1, b1)
p2 ~ dbeta(a2, b2)

## HYPOTHESIS TEST
## compute how often coin1's probability is greater than coin2's
diff <- p1 - p2 
p.value <- step(diff)
} 

## INPUT VALUES
list(n = 100, x1 = 70, x2 = 54, a1 = 1, b1 = 1, a2 = 1, b2 = 1)
...