rnorm2 <- function(n,mean,sd) { mean + sd * scale(rnorm(n)) } p1 <- rnorm2 (100000000, mean=70, sd=10) a1 <- qnorm(.95) a2 <- qnorm(.975) a3 <- qnorm(.995) s.size <- 16 s1 <- sample(p1, s.size, replace=T) m.s1 <- mean(s1) sd.s1 <- sd(s1) se <- sd.s1/sqrt(s.size) se2 <- se*2 se3 <- se*3 m.s1 - se2 m.s1 + se2 # or se2 <- se*a2 se3 <- se*a3 m.s1 - se2 m.s1 + se2