0

使用 JAGS，我将不同的模型拟合到数据中，并希望使用偏差信息标准 (DIC) 来比较它们的拟合度。我正在使用“run.jags”来拟合模型，然后在模型运行后“提取”来确定模型的 DIC。我的模型没有问题地收敛，但我只获得了 DIC 偏差部分的值。我所有的惩罚值都是 0 或 NA。我想我明白为什么我得到 NA - 这些是预测值和观察值都为 0 的场景。我不明白为什么我在其他实例中得到 0。有想法该怎么解决这个吗？

有人因处罚而获得 NA 的其他帖子建议更改先验（https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/2fcd66ea/），但他们使用的是 dic.samples()，而不是提炼（）。我尝试改变我的先验，但没有发现它改变了我的结果。

这是一些重现这种情况的代码（运行时间 < 1 分钟）：

# install.packages("remotes")
# remotes::install_github("gilesjohnr/hmob")

library(hmob) 
library(dplyr)
library(stringr)
library(foreach)
library(parallel)
library(doParallel)
library(zoo)
library(sp)
library(rgdal)
library(rgeos)
library(abind)
library(rjags)
library(coda)
library(runjags)
library(truncnorm)
library(rmutil)
library(dclone)
library(R2WinBUGS)

# subset of data to run
M <- matrix(c(0, 5514, 5290, 88, 5501, 0, 10868, 392, 5388, 10830, 0, 6641, 91, 400, 6660, 0),
   nrow = 4, ncol = 4, byrow = TRUE)

D <- matrix(c(0, 38, 58, 162, 38, 0, 35, 125, 58, 35, 0, 111, 162, 125, 111, 0),
   nrow = 4, ncol = 4, byrow = TRUE) 

N <-  c(15350, 17803, 29825, 5772) 

n.districts <- nrow(M)

jags.data <- list(
     M=M,                                          
     D=D,                                          
     N=N,                                           
     n.districts=n.districts)

# JAGS model
model.test <- "
model {

     for (i in 1:n.districts) {
          for (j in 1:n.districts) {
               M[i,j] ~ dpois(pi[i,j]*N[i])      
          }
          pi[i,1:n.districts] <- c[i,]/sum(c[i,]) 
     }

     for (i in 1:n.districts) {
          for (j in 1:n.districts) {
               c[i,j] <- ifelse(
                    i == j,
                    0,
                    exp(log(theta) + (omega.1*log(N[i]) + omega.2*log(N[j]) - log(f.d[i,j])))
               )
               f.d[i,j] <- D[i,j]^gamma
          }
     }

     ### Priors ###
     theta ~ dgamma(1, 1)
     omega.1 ~ dgamma(1, 1)
     omega.2 ~ dgamma(1, 1)
     gamma ~ dgamma(1, 1)

}"

params <- c('omega.1', 'omega.2', 'theta', 'gamma')

# Burnin and samples are intentionally low when troubleshooting
nc <- 4          # number of chains
na <- 1000       # adaptations
nb <- 4000       # burn in
ni <- 10000      # samples
nt <- 5          # thin

init.list <- replicate(nc,
                       list(.RNG.name='lecuyer::RngStream',
                            .RNG.seed= 423486),  #sample(1:1e6, 1)), uncomment for random sample
                       simplify=FALSE)

out <- run.jags(model=model.test,
                data=jags.data,
                monitor=params,
                n.chains=nc,
                adapt=na,
                burnin=nb,
                sample=ni,
                thin=nt,
                inits=init.list,
                modules=c('lecuyer'),
                method="parallel",
                summarise=FALSE)

dic.basic <- extract(out, what="dic")
4

0 回答 0