0

我正在使用nimble包来使用并行计算运行 MCMC 链。nimbleMCMCNimble 通常对于使用集成功能或runMCMC组件功能按顺序运行的每个链都有很好的进度条。

nimble您需要创建一个包含所有模型代码的自定义函数，包括该函数runMCMC，然后parLapply用于在集群上并行运行不同链的自定义函数。

https://r-nimble.org/nimbleExamples/parallelizing_NIMBLE.html

这计算得很好，但它会丢失原始进度条并且不显示进度。

这里有解决方案，但我认为它们不适用于 Nimble 框架：如何在 R 中显示并行计算中的代码进度？ 可以使用 来创建进度条pblapply，但这会将进度计算为基于已完成链/总链的百分比。例如，如果在 8 个内核上并行运行 8 个链，这基本上是没有信息的，因为它保持在 0%，直到所有链同时并行完成并跳转到 100%。

有没有办法从runMCMC自定义函数中打印原始进度条？或者，以其他方式显示进度？

这是一个玩具 CJS 隐藏马尔可夫模型来说明：

library(nimble)
library(parallel)
library(pbapply)

ymat <- matrix(sample(c(0,1),120,replace = TRUE),ncol = 6, nrow = 20)
y <- ymat

detectCores() # Number of cores on the machine = 8
this_cluster <- makeCluster(4) 

index <- c(1:4)

paraChain_code <- function(index, y){

library(nimble)  
first <- apply(y, 1, function(x) min(which(x!= 0),na.rm = TRUE))

hmm.model<- nimbleCode({
 
  # Priors
  phi ~ dunif(0,1)    # prior for survival φ
  p ~ dunif(0,1)      # prior for detection p
  
  # Vector of initial state probabilities δ
  delta[1] <- 1     # Alive
  delta[2] <- 0     # Dead
  
  # Transition matrix Γ
  gamma[1,1] <- phi      # transition Pr(alive t --> alive t+1)
  gamma[1,2] <- 1-phi    # transition Pr(alive t --> dead t+1)
  gamma[2,1] <- 0        # transition Pr(dead t --> alive t+1)
  gamma[2,2] <- 1        # transition Pr(dead t --> dead t+1)
  
  # Observation matrix Ω
  omega[1,1] <- 1-p   # observation Pr(alive t --> not detected t)
  omega[1,2] <- p     # observation Pr(alive t --> detected t)
  omega[2,1] <- 1     # observation Pr(dead t --> not detected t)
  omega[2,2] <- 0     # observation Pr(dead t --> detected t)
  
  # Likelihood
  for(i in 1:N){    # Loops over each individual (row)
    z[i,first[i]] ~ dcat(delta[1:2])  # Draws an initial state   
    for(j in (first[i]+1):T){    # Loops through each time point (column)
      z[i,j] ~ dcat(gamma[z[i,j-1], 1:2]) # Draws current state based on previous state
      y[i,j] ~ dcat(omega[z[i,j], 1:2])  # Draws an observation based on current state.
    }
  }
})

# Defines constants in the model
my.constants <- list(N = nrow(y),   # N individual
                     T = ncol(y),   # N of time points
                     first = first) # occasions of first capture

# Defines data in the model
# To remove 0s. 1 is non-detected, 2 is detected
my.data <- list(y = y+1)

# Defines initial values for states and probabilities
zinits <- y + 1
zinits[zinits == 2] <- 1
zinits[is.na(zinits) == TRUE] <- 1

initial.values <-  list(phi = runif(1,0,1),
                        p = runif(1,0,1),
                        z = zinits)

# Specfies the parameters of interest to record in output
parameters.to.save <- c("phi","p","z")

#MCMC model details
n.iter <- 1000
n.burnin <- 200

myModel <- nimbleModel(code = hmm.model,        
                               constants = my.constants,      
                               data = my.data,
                               inits = initial.values)
                               #monitors = parameters.to.save)
                               #niter = n.iter,
                               #nburnin = n.burnin)
                               #nchains = n.chains,
                               
CmyModel <- compileNimble(myModel)

myMCMC <- buildMCMC(CmyModel, monitors = parameters.to.save, enableWAIC = TRUE)

CmyMCMC <- compileNimble(myMCMC)

### This is the Nimble function that normally produces the progress bar. ###
results <- runMCMC(CmyMCMC, niter = n.iter, nburnin = n.burnin, WAIC = TRUE)
############################################################################

return(results)

}

# Function works correctly for multiple chains with parLapply
system.time({
chain_output <- parLapply(cl = this_cluster, index,
                          fun = paraChain_code,
                          y = ymat)
})

# The pblapply function includes a progress bar, but only for completed chains/total chains.
system.time({
  chain_output <- pblapply(X = index, cl = this_cluster, 
                            FUN = paraChain_code,
                            y = ymat)
 })

stopCluster(this_cluster)
4

0 回答 0