Cette question est un suivi technique de cette question .
J'ai du mal à comprendre et à reproduire le modèle présenté dans Raftery (1988): Inférence pour le paramètre binomial : une approche bayésienne hiérarchique dans WinBUGS / OpenBUGS / JAGS. Il ne s'agit pas seulement de code, donc il devrait être sur le sujet ici.
Contexte
Soit un ensemble de décomptes de succès d'une distribution binomiale avec et inconnus . De plus, je suppose que suit une distribution de Poisson avec le paramètre (comme discuté dans l'article). Ensuite, chaque a une distribution de Poisson avec une moyenne . Je veux spécifier les a priori en termes de et .
En supposant que je n'ai pas de bonnes connaissances préalables sur ou , je veux attribuer des a priori non informatifs à la fois à et à . Disons que mes prieurs sont et .
L'auteur utilise un a priori impropre de mais WinBUGS n'accepte pas les a priori incorrects.
Exemple
Dans le document (page 226), les décomptes de succès suivants des waterbucks observés sont fournis: . Je veux estimer , la taille de la population.
Voici comment j'ai essayé de travailler sur l'exemple dans WinBUGS ( mis à jour après le commentaire de @ 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)
Le modèle ne SILL pas convergé bien après 500'000 échantillons avec 20'000 burn-in échantillons. Voici la sortie d'une exécution 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
Des questions
De toute évidence, il me manque quelque chose, mais je ne vois pas quoi exactement. Je pense que ma formulation du modèle est erronée quelque part. Mes questions sont donc:
- Pourquoi mon modèle et sa mise en œuvre ne fonctionnent-ils pas?
- Comment le modèle donné par Raftery (1988) pourrait-il être formulé et mis en œuvre correctement?
Merci de votre aide.
la source
mu=lambda/theta
et remplacern ~ dpois(lambda)
parn ~ dpois(mu)
Réponses:
Eh bien, puisque vous avez fait fonctionner votre code, il semble que cette réponse soit un peu trop tardive. Mais j'ai déjà écrit le code, alors ...
Pour ce que ça vaut, c'est le même modèle * adapté(N,θ)
rstan
. Il est estimé en 11 secondes sur mon ordinateur portable grand public, atteignant une taille d'échantillon efficace plus élevée pour nos paramètres d'intérêt en moins d'itérations.Notez que je lance
theta
en 2 simplex. C'est juste pour la stabilité numérique. La quantité d'intérêt esttheta[1]
;theta[2]
est évidemment une information superflue.* Comme vous pouvez le voir, le résumé postérieur est pratiquement identique, et la promotion de en une quantité réelle ne semble pas avoir un impact substantiel sur nos inférences.N
Le quantile de 97,5% pour est 50% plus grand pour mon modèle, mais je pense que c'est parce que l'échantillonneur de Stan est meilleur pour explorer la gamme complète de la partie postérieure qu'une simple marche aléatoire, de sorte qu'il peut plus facilement en faire la queue. Je peux me tromper, cependant.N
En prenant les valeurs de générées à partir de stan, je les utilise pour tracer des valeurs prédictives postérieures ˜ y . Il ne faut pas s'étonner que la moyenne des prédictions postérieures ˜ y soit très proche de la moyenne des données de l'échantillon!N,θ y~ y~
rstan
Le code ci-dessous peut confirmer que nos résultats de stan ont du sens.
rstan
la source
n
ne peut pas être déclaré comme un entier et je ne connaissais pas de solution de contournement pour le problème.rstan
résultats sont les plus corrects. [1] stats.stackexchange.com/questions/114366/…Merci encore à @ StéphaneLaurent et @ user777 pour leur précieuse contribution dans les commentaires. Après quelques ajustements du prieur pourλ Je peux maintenant reproduire les résultats de l'article de Raftery (1988).
Voici mon script d'analyse et mes résultats avec JAGS et R:
Le calcul a pris environ 98 secondes sur mon ordinateur de bureau.
Les résultats sont:
La moyenne postérieure deN est 522 et la médiane postérieure est 233 . J'ai calculé le mode postérieur deN sur l'échelle logarithmique et retransformé l'estimation:
Et le 80% -HPD deN est:
Le mode postérieur pourN est 150 et le 80% -HPD est ( 78 ; 578 ) qui est très proche des limites données dans le document qui sont ( 80 ; 598 ) .
la source