在R中使用数据加速随机马尔可夫链.表或并行化

Question

我正在尝试使用数据加速离散时间非齐次马尔可夫链的蒙特卡罗模拟。表或某种形式的并行化。使用随机虚拟转换矩阵TM，我在每个N个模拟中模拟nSteps时间步，并从初始状态向量initialState开始记录currentState中的下一个更新状态。在每个时间步，I矩阵乘以当前状态与转移矩阵TM。

代码1与循环

nStates <- 5 #number of states
initialState <- c(rep(1/nStates, nStates)) #vector with uniform initial states
nSteps <- 10 #number of time steps
N <- 10000 #number of simulations
ind.arr <- matrix(1:(N*nSteps),ncol=nSteps, byrow=TRUE)
currentState <- vector("list",(N*(nSteps))) #collects the nSteps state vectors for each simulation
system.time(
  for (i in 1:N) {
    TM <- matrix(runif(nStates^2), ncol=nStates) #random transition matrix for each time step and each simulation
    currentState[[(ind.arr[i,1])]] <- initialState %*% (TM / rowSums(TM)) #/rowSums(TM) ensures that TM is a transition matrix
    for (t in 2:nSteps){
      TM <- matrix(runif(nStates^2), ncol=nStates)
      currentState[[(ind.arr[i,t])]] <- currentState[[(ind.arr[i,t-1])]] %*% (TM / rowSums(TM))
    }
  })

代码不是超级慢，但是我想知道避免N-loop是否可以加速代码。如果我把n循环的主体放到函数

中

statefun <- function(initialState, nSteps, nStates){
  TM <- matrix(runif(nStates^2), ncol=nStates) #random transition matrix for each time step and each simulation
  currentState <- matrix(rep(NA, nSteps*nStates), ncol=nStates)
  currentState[1,] <- initialState %*% (TM / rowSums(TM)) #/rowSums(TM) ensures that TM is a transition matrix
  for (t in 2:nSteps){
    TM <- matrix(runif(nStates^2), ncol=nStates)
    currentState[t,] <- currentState[t-1,] %*% (TM / rowSums(TM))
  }
  return(currentState)
}

和使用数据。表中，我得到一个错误，而不是期望的结果

library(data.table)
system.time(dt <- data.table(i=1:N)[, c("s1", "s2", "s3", "s4", "s5") := list(statefun(initialState, nSteps, nStates)), by=i])
#As each simulation run is independent and the call of statefun is expensive, I was hoping that parallelisation helps to accelerate the code, but trying foreach is actually slower than where I started.  
library(foreach)
system.time(res <- foreach(i=1:N, .combine='c') %do% statefun(initialState, nSteps, nStates))

我很欣赏任何关于如何制作数据的评论。表工作或在这种情况下使用并行化。非常感谢,蒂姆。

@编辑:这一个不拿起函数调用的十行输出…

system.time( #does not work 
  dt <- data.table(i=1:N)[,c("s1", "s2", "s3", "s4", "s5"):=as.list(statefun(initialState, nSteps, nStates)),by=i]
)

Answer 1

如果将外部for循环转换为具有10,000个任务的foreach循环，则性能不会很好，因为任务太小了。让任务的数量与工人的数量相等通常会更好。这里有一个简单的方法，使用iterators包中的idiv函数:

library(doParallel)
nw <- 4
cl <- makePSOCKcluster(nw)
registerDoParallel(cl)
nStates <- 5
initialState <- c(rep(1/nStates, nStates))
nSteps <- 10
N <- 10000
currentState <- foreach(n=idiv(N, chunks=nw), .combine='c') %dopar% {
  ind.arr <- matrix(1:(n * nSteps), ncol=nSteps, byrow=TRUE)
  cur <- vector("list", n * nSteps)
  for (i in 1:n) {
    TM <- matrix(runif(nStates^2), ncol=nStates)
    cur[[ind.arr[i,1]]] <- initialState %*% (TM / rowSums(TM))
    for (t in 2:nSteps) {
      TM <- matrix(runif(nStates^2), ncol=nStates)
      cur[[(ind.arr[i,t])]] <-
          cur[[(ind.arr[i,t-1])]] %*% (TM / rowSums(TM))
    }
  }
  cur
}

不是简单地并行化外部for循环，而是在较小版本的顺序代码周围添加foreach循环。因此，如果你找到一种改进顺序代码的方法，你就可以很容易地在并行版本中使用它。如果不返回所有的中间状态，还可以获得更好的性能。

Answer 2

这个线程中有一个例子可以满足您的需求。您需要使用复制，从base.

中的lapply函数。

 replicate(N, statefun(initialState, nSteps, nStates))