##### ##### NMST440: Advanced Aspects of the R Environment ##### ##### --------------------------------------------------------------------- ##### ##### Simulation script towards distribution ##### of a test statistic of a certain test of independence. ##### ##### --------------------------------------------------------------------- ##### ##### Arnošt Komárek ##### https://www2.karlin.mff.cuni.cz/~komarek ##### komarek@karlin.mff.cuni.cz ##### ##### ====================================================================== args <- commandArgs(TRUE) #args <- c("1", "0.5", "20") a <- as.numeric(args[1]) rho <- as.numeric(args[2]) n <- as.numeric(args[3]) ### For real simulation on cluster # ROOT <- "/usr/work/komarek/NMST440/" # SRCDIR <- ROOT # RESDIR <- paste(ROOT, "SimulIndTest/", sep = "") # Mmax <- 10000 ## should be a multiple of 8 ### For developing purposes on a local PC # a <- 1; rho <- 0.5; n <- 20 ROOT <- "/home/komarek/teach/mff_2021/nmst440_AdvRko/Tutorial11/" SRCDIR <- "/home/komarek/teach/mff_2021/nmst440_AdvRko/Tutorial11/" RESDIR <- paste(ROOT, "SimulIndTest/", sep = "") Mmax <- 1000 ### RESDIR is assumed to have subdirectories "Log" and "RData" ### Main name of the file to store results etc. NAME <- paste("I_", a*10, "_", rho*100, "_", n, sep = "") # print(NAME) ### Function to run simulation (M replicates) ### Returns: simulated values of the test statistic # nrun <- 1 simFun <- function(nrun = 1, Mmax = 100, a, rho, n, SRCDIR, RESDIR, NAME) { ### Load functions and compiled objects source(paste(SRCDIR, "Rfun/rMVN2.R", sep = "")) source(paste(SRCDIR, "Rfun/indTest.R", sep = "")) dyn.load(paste(SRCDIR, "C/indTest.so", sep = "")) ### Open log file sink(paste(RESDIR, "Log/", NAME, "-", nrun, ".log", sep = "")) cat("Started: ", date(), "\n\n", sep = "") ### Main simulation loop T <- as.numeric(rep(NA, Mmax)) M <- 0 while(M < Mmax){ ### Generate data xe <- rMVN2(n, mean = c(0, 0), sd = c(1, 1), rho = rho) colnames(xe) <- c("x", "e") ### Calculate value of the test statistic Tobj <- try(indTestC2(x = xe[, "x"], e = xe[, "e"], a = a)) if (class(Tobj)[1] == "try-error") next ### Keep value of the test statistic M <- M + 1 T[M] <- Tobj$T ### Print information on progress (at every 100th step) if (M %% 100 == 0) cat(" *** M = ", M, ": ", date(), "\n", sep = "") } ### Close log file cat("\nFinished: ", date(), "\n\n", sep = "") sink() ### Return simulated values of the test statistic return(T) } ### Test one run # simFun(nrun = 1, Mmax = 100, a = a, rho = rho, n = n, SRCDIR = SRCDIR, RESDIR = RESDIR, NAME = NAME) ### Run simulation (in parallel on 8 cores) ### On Sněhurka: we ask 4 cores for each simulation setting, ### since each Sněhurka core has 2 CPUs, we may ask for 8 cores library("parallel") myCluster <- makeCluster(8, type = "FORK") ## type = "FORK" must be specified, otherwise problems on the cluster ## with releasing slaves clusterSetRNGStream(myCluster, 221913282) Ts <- parLapply(myCluster, 1:8, simFun, M = Mmax/8, a = a, rho = rho, n = n, SRCDIR = SRCDIR, RESDIR = RESDIR, NAME = NAME) stopCluster(myCluster) ### Merge results T <- c(Ts[[1]], Ts[[2]], Ts[[3]], Ts[[4]], Ts[[5]], Ts[[6]], Ts[[7]], Ts[[8]]) ### Save results save(list = "T", file = paste(RESDIR, "RData/", NAME, ".RData", sep = ""))