library(mvtnorm) makedat = function(n,mean,var,delta){ rmvnorm(n,c(mean,mean+delta),matrix(c(var,.6*var,.6*var,var),ncol=2)) } runsim = function(nsim,n,mean,var,delta){ nsig = 0 for(i in 1:nsim){ dat = makedat(n,mean,var,delta) if(t.test(dat[,1],dat[,2],paired=TRUE)$p.value < 0.05)nsig = nsig + 1 } nsig / nsim } print(system.time(res <- sapply(seq(20,200,by=10),function(n)runsim(1000,n,.3,.3^2,.1)))) library(multicore) print(system.time(res1 <- unlist(mclapply(seq(20,200,by=10),function(n)runsim(1000,n,.3,.3^2,.1)))))