Estimação bayesiana de

Esta pergunta é um acompanhamento técnico desta pergunta .

Tenho dificuldade em entender e replicar o modelo apresentado em Raftery (1988): Inferência para o binômio $N$ parâmetro : uma abordagem hierárquica de Bayes no WinBUGS / OpenBUGS / JAGS. Não se trata apenas de código, portanto, ele deve estar no tópico aqui.

fundo

Seja um conjunto de contagens de sucesso de uma distribuição binomial com N e θ desconhecidos . Além disso, suponho que N siga uma distribuição de Poisson com o parâmetro μ (como discutido no artigo). Então, cada x i tem uma distribuição de Poisson com média λ = μ θ . Quero especificar os anteriores em termos de λ e θ .$x=(x_{1},\ldots,x_{n})$$N$$\theta$$N$$\mu$$x_{i}$$\lambda = \mu \theta$$\lambda$$\theta$

Supondo que eu não tenha nenhum conhecimento prévio bom sobre ou θ , desejo atribuir anteriores não informativos a λ e θ . Por exemplo, os antecedentes são λ ~ L um m m um ( 0,001 , 0,001 ) e θ ~ L n i f o r m ( 0 , 1 ) .$N$$\theta$$\lambda$$\theta$$\lambda\sim \mathrm{Gamma}(0.001, 0.001)$$\theta\sim \mathrm{Uniform}(0, 1)$

O autor usa um prior inadequado de mas o WinBUGS não aceita antecedentes impróprios.$p(N,\theta)\propto N^{-1}$

Exemplo

No artigo (página 226), são fornecidas as seguintes contagens de sucesso dos waterbucks observados: . Quero estimar N , o tamanho da população.$53, 57, 66, 67, 72$$N$

Aqui está como eu tentei descobrir o exemplo no WinBUGS ( atualizado após o comentário 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)

O modelo não Sill não convergem bem após 500.000 foi amostras com 20.000 burn-in amostras. Aqui está a saída de uma execução 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

Questões

Claramente, estou perdendo alguma coisa, mas não consigo ver exatamente o que. Eu acho que minha formulação do modelo está errada em algum lugar. Então, minhas perguntas são:

• Por que meu modelo e sua implementação não funcionam?
• Como o modelo dado por Raftery (1988) pode ser formulado e implementado corretamente?

Obrigado pela ajuda.

Bem, como você conseguiu que seu código funcionasse, parece que essa resposta é um pouco tarde demais. Mas eu já escrevi o código, então ...

Pelo que vale, esse é o mesmo modelo * rstan. É estimado em 11 segundos no meu laptop consumidor, atingindo um tamanho de amostra efetivo mais alto para nossos parâmetros de interesse em menos iterações.$(N, \theta)$

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

Observe que eu faço a conversão thetacomo 2 simplex. Isso é apenas para estabilidade numérica. A quantidade de interesse é theta[1]; obviamentetheta[2] é informação supérflua.

* Como você pode ver, o resumo posterior é praticamente idêntico e promove $N$ para uma quantidade real não parece ter um impacto substancial em nossas inferências.

O quantil de 97,5% para é 50% maior para o meu modelo, mas acho que é porque o amostrador de stan é melhor em explorar toda a extensão da parte posterior do que uma simples caminhada aleatória, para que possa ser mais fácil entrar nas caudas. Eu posso estar errado, no entanto.$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

Tomando os valores de gerados a partir stan, eu usá-los para desenhar posteriores valores preditivos ~ y . Nós não devemos nos surpreender que média das previsões posteriores ~ y está muito próximo da média dos dados da amostra!$N, \theta$$\tilde{y}$$\tilde{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 

Para verificar se o rstanamostrador é um problema ou não, calculei a parte posterior sobre uma grade. Podemos ver que o posterior é em forma de banana; esse tipo de posterior pode ser problemático para o HMC métrico euclidiano. Mas vamos verificar os resultados numéricos. (A gravidade da forma de banana é realmente suprimida aqui desde é na escala log). Se você pensar sobre a forma de banana por um minuto, você vai perceber que ele deve estar na linha ˉ y = θ N .$N$$\bar{y}=\theta N$

O código abaixo pode confirmar que nossos resultados da stan fazem 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

rstan$(0,1)\times\left\{N|N\in\mathbb{Z}\land N\ge72)\right\}$

$\lambda$ , agora posso replicar os resultados do artigo de Raftery (1988).

Aqui está o meu script de análise e resultados usando JAGS e 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")

A computação levou cerca de 98 segundos no meu PC de mesa.

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

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

Os resultados são:

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

$N$$522$$233$$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

$N$

(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

$N$$150$$(78; 578)$$(80; 598)$