Created
November 1, 2022 01:26
-
-
Save chenpanliao/6e4d99580e9f58b710f34adb368df2a0 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 兩個獨立常態母體,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)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment