The following code essentially redoes what the code up to the first conf_interval block does there. Which one is more clear may be debatable but it's shorter by

The following code essentially redoes what the code up to the first conf_interval block does there. Which one is more clear may be debatable but it's shorter by a factor of two and faster by a factor of ten (45 seconds vs 4 for me). sample_size <- 60 sample_meansB <- lapply(population_dataB, function(x){ t(apply(replicate(20000, sample(x, sample_size)), 2, function(x) c(sample_mean=mean(x), sample_sd=sd(x)))) }) lapply(sample_meansB, head) ## check first rows population_data_statsB <- lapply(population_dataB, function(x) c(population_mean=mean(x), population_sd=sd(x), n=length(x))) do.call(rbind, population_data_statsB) ## stats table cltB <- mapply(function(s, p) (s[,"sample_mean"]-p["population_mean"])/(p["population_sd"]/sqrt(sample_size)), sample_meansB, population_data_statsB) head(cltB) ## check first rows small_sample_size <- 6 repeated_samplesB <- lapply(population_dataB, function(x){ t(apply(replicate(10000, sample(x, small_sample_size)), 2, function(x) c(sample_mean=mean(x), sample_sd=sd(x)))) }) conf_intervalsB <- lapply(repeated_samplesB, function(x){ sapply(c(lower=0.025, upper=0.975), function(q){ x[,"sample_mean"]+qnorm(q)*x[,"sample_sd"]/sqrt(small_sample_size) })}) within_ci <- mapply(function(ci, p) (p["population_mean"]>ci[,"lower"]&p["population_mean"]<ci[,"upper"]), conf_intervalsB, population_data_statsB) apply(within_ci, 2, mean) ## coverage One can do simple plots similar to the ones in that page as follows: par(mfrow=c(2,3), mex=0.8) for (d in colnames(population_dataB)) plot(density(population_dataB[,d], bw="SJ"), main=d, ylab="", xlab="", las=1, bty="n") for (d in colnames(cltB)) plot(density(cltB[,d], bw="SJ"), main=d, ylab="", xlab="", las=1, bty="n") for (d in colnames(cltB)) { qqnorm(cltB[,d], main=d, ylab="", xlab="", las=1, bty="n"); qqline(cltB[,d], col="red") }