Параллельные RJAGS с тестированием сходимости

Question

Я модифицирую существующую модель, используя RJAGS. Я хотел бы запускать цепочки параллельно и периодически проверять диагностику сходимости Гельмана-Рубина, чтобы узнать, нужно ли мне продолжать работать. Проблема в том, что если мне нужно возобновить работу на основе диагностического значения, перекомпилированные цепочки перезапускаются с первых инициализированных предыдущих значений, а не с места в пространстве параметров, где цепочка остановилась. Если я не перекомпилирую модель, RJAGS жалуется. Есть ли способ сохранить позиции цепочек, когда они останавливаются, чтобы я мог повторно инициализировать с того места, где остановился? Здесь я приведу очень упрощенный пример.

example1.bug:

model {
  for (i in 1:N) {
      x[i] ~ dnorm(mu,tau)
  }
  mu ~ dnorm(0,0.0001)
  tau <- pow(sigma,-2)
  sigma ~ dunif(0,100)
}

parallel_test.R:

#Make some fake data
N <- 1000
x <- rnorm(N,0,5)
write.table(x,
        file='example1.data',
        row.names=FALSE,
        col.names=FALSE)

library('rjags')
library('doParallel')
library('random')

nchains <- 4
c1 <- makeCluster(nchains)
registerDoParallel(c1)

jags=list()
for (i in 1:getDoParWorkers()){
  jags[[i]] <- jags.model('example1.bug',
                          data=list('x'=x,'N'=N))
}

# Function to combine multiple mcmc lists into a single one
mcmc.combine <- function( ... ){
  return( as.mcmc.list( sapply( list( ... ),mcmc ) ) )
}

#Start with some burn-in
jags.parsamples <- foreach( i=1:getDoParWorkers(),
                           .inorder=FALSE,
                           .packages=c('rjags','random'),
                           .combine='mcmc.combine',
                           .multicombine=TRUE) %dopar%
{
  jags[[i]]$recompile()

  update(jags[[i]],100)
  jags.samples <- coda.samples(jags[[i]],c('mu','tau'),100)

  return(jags.samples)
}   

#Check the diagnostic output
print(gelman.diag(jags.parsamples[,'mu']))

counter <- 0

#my model doesn't converge so quickly, so let's simulate doing
#this updating 5 times:
#while(gelman.diag(jags.parsamples[,'mu'])[[1]][[2]] > 1.04)
while(counter < 5)
{
counter <- counter + 1
jags.parsamples <- foreach(i=1:getDoParWorkers(),
                             .inorder=FALSE,
                             .packages=c('rjags','random'),
                             .combine='mcmc.combine',
                             .multicombine=TRUE) %dopar%
  {
    #Here I lose the progress I've made
    jags[[i]]$recompile()
    jags.samples <- coda.samples(jags[[i]],c('mu','tau'),100)
    return(jags.samples)
  }
}

print(gelman.diag(jags.parsamples[,'mu']))
print(summary(jags.parsamples))
stopCluster(c1)

В выводе я вижу:

Iterations = 1001:2000

где я знаю, должно быть> 5000 итераций. (кросс-пост на stats.stackexchange.com, который может быть более подходящим местом)

Answer 1

Каждый раз, когда ваша модель JAGS запускается на рабочих узлах, возвращаются образцы кода, но состояние модели теряется. Поэтому в следующий раз, когда он перекомпилируется, он перезапускается с самого начала, как вы видите. Чтобы обойти это, вам нужно получить и вернуть состояние модели в вашей функции (на рабочих узлах) следующим образом:

 endstate <- jags[[i]]$state(internal=TRUE)

Затем вам нужно передать это обратно на рабочий узел и заново сгенерировать модель в рабочей функции, используя jags.model() с inits=endtate (для соответствующей цепочки).

Я бы порекомендовал посмотреть на пакет runjags, который сделает все это за вас. Например:

library('runjags')
parsamples <- run.jags('example1.bug', data=list('x'=x,'N'=N), monitor=c('mu','tau'), sample=100, method='rjparallel')
summary(parsamples)
newparsamples <- extend.jags(parsamples, sample=100)
summary(parsamples)
# etc

Или даже:

parsamples <- autorun.jags('example1.bug', data=list('x'=x,'N'=N), monitor=c('mu','tau'), method='rjparallel')

Надеемся, что версия 2 runjags скоро будет загружена в CRAN, но сейчас вы можете скачать двоичные файлы по адресу: https://sourceforge.net/projects/runjags/files/runjags/

Matt