J'ai un modèle pour obtenir des estimations bayésiennes de la taille de la population et probabilité de détection dans une distribution binomiale uniquement basée sur le nombre observé d’objets observés : pour . Pour simplifier, nous supposons que N est fixé à la même valeur inconnue pour chaque y_i . Dans cet exemple, y = 53,57,66,67,73 .
Ce modèle, lorsqu'il est estimé en rstan
, diverge des résultats obtenus à partir d'une approximation en grille de la partie postérieure. J'essaie de comprendre pourquoi. (Les lecteurs intéressés pourraient trouver que cette question fait suite à ma réponse ici .)
rstan
Approximation
Pour référence, il s'agit du code rstan.
raftery.model <- "
data{
int I;
int y[I];
}
parameters{
real<lower=max(y)> N;
simplex[2] theta;
}
transformed parameters{
}
model{
vector[I] Pr_y;
for(i in 1:I){
Pr_y[i] <- binomial_coefficient_log(N, y[i])
+multiply_log(y[i], theta[1])
+multiply_log((N-y[i]), theta[2]);
}
increment_log_prob(sum(Pr_y));
increment_log_prob(-log(N));
}
"
raft.data <- list(y=c(53,57,66,67,72), I=5)
system.time(fit.test <- stan(model_code=raftery.model, data=raft.data,iter=10))
system.time(fit <- stan(fit=fit.test, data=raft.data,iter=10000,chains=5))
Notez que je lance theta
en 2 simplex. C'est juste pour la simplicité. La quantité d'intérêt est theta[1]
; theta[2]
est évidemment une information superflue.
De plus, est une valeur réelle ( rstan
accepte uniquement les paramètres à valeur réelle car c'est une méthode de gradient), j'ai donc écrit une distribution binomiale à valeur réelle.
Résultats Rstan
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
N 1078.75 256.72 15159.79 94.44 148.28 230.61 461.63 4575.49 3487 1
theta[1] 0.29 0.00 0.19 0.01 0.14 0.27 0.42 0.67 2519 1
theta[2] 0.71 0.00 0.19 0.33 0.58 0.73 0.86 0.99 2519 1
lp__ -19.88 0.02 1.11 -22.89 -20.31 -19.54 -19.09 -18.82 3339 1
Approximation de la grille
L'approximation de la grille a été produite comme ci-dessous. Les contraintes de mémoire m'empêchent de faire une grille plus fine sur mon ordinateur portable.
theta <- seq(0+1e-10,1-1e-10, len=1e3)
N <- round(seq(72, 5000, len=1e3)); N[2]-N[1]
grid <- expand.grid(N,theta)
y <- c(53,57,66,67,72)
raftery.prob <- function(x, z=y){
N <- x[1]
theta <- x[2]
exp(sum(dbinom(z, size=N, prob=theta, log=T)))/N
}
post <- matrix(apply(grid, 1, raftery.prob), nrow=length(N), ncol=length(theta),byrow=F)
post.norm <- post/sum(post)
J'ai utilisé l'approximation de la grille pour produire cet affichage de la densité postérieure. On voit que le postérieur est en forme de banane; ce type de postérieur peut être problématique pour la métrique euclidienne HMC. (La gravité de la forme de la banane est en fait supprimée ici car est sur l'échelle logarithmique.) Si vous pensez à la forme de la banane pendant une minute, vous vous rendrez compte qu'elle doit se trouver sur la ligne . (De plus, l'approximation de la grille affichée dans ce graphique n'est pas normalisée pour des raisons de clarté - sinon la banane est un peu trop étroite pour être clairement identifiée.)
Résultats d'approximation de la grille
do.call(cbind, lapply(c(0.025, .25, .5, .75, .975), function(quantile){
approx(y=N, x=cumsum(rowSums(post.norm))/sum(post.norm), xout=quantile)
}))
[,1] [,2] [,3] [,4] [,5]
x 0.025 0.25 0.5 0.75 0.975
y 92.55068 144.7091 226.7845 443.6359 2475.398
Discussion
Le quantile de 97,5% pour est beaucoup plus grand dans mon modèle que pour l'approximation de la grille, mais ses quantiles sont similaires à l'approximation de la grille autrement. J'interprète cela comme indiquant que les deux méthodes sont généralement d'accord. Cependant, je ne sais pas comment interpréter l'écart dans le quantile de 97,5%.rstan
J'ai développé plusieurs explications possibles pour ce qui pourrait expliquer la divergence entre l'approximation de la grille et les résultats de l' rstan
échantillonnage HMC-NUTS, mais je ne sais pas comment comprendre si une, les deux ou aucune explication n'est correcte.
- Rstan a tort et la grille est correcte. La densité en forme de banane est problématique
rstan
, d'autant que dérive vers , donc ces quantités de queue ne sont pas fiables. Nous pouvons voir de la parcelle de la partie postérieure sur la grille que la queue est très forte à des valeurs plus grandes . - Rstan a raison et la grille est fausse. La grille fait deux approximations qui peuvent saper les résultats. Premièrement, la grille n'est qu'un ensemble fini de points sur un sous-espace postérieur, c'est donc une approximation approximative. En second lieu , parce que c'est un sous - espace fini, nous avoir faussement déclaré qu'il y ait 0 probabilité a posteriori sur les valeurs plus grand que notre plus grande valeur de grille pour . De même, est meilleur pour entrer dans les queues de la grille, donc ses quanitles de queue sont corrects.
rstan
J'avais besoin de plus d'espace pour clarifier un point de la réponse de Juho. Si je comprends bien, nous pouvons intégrer hors de la partie postérieure pour obtenir la distribution bêta-binomiale:
Dans notre cas, et car nous avons un préalable uniforme sur . Je crois que le postérieur devrait alors être où car . Mais cela semble diverger énormément de la réponse de Juho. Où ai-je mal tourné?
la source