Estimation bayésienne de

16

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 N : 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 x=(x1,,xn) un ensemble de décomptes de succès d'une distribution binomiale avec N et inconnus θ. De plus, je suppose que N suit une distribution de Poisson avec le paramètre μ (comme discuté dans l'article). Ensuite, chaque xi 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 N ou θ , je veux attribuer des a priori non informatifs à la fois à λ et à θ . Disons que mes prieurs sont λGamma(0.001,0.001) et θUniform(0,1) .

L'auteur utilise un a priori impropre de p(N,θ)N1 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: 53,57,66,67,72 . Je veux estimer N , 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.

COOLSerdash
la source
2
Après le papier, vous devez ajouter mu=lambda/thetaet remplacer n ~ dpois(lambda)parn ~ dpois(mu)
Stéphane Laurent
@ StéphaneLaurent Merci pour la suggestion. J'ai changé le code en conséquence. Malheureusement, le modèle ne converge toujours pas.
COOLSerdash
1
Que se passe-t-il lorsque vous échantillonnez ? N<72
Sycorax dit Réintégrer Monica
1
Si , la probabilité est nulle, car votre modèle suppose qu'il y a au moins 72 waterbuck. Je me demande si cela cause des problèmes pour l'échantillonneur. N<72
Sycorax dit Réintégrer Monica
3
Je ne pense pas que le problème soit la convergence. Je pense que le problème est que l'échantillonneur est peu performant en raison du degré élevé de corrélation au niveau des multiples niveaux du est faible, tandis que n e f f est faible par rapport au nombre total d'itérations. Je suggère que le calcul du postérieur directement, par exemple, sur une grille θ , N . R^neffθ,N
Sycorax dit Réintégrer Monica

Réponses:

7

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é 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.(N,θ)

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 thetaen 2 simplex. C'est juste pour la stabilité numérique. La quantité d'intérêt est theta[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

            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

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~

N.samples   <- round(extract(fit, "N")[[1]])
theta.samples   <- extract(fit, "theta")[[1]]
y_pred  <- rbinom(50000, size=N.samples, prob=theta.samples[,1])
mean(y_pred)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  32.00   58.00   63.00   63.04   68.00  102.00 

rstanNy¯=θN

posterior over a grid

Le code ci-dessous peut confirmer que nos résultats de stan ont du sens.

theta   <- seq(0+1e-10,1-1e-10, len=1e2)
N       <- round(seq(72, 5e5, len=1e5)); 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)    
approx(y=N, x=cumsum(rowSums(post))/sum(rowSums(post)), xout=0.975)
$x
[1] 0.975

$y
[1] 3236.665

rstan(0,1)×{N|NZN72)}

Sycorax dit de réintégrer Monica
la source
+1 et accepté. Je suis impressionné! J'ai également essayé d'utiliser Stan pour une comparaison, mais je n'ai pas pu transférer le modèle. Mon modèle prend environ 2 minutes pour estimer.
COOLSerdash
Le seul problème avec stan pour ce problème est que tous les paramètres doivent être réels, ce qui le rend un peu gênant. Mais comme vous pouvez pénaliser la probabilité de logarithme par n'importe quelle fonction arbitraire, il vous suffit de passer par la difficulté de la programmer ... Et de creuser les fonctions composées pour le faire ...
Sycorax dit Reinstate Monica
Oui! C'était exactement mon problème. nne peut pas être déclaré comme un entier et je ne connaissais pas de solution de contournement pour le problème.
COOLSerdash
Environ 2 minutes sur mon bureau.
COOLSerdash
1
@COOLSerdash Vous pouvez être intéressé par [cette] [1] question, où je demande lequel des résultats de la grille ou les rstanrésultats sont les plus corrects. [1] stats.stackexchange.com/questions/114366/…
Sycorax dit
3

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:

#===============================================================================================================
# Load packages
#===============================================================================================================

