帶收斂測試的並行RJAGS

Question

我正在使用RJAGS修改現有模型。我想平行運行鏈，偶爾檢查Gelman-Rubin收斂診斷，看看是否需要繼續跑步。問題是，如果需要根據診斷值繼續運行，則重新編譯的鏈將從第一個初始化的先前值重新啓動，而不是在鏈停止的參數空間中的位置。如果我不重新編譯模型，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 = endstate（適用於相應的鏈）的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/

馬特