作者toto01 (ur ...)
看板Statistics
标题[程式] 使用Monte Carlo估计power
时间Sat Jun 4 22:16:15 2011
小弟使用R去估计power,估计的方式是Monte Carlo,
mu <- 500 # The null hypothesis mean rqual 500
sigma <- 100 # Standard deviation
n <- 20 # The sample size of X is 20
nsim <- 10000 # The number of simulations
pv <- rep(NA, nsim) # empty vector to store p-values from t-test
for (i in 1:nsim){
set.seed(i)
x <- rnorm(n, mu, sigma)
pv[i] <- t.test(x, alternative = "two.sided",mu = 500)$p.value
}
mean(pv < 0.05) # The significat level is alpha = 0.05
hist(pv,prob=T)
最後算出power!
小弟想请教如何对这模拟10000次的t-test画出empirical power curve
PS:模拟10000次的结果用长条图呈现似乎不是常态分布
诚心求教
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 68.81.95.37
1F：→ bmka:如果你是simulate under the null,那麽这10000个p-value值 06/04 22:49
2F：→ bmka:画histogram应该长得像uniform distribution,不会是normal 06/04 22:50
3F：→ bmka:(就算不是simulated under the null,也通常不会是normal) 06/04 22:51
4F：→ toto01:那如何对这模拟10000次的t-test画出empirical power curve 06/04 22:56
5F：→ toto01:? 06/04 22:57
6F：→ toto01:没错!的确长得像uniform distribution 06/04 22:57
7F：→ bmka:你先搞清楚power 是什麽再说 06/05 10:36