# Monte Carlo simulation out <- t.test(extra ~ group, data = sleep, var.equal = TRUE) names(out) out[["p.value"]] out[["p.value"]] < 0.05 set.seed(127777) y1 <- rnorm(20, 10, 2) y2 <- rnorm(20, 10, 2) y <- c(y1, y2) g <- rep(1:2, each = 20) dat <- data.frame(g, y) out <- t.test(y ~ g, data = dat, var.equal = TRUE) out[["p.value"]] < 0.05 for(i in 1:1000) { print(i) } for(i in c("orange", "apple", "coconut", "banana")) { print(i) } for(i in seq(1, 1000, 129)) { print(i) } for(i in 1:50) { print(mean(rnorm(3, 0, 1))) } par(mfrow = c(3, 1)) result <- rep(NA, 1000) for(i in 1:1000) { result[i] <- mean(rnorm(100, 0, 1)) } hist(result, xlim = c(-5, 5), ylim = c(0,1), prob = TRUE) set.seed(123321) pvalues <- rep(NA, 10000) for(i in 1:10000) { y1 <- rnorm(20, 10, 2) y2 <- rnorm(20, 10, 2) y <- c(y1, y2) g <- rep(1:2, each = 20) dat <- data.frame(g, y) out <- t.test(y ~ g, data = dat, var.equal = TRUE) result[i] <- out[["p.value"]] } mean(result < 0.05) set.seed(123321) pvalues <- rep(NA, 10000) ns <- c(100, 20) vs <- c(1, 10) for(i in 1:10000) { y1 <- rnorm(ns[1], 10, sqrt(vs[1])) y2 <- rnorm(ns[2], 10, sqrt(vs[2])) y <- c(y1, y2) g <- rep(1:2, ns) dat <- data.frame(g, y) out <- t.test(y ~ g, data = dat, var.equal = TRUE) result[i] <- out[["p.value"]] } mean(result < 0.05) ############## sapply set.seed(123321) ns <- c(100, 20) vs <- c(1, 10) dat.l <- list() for(i in 1:1000) { y1 <- rnorm(ns[1], 10, sqrt(vs[1])) y2 <- rnorm(ns[2], 10, sqrt(vs[2])) y <- c(y1, y2) g <- rep(1:2, ns) dat.l[[i]] <- data.frame(g, y) } out.l <- list() for(i in 1:1000) { out.l[[i]] <- t.test(y ~ g, data = dat.l[[i]], var.equal = TRUE) } pvalues <- rep(NA, 1000) for(i in 1:1000) { pvalues[i] <- out.l[[i]][["p.value"]] } mean(pvalues < 0.05) nrow(dat.l[[3]]) sapply(dat.l, nrow) aggregate(y ~ g, data = dat.l[[3]], FUN = mean)[,2] sapply(dat.l, function(u) aggregate(y ~ g, data = u, FUN = mean)[,2]) tempfun <- function(u) { aggregate(y ~ g, data = u, FUN = mean)[,2] } sapply(dat.l, tempfun) tempfun2 <- function(u) { u[["p.value"]] } sapply(out.l, tempfun2) sapply(out.l, function(u) u[["p.value"]]) sapply(out.l, "[[", "p.value") lapply(out.l, "[[", "p.value") lapply(dat.l, function(u) tapply(u$y, u$g, mean)) lapply(dat.l, function(u) aggregate(y ~ g, u, FUN = mean)) lapply(dat.l, cor) analysisfun <- function(dat) { out <- t.test(y ~ g, data = dat, var.equal = TRUE) } lapply(dat.l, analysisfun) ################################### set.seed(123321) pvalues <- rep(NA, 10000) ns <- c(100, 20) vs <- c(1, 10) for(i in 1:10000) { y1 <- rnorm(ns[1], 10, sqrt(vs[1])) y2 <- rnorm(ns[2], 10, sqrt(vs[2])) y <- c(y1, y2) g <- rep(1:2, ns) dat <- data.frame(g, y) out <- t.test(y ~ g, data = dat, var.equal = TRUE) result[i] <- out[["p.value"]] } mean(result < 0.05) gendata <- function(index, ns, vs) { y1 <- rnorm(ns[1], 10, sqrt(vs[1])) y2 <- rnorm(ns[2], 10, sqrt(vs[2])) y <- c(y1, y2) g <- rep(1:2, ns) data.frame(g, y) } analysisdata <- function(dat) { t.test(y ~ g, data = dat, var.equal = TRUE) } set.seed(123321) dat.l <- lapply(1:1000, gendata, ns = c(100, 20), vs = c(1, 10)) out.l <- lapply(dat.l, analysisdata) pvalues <- sapply(out.l, "[[", "p.value") mean(pvalues < 0.05) ############# vec <- rbinom(100, 1, 0.2) mean(vec) out <- chisq.test(table(vec), p = c(0.8, 0.2)) out[["p.value"]] ### Parallel Processing library(MASS) ns <- rep(100, 1000) S <- matrix(c(1, 0.2, 0.3, 0.2, 1, 0.5, 0.3, 0.5, 1), 3, 3, byrow = TRUE) mvrnorm(100, mu = c(0, 0, 0), Sigma = S) dat2.l <- lapply(ns, mvrnorm, mu = c(0, 0, 0), Sigma = S) tempfun <- function(u) { colnames(u) <- c("y", "x1", "x2") return(data.frame(u)) } dat2.l <- lapply(dat2.l, tempfun) output.l <- lapply(dat2.l, lm, formula = y ~ x1 + x2) library(parallel) detectCores() # Function for Macs and Linux output.l <- mclapply(dat2.l, lm, formula = y ~ x1 + x2, mc.cores = 2) # Function for Windows cl <- makeCluster(rep("localhost", 2), type = "PSOCK") output.l <- parLapply(cl, dat2.l, lm, formula = y ~ x1 + x2) stopCluster(cl)