CLT simulation

[bensoltoff]

library(tidyverse)
set.seed(1234)

samp_n <- 1000

# draw population data
mu <- 168
pop <- rnorm(n = 1e06, mean = mu, sd = 1)
head(pop)
#> [1] 166.7929 168.2774 169.0844 165.6543 168.4291 168.5061

# draw a single sample
samp_1 <- sample(pop, size = samp_n)

(samp_1_mean <- mean(samp_1))
#> [1] 167.9602
(samp_1_sd <- sd(samp_1))
#> [1] 0.9835359

# 95% CI
(samp_1_CI <- qnorm(c(.025, .975), mean = samp_1_mean, sd = samp_1_sd / sqrt(samp_n)))
#> [1] 167.8992 168.0211

# draw 1000 samples and calculate the mean
samples <- data_frame(samp = rerun(1000, sample(pop, size = samp_n))) %>%
  mutate(mean = map_dbl(samp, mean),
         se = map_dbl(samp, sd) / sqrt(samp_n),
         CI = map2(mean, se, ~ qnorm(c(.025, .975), mean = .x, sd = .y)),
         pop_in_CI = map_lgl(CI, ~ mu >= .x[[1]] & mu <= .x[[2]]),
         mean_in_samp1 = map_lgl(mean, ~ .x >= samp_1_CI[[1]] & .x <= samp_1_CI[[2]]))
samples
#> # A tibble: 1,000 x 6
#>    samp           mean     se CI        pop_in_CI mean_in_samp1
#>    <list>        <dbl>  <dbl> <list>    <lgl>     <lgl>        
#>  1 <dbl [1,000]>  168. 0.0320 <dbl [2]> TRUE      TRUE         
#>  2 <dbl [1,000]>  168. 0.0314 <dbl [2]> TRUE      TRUE         
#>  3 <dbl [1,000]>  168. 0.0315 <dbl [2]> TRUE      TRUE         
#>  4 <dbl [1,000]>  168. 0.0322 <dbl [2]> TRUE      TRUE         
#>  5 <dbl [1,000]>  168. 0.0313 <dbl [2]> TRUE      FALSE        
#>  6 <dbl [1,000]>  168. 0.0312 <dbl [2]> TRUE      TRUE         
#>  7 <dbl [1,000]>  168. 0.0314 <dbl [2]> TRUE      TRUE         
#>  8 <dbl [1,000]>  168. 0.0325 <dbl [2]> TRUE      TRUE         
#>  9 <dbl [1,000]>  168. 0.0314 <dbl [2]> TRUE      TRUE         
#> 10 <dbl [1,000]>  168. 0.0324 <dbl [2]> TRUE      TRUE         
#> # ... with 990 more rows

# proportion of sample 95% CI which contain the population mean
mean(samples$pop_in_CI)
#> [1] 0.951

# proportion of sample means which fall in the 95% confidence interval for sample 1
mean(samples$mean_in_samp1)
#> [1] 0.759

Created on 2018-09-18 by the reprex
package (v0.2.0).