生成二项分布 N 的年度估计

Question

我正在尝试根据 https://stats.stackexchange.com/questions/113851/bayesian-estimation-of-n-of-a-binomial-distribution 中的二项分布估计 N。我有 40 年的计数和估计的成功率。我想生成每年 N 的估计值。我是 JAGS 的新手，我认为我可以通过添加索引 (n[i]) 来获得 N 的估计值，但我收到以下错误：“Index out of range taking subset of n.”我在这里所做的事情是否完全偏离了基地？下面是一个模型和随机数据集，我认为它会生成 N 的年度估计值。我确实有一个模型可以生成 N 的单个估计值，它只是从 n 中删除索引 ([i]) 的可能性 (dbin( θ, n)).预先感谢您的帮助。

sink("file.jags")
cat("
    model {

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

    ## Priors    
    lambda ~ dgamma(0.005, 0.005)    
    theta ~ dbeta(22, 37)
    mu <- lambda/theta
    n ~ dpois(mu)
    
}    
", fill = TRUE)

sink()

# Initial values
inits <- function() { 
  list(
    n = sample(1:100
               , size = 1) 
    , lambda = runif(1, 1, 100)
  )
}

# Parameters monitored
params <- c('theta'
            , 'n'
)

# MCMC settings
ni <- 2000
nb <- 200
nc <- 3

# Bundle data
jags.data <- list(n = c(round(runif(40, 1, 100)))
                  , nyear = 40)

# Fit model
out <- jagsUI(data = jags.data
              , inits = inits
              , parameters.to.save = params
              , model.file = "file.jags"
              , n.chains = nc
              , n.iter = ni
              , n.burnin = nb
              , verbose = TRUE)

# Summarize posteriors
print(out$summary, dig = 2)

编辑 1：正如@DaveArmstrong 所建议的，我需要在 for 循环中包含 n 的先验以生成年度先验。此外，我从初始值列表中删除了 n，因为我在“setParameters”中收到维度不匹配错误。一旦我从“inits”中删除 n，模型就会运行。

我还根据我对 Dave 的观点的解释将 lambda 上的先验从伽玛切换为制服：模型无法初始化。模型运行后，我检查了 lambda 上的不同先验是否有影响，但没有。该模型仍然使用“dgamma(0.005, 0.005)”运行，但我要先离开 dunif。

"model {

    # Likelihood    
    for (i in 1:nyear) {
      
      reported_morts[i] ~ dbin(report_rate, total_morts[i])
      
      # Annual prior for total mortality. Necessary to get annual estimates.
      total_morts[i] ~ dpois(mu)
    }

    # Priors
    lambda ~ dunif(1, 100)
    report_rate ~ dbeta(22, 37)
    mu <- round(lambda/report_rate)
    
}"  

# Initial values. Removed "n = sample(1:100, size = 1)".
inits <- function() { 
  list(
    lambda = runif(1, 1, 100)
  )
}
  

Answer 1

你应该能够通过将模型更改为这样的东西来做到这一点：

mod <- "model {

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

    ## Priors    
    lambda ~ dgamma(0.005, 0.005)    
    theta ~ dbeta(22, 37)
    mu <- lambda/theta
}    
"

之前的问题是

n

在先验中是标量，在模型中是向量。您必须先对

n

的每个值进行设置。该模型在语法上是正确的，但它不会初始化，因为您设置的先验与数据相距甚远。如果你只是得到

x

的先验预测分布，它们的后验均值都小于0.1，而

x

的数据是1到100之间的整数值。这会产生“节点与父节点不一致”的错误。我假设这些只是虚拟数据，所以也许你真正的问题会解决。否则，先验必须更现实一些才能工作。