Chapter 7.5

Question 4

library(tidyverse)
library(purrrfect)

(sim_study <- parameters(~n, ~mu, ~sigma,
                         c(4,8,15), c(-2,0,2), c(1,2,3))
  %>% add_trials(10000)
  %>% mutate(ysample = pmap(list(n, mu, sigma), \(n,m,s) rnorm(n, mean = m, sd = s)))
  %>% mutate(ybar = map_dbl(ysample, mean),
             s = map_dbl(ysample, sd))
  %>% mutate(simulated = (ybar - mu)/(s/sqrt(n)))
  %>% mutate(analytic = dt(simulated, df = n-1))
)

# A tibble: 270,000 × 9
       n    mu sigma .trial ysample     ybar     s simulated analytic
   <dbl> <dbl> <dbl>  <dbl> <list>     <dbl> <dbl>     <dbl>    <dbl>
 1     4    -2     1      1 <dbl [4]> -2.49  0.568    -1.71    0.0941
 2     4    -2     1      2 <dbl [4]> -2.15  1.04     -0.281   0.349 
 3     4    -2     1      3 <dbl [4]> -1.76  1.20      0.409   0.330 
 4     4    -2     1      4 <dbl [4]> -0.889 0.612     3.63    0.0126
 5     4    -2     1      5 <dbl [4]> -1.21  1.20      1.32    0.147 
 6     4    -2     1      6 <dbl [4]> -2.08  0.741    -0.218   0.356 
 7     4    -2     1      7 <dbl [4]> -2.13  0.506    -0.525   0.308 
 8     4    -2     1      8 <dbl [4]> -2.90  0.829    -2.17    0.0555
 9     4    -2     1      9 <dbl [4]> -1.34  1.34      0.983   0.210 
10     4    -2     1     10 <dbl [4]> -1.86  1.24      0.221   0.356 
# ℹ 269,990 more rows

(ggplot(data = sim_study, aes(x = simulated))
 + geom_histogram(aes(y = after_stat(density)),
                  fill = 'goldenrod',
                  center = 0.1, binwidth = 0.2)
 + geom_line(aes(y = analytic), col = 'cornflowerblue')
 + facet_grid(mu~sigma~n, labeller = label_both, scale = 'free_x')
 + labs(title = 'Verifying t distribution')
 + theme_classic()
 )