извлечение (модель, что = di c) из модели JAGS возвращает NA для штрафа - PullRequest
Новая
       94

извлечение (модель, что = di c) из модели JAGS возвращает NA для штрафа

0 голосов
/ 04 марта 2020

Используя JAGS, я подгоняю различные модели к данным и хотел бы сравнить их соответствия, используя информационный критерий отклонения (DI C). Я использую "run.jags", чтобы соответствовать модели, а затем "извлекать", чтобы определить DI C для модели после ее запуска. Мои модели сходятся без проблем, но я получаю значения только для части отклонения DI C. Все мои штрафные значения равны 0 или NA. Я думаю, я понимаю, почему я получаю NA - это сценарий ios, где прогнозируемое значение и наблюдаемое значение равны 0. Я не понимаю, почему я получаю 0 для других экземпляров. Есть какие нибудь идеи как это починить?

Другие посты, где кто-то получал NA за штраф, предлагали изменить приоры (https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/2fcd66ea/), но они использовали di c .samples (), а не extract (). Я попытался изменить свои приоры, но не нашел, что это изменило мой результат.

Вот код, который воспроизводит ситуацию (время выполнения <1 мин): </p> 

# 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")
...