Этот вопрос является техническим продолжением этого вопроса .
У меня проблемы с пониманием и тиражированием модели, представленной в Raftery (1988): Вывод для биномиального параметра : иерархический байесовский подход в WinBUGS / OpenBUGS / JAGS. Речь идет не только о коде, поэтому он должен быть здесь по теме.
Фон
Пусть - набор отсчетов успеха из биномиального распределения с неизвестными N и θ . Далее, я предполагаю, что N следует распределению Пуассона с параметром μ (как обсуждалось в статье). Тогда каждый x i имеет распределение Пуассона со средним λ = µ θ . Я хочу указать приоры в терминах λ и θ .
Предполагая, что у меня нет хороших предварительных знаний о или θ , я хочу назначить неинформативные априорные значения как λ, так и θ . Скажем, мои априорные являются λ ~ G м м ( 0,001 , 0,001 ) и & thetas ; ~ U н я е ö г м ( 0 , 1 ) .
Автор использует неправильный априор но WinBUGS не принимает неправильные априорные значения.
пример
В статье (стр. 226) приводятся следующие показатели успеха наблюдаемых водяных козлов: . Я хочу оценить N , численность населения.
Вот как я пытался отработать пример в WinBUGS ( обновлено после комментария @ Stéphane Laurent):
model {
# Likelihood
for (i in 1:N) {
x[i] ~ dbin(theta, n)
}
# Priors
n ~ dpois(mu)
lambda ~ dgamma(0.001, 0.001)
theta ~ dunif(0, 1)
mu <- lambda/theta
}
# Data
list(x = c(53, 57, 66, 67, 72), N = 5)
# Initial values
list(n = 100, lambda = 100, theta = 0.5)
list(n = 1000, lambda = 1000, theta = 0.8)
list(n = 5000, lambda = 10, theta = 0.2)
Модель не подоконник не сходится хорошо после 500'000 образцов с 20'000 ожогом в образцах. Вот результат запуска JAGS:
Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 5
n.sims = 480000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
lambda 63.081 5.222 53.135 59.609 62.938 66.385 73.856 1.001 480000
mu 542.917 1040.975 91.322 147.231 231.805 462.539 3484.324 1.018 300
n 542.906 1040.762 95.000 147.000 231.000 462.000 3484.000 1.018 300
theta 0.292 0.185 0.018 0.136 0.272 0.428 0.668 1.018 300
deviance 34.907 1.554 33.633 33.859 34.354 35.376 39.213 1.001 43000
Вопросов
Я явно что-то упускаю, но не вижу, что именно. Я думаю, что моя формулировка модели где-то не так. Итак, мои вопросы:
- Почему моя модель и ее реализация не работают?
- Как правильно сформулировать и реализовать модель, предложенную Рэфтери (1988)?
Спасибо за вашу помощь.
источник
mu=lambda/theta
и заменитьn ~ dpois(lambda)
наn ~ dpois(mu)
Ответы:
Ну, так как ваш код работает, похоже, что этот ответ уже слишком поздний. Но я уже написал код, так что ...
Для чего это стоит, это та же * модель подходит с(N,θ)
rstan
. Он оценивается в 11 секунд на моем потребительском ноутбуке, что позволяет получить более эффективный размер выборки для наших параметров за меньшее количество итераций.Обратите внимание, что я разыгрываю
theta
как 2-симплекс. Это только для численной стабильности. Количество процентов составляетtheta[1]
; очевидноtheta[2]
это лишняя информация.* Как видите, апостериорная сводка практически идентична, и повышение до реального количества, по-видимому, не оказывает существенного влияния на наши выводы.N
Взяв значения сгенерированные из stan, я использую их, чтобы нарисовать апостериорные прогностические значенияN,θ y~ y~
rstan
Код ниже может подтвердить, что наши результаты из stan имеют смысл.
Гектометр Это не совсем то, что я ожидал. Оценка сетки для квантиля 97,5% ближе к результатам JAGS, чем к(0,1)×{N|N∈Z∧N≥72)}
rstan
результатам. В то же время я не считаю, что результаты сетки следует воспринимать как Евангелие, потому что оценка сетки делает несколько довольно грубых упрощений: разрешение сетки не слишком хорошее, с одной стороны, а с другой - ( ) утверждая, что полная вероятность в сетке должна быть 1, поскольку мы должны нарисовать границу (и конечные точки сетки), чтобы задача была вычислимой (я все еще жду на бесконечной ОЗУ). По правде говоря, наша модель имеет положительную вероятностьисточник
n
не может быть объявлено как целое число, и я не знаю обходной путь для этой проблемы.rstan
результаты являются более правильными. [1] stats.stackexchange.com/questions/114366/…Вот мой сценарий анализа и результаты с использованием JAGS и R:
На моем настольном компьютере вычисления заняли около 98 секунд.
Результаты:
И 80% -HPDN
источник