Skip to content

Instantly share code, notes, and snippets.

@jepusto
Last active December 31, 2015 23:28
Show Gist options
  • Save jepusto/8059893 to your computer and use it in GitHub Desktop.
Save jepusto/8059893 to your computer and use it in GitHub Desktop.
An adaptation of the Welch t-test simulation, for running in parallel on the Stampede server of the Texas Advanced Computing Cluster.
#----------------------------------------------
# data-generating model
#----------------------------------------------
two_group_data <- function(iterations, n, p, var_ratio, delta) {
Group <- c(rep("C", n * p), rep("T", n * (1 - p)))
Y_C <- matrix(rnorm(iterations * n * p, mean = 0, sd = 1), n * p, iterations)
Y_T <- matrix(rnorm(iterations * n * (1 - p), mean = delta, sd = sqrt(var_ratio)), n * (1 - p), iterations)
dat <- data.frame(Group, rbind(Y_C, Y_T))
return(dat)
}
#----------------------------------------------
# estimation procedures
#----------------------------------------------
CI_welch <- function(dat, alpha = 0.05) {
CI <- apply(dat[,-1], 2, function(X) t.test(X ~ dat$Group, var.equal=FALSE, conf.level = 1 - alpha)$conf.int)
return(t(CI))
}
#----------------------------------------------
# performance statistics
#----------------------------------------------
coverage <- function(CI, delta) {
covered <- (CI[,1] < delta) & (delta < CI[,2])
return(mean(covered))
}
#----------------------------------------------
# Simulation driver
#----------------------------------------------
run_sim <- function(iterations, n, p, var_ratio, delta, seed=NULL) {
if (!is.null(seed)) set.seed(seed)
dat <- two_group_data(iterations, n, p, var_ratio, delta)
Welch <- coverage(CI_welch(dat), delta)
return(c(Welch = Welch))
}
source_func <- ls()
#------------------------------------
# Design parameters
#------------------------------------
iterations <- 1000
n <- 12 * (1:5)
p <- 1 / c(2,3,4,6)
R <- 2^(-2:2)
delta <- 0
parms <- expand.grid(iterations=iterations, n=n, p=p, var_ratio=R, delta=delta)
parms$seed <- round(runif(nrow(parms)) * 2^30)
#--------------------------------------
# run simulations in parallel
#--------------------------------------
library(Rmpi)
library(snow)
library(foreach)
library(iterators)
library(doSNOW)
library(plyr)
# set up parallel processing
cluster <- getMPIcluster()
registerDoSNOW(cluster)
# export source functions
clusterExport(cluster, source_func)
# execute simulations
BFtime <- system.time(BFresults <- mdply(parms, .fun = run_sim, .parallel=TRUE))
stopCluster(cluster)
print(BFtime)
save(BFresults, BFtime, file="BFresults.Rdata")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment