Estimación bayesiana de $N$ de una distribución binomial

Question

Esta pregunta es una técnica de seguimiento de esta cuestión.

Tengo problemas para entender y replicar el modelo presentado en Raftery (1988): la Inferencia para el binomio $N$ parámetro: un enfoque Bayesiano jerárquico en WinBUGS/OpenBUGS/ENTRECORTADO. No es sólo acerca de código, aunque por lo que debe ser un tema aquí.

De fondo

Deje $x=(x_{1},\ldots,x_{n})$ ser una serie de éxito cuenta de una distribución binomial con desconocidos $N$$\theta$. Además, supongo que $N$ sigue una distribución de Poisson con parámetro de $\mu$ (como se discutió en el papel). A continuación, cada una de las $x_{i}$ tiene una distribución de Poisson con una media de $\lambda = \mu \theta$. Quiero especificar las prioridades en términos de$\lambda$$\theta$.

Suponiendo que no tengo buen conocimiento previo acerca de la $N$ o $\theta$, quiero asignar los informativos de los priores tanto $\lambda$$\theta$. Decir, mi priores se $\lambda\sim \mathrm{Gamma}(0.001, 0.001)$$\theta\sim \mathrm{Uniform}(0, 1)$.

El autor utiliza un inadecuado antes de $p(N,\theta)\propto N^{-1}$ pero WinBUGS no acepta inadecuado de los priores.

Ejemplo

En el papel (página 226), el siguiente éxito de los recuentos observados antílopes acuáticos: $53, 57, 66, 67, 72$. Quiero estimación $N$, el tamaño de la población.

Aquí es cómo he intentado averiguar el ejemplo en WinBUGS (actualizado después de @Stéphane Laurent comentario):

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)

El modelo de umbral no convergen muy bien después de 500'000 muestras con 20'000 quemar-en las muestras. Aquí está la salida de un ENTRECORTADO ejecutar:

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

Preguntas

Está claro que me estoy perdiendo algo, pero no puedo ver exactamente qué. Creo que mi formulación del modelo que está mal en alguna parte. Así que mis preguntas son:

• ¿Por qué mi modelo y su aplicación no funciona?
• ¿Cómo podría el modelo dado por Raftery (1988) ser formulado y aplicado correctamente?

Gracias por tu ayuda.

Answer 1

Bien, ya tienes tu código de trabajo, parece que esta respuesta es un poco demasiado tarde. Pero ya he escrito el código, así que...

Para lo que vale, este es el mismo* el ajuste del modelo con rstan. Se estima en 11 segundos en mi portátil consumer, el logro de un mayor tamaño de muestra efectiva de nuestros parámetros de interés $(N, \theta)$ en menos iteraciones.

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))

Tenga en cuenta que yo echo theta como un 2-simplex. Esto es solo para la estabilidad numérica. La cantidad de interés es theta[1]; obviamente theta[2] es información superflua.

*Como se puede ver, el posterior resumen es prácticamente idéntico, y la promoción de la $N$ a una cantidad real no parecen tener un impacto sustantivo en nuestras inferencias.

El 97.5% cuantil $N$ es 50% más grande para mi modelo, pero creo que es porque stan sampler es mejor para explorar la gama completa de la parte posterior de un simple paseo aleatorio, por lo que es más fácil hacerlo en las colas. Puedo estar equivocado, sin embargo.

            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

Tomando los valores de $N, \theta$ generados a partir de stan, yo uso estos para dibujar posterior valores predictivos $\tilde{y}$. No debemos sorprendernos de que la media de la parte posterior de las predicciones $\tilde{y}$ está muy cerca de la media de los datos de la muestra!

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 

Para comprobar si la rstan sampler es un problema o no, me calcula la parte posterior a través de una red. Podemos ver que la parte posterior es en forma de banano; este tipo de posterior puede ser problemático para la métrica euclidiana HMC. Pero vamos a ver los resultados numéricos. (La gravedad de la banana, la forma es en realidad suprimida aquí desde $N$ es en la escala logarítmica.) Si usted piensa acerca de la forma de plátano por un minuto, te das cuenta de que debía de estar en la línea de $\bar{y}=\theta N$.

El código de abajo, se puede confirmar que los resultados de stan sentido.

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

Hm. Esto no es exactamente lo que yo habría esperado. Rejilla de evaluación para el 97.5% de los cuantiles es más cercana a las EXIGENCIAS de los resultados de la rstan de resultados. Al mismo tiempo, no creo que los resultados de la cuadrícula debe ser tomado como el evangelio de la causa de la cuadrícula de la evaluación es hacer varias bastante gruesa simplificaciones: resolución de la rejilla no es demasiado fino, por un lado, mientras que por el otro, estamos (falsamente) la afirmación de que el total de la probabilidad en la red debe ser 1, ya que debemos llegar a un límite (y finito de puntos de cuadrícula) para que el problema sea computable (todavía estoy a la espera infinita de RAM). En verdad, nuestro modelo tiene probabilidad positiva sobre $(0,1)\times\left\{N|N\in\mathbb{Z}\land N\ge72)\right\}$. Pero tal vez algo más sutil que está en juego aquí.

Answer 2

Gracias de nuevo a @StéphaneLaurent y @user777 por su valioso aporte en los comentarios. Después de algunos ajustes de la previa para $\lambda$ ahora puedo replicar los resultados del papel de Raftery (1988).

Aquí está mi análisis de la secuencia de comandos y resultados utilizando PUNTAS y 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")

Cálculo tomó alrededor de 98 segundos en mi pc de escritorio.

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

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

Los resultados son:

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 parte posterior de la media de $N$ $522$ y la parte posterior de la mediana es $233$. He calculado el posterior modo de $N$ en el registro de escala y de nuevo transformado en la estimación:

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

Y el 80%-HPD de $N$ es:

(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

La posterior modo de $N$ $150$ y el 80%-el DEPARTAMENTO es $(78; 578)$, que está muy cerca de los límites indicados en el papel que se $(80; 598)$.