r-带收敛测试的并行RJAGS

Question

我正在使用RJAGS修改现有模型。我想并行运行链，偶尔检查Gelman-Rubin收敛诊断，看看是否需要继续运行。问题是，如果我需要根据诊断值恢复运行，重新编译的链会从第一个初始化的先前值重新启动，而不是从参数空间中链停止的位置重新启动。如果我不重新编译模型，RJAGS会抱怨。有没有一种方法可以在链条停止时存储它们的位置，这样我就可以从停止的地方重新初始化？这里我将举一个非常简单的例子。

示例1.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)
}

平行测试。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模型在worker节点上运行时，都会返回coda样本，但模型的状态会丢失。因此，下次重新编译时，它将从头开始重新启动，正如您所看到的那样。要绕过这一点，您需要在函数（在工作节点上）中获取并返回模型的状态，如下所示：

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

然后，您需要将其传递回worker节点，并使用inits=endstate的jags.model（）在worker函数中重新生成模型（对于适当的链）。

实际上，我建议你看看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')

runjags的第2版有望很快上传到CRAN，但目前您可以从以下位置下载二进制文件：https://sourceforge.net/projects/runjags/files/runjags/

Matt