一个关于R中正态性测试的小型模拟研究

我正在进行一个小型模拟研究来判断两个正态性测试真的有多好。我的计划是生成多个不太多观测值的正态分布样本，并确定每次测试拒绝正常性假设的频率。一个关于R中正态性测试的小型模拟研究

的（不完全的）代码我迄今为止是

library(nortest) 
    y<-replicate(10000,{ 
    x<-rnorm(50) 
    ad.test(x)$p.value 
    ks.test(x,y=pnorm)$p.value 
    } 
    )

现在我想算这些p值是比0.05为每个测试更小的比例。你能告诉我我该怎么做吗？我很抱歉，如果这是一个新手问题，但我自己是新来的R.

谢谢。

来源

2014-12-04 JohnK

library(nortest) 
nsim <- 10000 
nx <- 50 

set.seed(101) 
y <- replicate(nsim,{ 
    x<-rnorm(nx) 
    c(ad=ad.test(x)$p.value, 
     ks=ks.test(x,y=pnorm)$p.value) 
    } 
) 
apply(y<0.05,MARGIN=1,mean) 
##  ad  ks 
## 0.0534 0.0480

使用MARGIN=1告诉apply取平均值跨行，而不是列 - 这是明智的考虑到replicate()的内置简化生产顺序。

对于这种类型的例子，任何标准测试的类型I错误率将接近其标称值（在该示例中为0.05）为极端。

来源

2014-12-04 22:30:06

PS你可以通过一次选取所有Normal值并将它们放入一个矩阵中来加速这一点... – 2014-12-04 22:31:17

谢谢，就是这样。如果我可能会问，为什么你把1作为第二个参数在apply函数中？ – JohnK 2014-12-04 22:37:36

如果分别运行每个测试，那么您可以简单地计算y中存储的哪些值小于0.05。

y<-replicate(1000,{ 
    x<-rnorm(50) 
    ks.test(x,y=pnorm)$p.value}) 
length(which(y<0.05))

来源

2014-12-04 22:26:16 darwin

谢谢，但我希望每个测试都在同一个样本上。 – JohnK 2014-12-04 22:27:37

您的代码不输出p值。你可以这样做：

rep_test <- function(reps=10000) { 

    p_ks <- rep(NA, reps) 
    p_ad <- rep(NA, reps) 

    for (i in 1:reps) { 
    x <- rnorm(50) 
    p_ks[i] <- ks.test(x, y=pnorm)$p.value 
    p_ad[i] <- ad.test(x)$p.value 
    } 

    return(data.frame(cbind(p_ks, p_ad))) 
} 

sum(test$p_ks<.05)/10000 
sum(test$p_ad<.05)/10000

来源

2014-12-04 22:30:40 eamcvey

一个关于R中正态性测试的小型模拟研究

回答

相关问题