7.5 code

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.3.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(purrrfect)

Attaching package: 'purrrfect'

The following objects are masked from 'package:base':

    replicate, tabulate

N <- 10000

(sim_grid <- parameters(~n, ~mu, ~sigma,
                       c(4, 8, 15),
                       c(-2, 0, 2),
                       c(1, 2, 3)) %>%
  add_trials(N) %>%
  mutate(ysample = pmap(list(n, mu, sigma),
                        \(nn, mm, ss) rnorm(nn, mean = mm, sd = ss))) %>%
  mutate(ybar = map_dbl(ysample, mean),
         s    = map_dbl(ysample, sd),
         T    = (ybar - mu) / (s / sqrt(n))))

# A tibble: 270,000 × 8
       n    mu sigma .trial ysample    ybar     s       T
   <dbl> <dbl> <dbl>  <dbl> <list>    <dbl> <dbl>   <dbl>
 1     4    -2     1      1 <dbl [4]> -2.13 1.55  -0.173 
 2     4    -2     1      2 <dbl [4]> -2.14 1.11  -0.260 
 3     4    -2     1      3 <dbl [4]> -1.71 0.683  0.838 
 4     4    -2     1      4 <dbl [4]> -1.97 1.00   0.0659
 5     4    -2     1      5 <dbl [4]> -1.31 0.996  1.39  
 6     4    -2     1      6 <dbl [4]> -1.98 0.895  0.0432
 7     4    -2     1      7 <dbl [4]> -1.73 0.964  0.558 
 8     4    -2     1      8 <dbl [4]> -2.38 1.31  -0.580 
 9     4    -2     1      9 <dbl [4]> -2.21 1.71  -0.247 
10     4    -2     1     10 <dbl [4]> -2.72 0.947 -1.52  
# ℹ 269,990 more rows

library(dplyr)
library(tidyr)
library(ggplot2)

t_curve <- sim_grid %>%
  distinct(n, mu, sigma) %>%
  mutate(df = n - 1) %>%
  crossing(T = seq(-6, 6, length.out = 400)) %>%
  mutate(dens = dt(T, df = df))

ggplot(sim_grid, aes(x = T)) +
  geom_density() +
  geom_line(data = t_curve, aes(x = T, y = dens), inherit.aes = FALSE) +
  facet_grid(mu + sigma ~ n) +
  coord_cartesian(xlim = c(-6,6)) +
  theme_classic()