sims.R

# 兩個獨立常態母體，4個參數 mu1 mu2 sigma1 sigma2 未知。各自隨機抽n1和n2個樣本，想知道(mu2 - mu1) / sigma1的抽樣分布。

library(magrittr)
library(data.table)
library(ggpubr)

# 單次模擬的function
sim <-
  function(mu1 = 10,
           mu2 = 10,
           sigma1 = 1,
           sigma2 = 1,
           n1 = 20,
           n2 = 20)
  {
    sample1 <- rnorm(n1, mu1, sigma1)
    sample2 <- rnorm(n2, mu2, sigma2)
    # stat是要求的統計量
    stat <- (mean(sample2) - mean(sample1)) / sd(sample1)
    stat
  }

# 同參數多次模擬的物件及其plot()與print()方法
sims <-
  function(nsim = 10000, mu1, mu2, sigma1, sigma2, n1, n2) {
    stat <- replicate(nsim, sim(mu1, mu2, sigma1, sigma2, n1, n2))
    class(stat) <- "sims"
    attr(stat, "pars") <- c(mu1, mu2, sigma1, sigma2, n1, n2)
    stat
  }
plot.sims <- function(x) {
  qplot(x %>% as.double, geom = "histogram", bins = 100) + 
    xlab("stat") + ylab("freq") + labs(subtitle = attr(x, "pars") %>% paste(collapse = ", "))
}
print.sims <- function(x) {
  c("xbar" = mean(x), "sd" = sd(x)) %>% print
  attr(x, "pars") %>% print
}

# 建立模擬的參數表
pars <- 
  expand.grid(
    mu1 = 50,
    mu2 = c(30, 50, 70),
    sigma1 = c(1, 5),
    sigma2 = c(1, 5),
    n1 = c(20, 50),
    n2 = c(20, 50)
  )
res <- vector("list", nrow(pars))
for(i in 1:length(res)) {
  res[[i]] <-
    sims(10000, pars[i, 1], pars[i, 2], pars[i, 3], pars[i, 4], pars[i, 5], pars[i, 6])
}
ggarrange(plotlist = lapply(res, plot))