チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

Rでforeachをつかって楽にシミュレーションを分割してみる雛形

この雛形を使うと並列計算ができる。個人的には返り値などは基本的に全部書きだしてあとで分析する形にしておくと安心。
シミュレーションを関数にしてしまうことと,foreach関数の指定のときに必要なパッケージや環境の変数(自作の関数とか)を全部持っていてもらうように指定するのがポイント。

#
# シミュレーションの雛形
# 
rm(list=ls())
library("MASS")
sim_times = 1:100

simulation_fun <- function(sim_times){
  setwd("hogedir")
  for(i in c(500, 5000)){
    for(pop_cor in c(0, 0.5, 0.8)){
      #結果保存用フォルダ  
      dir.create(paste0("sim_",i,": pop_cor = ",pop_cor)) 
      
      for(n in sim_times){
        set.seed(i*10 + 100*pop_cor)
        cat(format(Sys.time(), "%Y %b%a %d %X"),  ": i = ",i,": pop_cor = ",pop_cor,": n = ",n,"\n")

        Sig <- matrix(pop_cor, 10, 10 )
        diag(Sig) <- 1
        mu = rep(0,K)
        
        x = MASS::mvrnorm(I,mu = mu, Sigma = Sig)
        
        time1 = as.POSIXlt(Sys.time())
        res <- list(mean=colMeans(x),cov = cov(x))
        time2 = as.POSIXlt(Sys.time())
        use.time = difftime(time2, time1, units="secs")

        #use.time
        save(list = c("res", "x","time1","time2","use.time"),file = paste0("sim_",i,"_",pop_cor,"/sim_",i,"_",pop_cor,"_",n,".Rdata")) 
      } # sim times 
    } # pop cor
  } # sample size
  return()
}


library("foreach")
library(doParallel)

package_list <- c("MASS")

#
# 指定したいコアの数
# 
n_div <- 4
# 
# いくつからいくつのシミュレーションを分割するか指定
# 
start_ind <- seq(from=1, to = max(sim_times),by = max(sim_times)/n_div)
end_ind <- seq(from=max(sim_times), to = 1,by = -max(sim_times)/n_div)
end_ind <- end_ind[order(end_ind )]
divide_indicator <- mapply(function(x,y)x:y, start_ind, end_ind,SIMPLIFY = F)
divide_indicator

cores = n_div 
cluster = makeCluster(cores)
registerDoParallel(cluster)

#
# 環境を全部持っていってもらう
# ls(envir=parent.frame())
# 
foreach (sim_times = divide_indicator, .combine=c,.packages = package_list,.export = ls(envir=parent.frame())) %dopar% {
  simulation_fun(sim_times)
}
stopCluster(cluster)

>||