nooutput=sapply(1:3,function(i){ n=10 mu=10 sigma=4 sem=sigma/sqrt(n) # Compute 1000 times the mean a random sample and get the standard deviation of the mean means_sd=sd(sapply(1:1000,function(dummy){mean(rnorm(n,mean=mu,sd=sigma))})) # See how SEM is a good approximation of the mean standard deviation message(paste("----- Experiment",i,"-----")) message(paste0("Means SD: ", round(means_sd,digits=2))) message(paste0("SEM ", round(sem,digits=2))) })