sapply(c("ggplot2"
         , "rjags"
         , "R2jags"
         , "hdrcde"
         , "runjags"
         , "mcmcplots"
         , "KernSmooth"), library, character.only = TRUE)

#===============================================================================================================
# Model file
#===============================================================================================================

cat("
    model {

    # Likelihood    
    for (i in 1:N) {
      x[i] ~ dbin(theta, n)
    }

    # Prior       
    n ~ dpois(mu)
    lambda ~ dgamma(0.005, 0.005)
#     lambda ~ dunif(0, 1000)
    mu <- lambda/theta
    theta ~ dunif(0, 1)    
}    
", file="jags_model_binomial.txt")


#===============================================================================================================
# Data
#===============================================================================================================

data.list <- list(x = c(53, 57, 66, 67, 72, NA), N = 6) # Waterbuck example from Raftery (1988)

#===============================================================================================================
# Inits
#===============================================================================================================

jags.inits <- function() { 
  list(
    n = sample(max(data.list$x, na.rm = TRUE):1000, size = 1) 
    , theta = runif(1, 0, 1)
    , lambda = runif(1, 1, 10)
#     , cauchy  = runif(1, 1, 1000)
    #     , mu = runif(1, 0, 5)
  )
}

#===============================================================================================================
# Run the chains
#===============================================================================================================

# Parameters to store

params <- c("n"
            , "theta"
            , "lambda"
            , "mu"
            , paste("x[", which(is.na(data.list[["x"]])), "]", sep = "")
)

# MCMC settings

niter <- 500000 # number of iterations
nburn <- 20000  # number of iterations to discard (the burn-in-period)
nchains <- 5    # number of chains

# Run JAGS

out <- jags(
  data                 = data.list
  , parameters.to.save = params
  , model.file         = "jags_model_binomial.txt"
  , n.chains           = nchains
  , n.iter             = niter
  , n.burnin           = nburn
  , n.thin             = 50
  , inits              = jags.inits
  , progress.bar       = "text")

Le calcul a pris environ 98 secondes sur mon ordinateur de bureau.

#===============================================================================================================
# Inspect results
#===============================================================================================================

print(out
      , digits = 2
      , intervals = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9,  0.975))

Les résultats sont:

Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
 5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 50
 n.sims = 48000 iterations saved
         mu.vect sd.vect  2.5%    10%    25%    50%    75%     90%   97.5% Rhat n.eff
lambda     62.90    5.18 53.09  56.47  59.45  62.74  66.19   69.49   73.49    1 48000
mu        521.28  968.41 92.31 113.02 148.00 232.87 467.10 1058.17 3014.82    1  1600
n         521.73  968.54 95.00 114.00 148.00 233.00 467.00 1060.10 3028.00    1  1600
theta       0.29    0.18  0.02   0.06   0.13   0.27   0.42    0.55    0.66    1  1600
x[6]       63.03    7.33 49.00  54.00  58.00  63.00  68.00   72.00   78.00    1 36000
deviance   34.88    1.53 33.63  33.70  33.85  34.34  35.34   36.81   39.07    1 48000

La moyenne postérieure de N 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:

jagsfit.mcmc <- as.mcmc(out)
jagsfit.mcmc <- combine.mcmc(jagsfit.mcmc)

hpd.80 <- hdr.den(log(as.vector(jagsfit.mcmc[, "n"])), prob = c(80), den = bkde(log(as.vector(jagsfit.mcmc[, "n"])), gridsize = 10000))

exp(hpd.80$mode)

[1] 149.8161

Et le 80% -HPD de N est:

(hpd.ints <- HPDinterval(jagsfit.mcmc, prob = c(0.8)))

               lower      upper
deviance 33.61011007  35.677810
lambda   56.08842502  69.089507
mu       72.42307587 580.027182
n        78.00000000 578.000000
theta     0.01026193   0.465714
x[6]     53.00000000  71.000000

Le mode postérieur pour N 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).

COOLSerdash
la source