How different from normal are the numbers generated by MASS::mvrnorm( empirical = TRUE)?
suppressMessages(library(tidyverse))
sim <- expand_grid(empirical=c(TRUE, FALSE),
n_sample = c(2, 4, 10, 20)) %>%
mutate(data = map2(empirical, n_sample, ~rerun(5000, enframe(MASS::mvrnorm(n = .y, rep(0,1), 1, tol = 1e-6, empirical = .x)[1,1], name = "n_obs")) %>%
bind_rows())) %>%
unnest(data)
## Plot
sim %>%
mutate(empirical = paste("Empirical =", empirical),
n_sample = paste("Sample size:", n_sample) %>%
factor(levels = paste("Sample size:", c(2,4,10,20)))) %>%
filter(n_obs==1) %>%
ggplot(aes(x=value)) +
geom_density() +
facet_grid(n_sample~empirical, scales="free_y")