2

如果 p 变量正态随机向量的样本具有理论平均值，我希望我可以通过 hotellings 测试证明。但是与 ks.test 的交叉检查，如果来自 HottelingsT2 函数的分布与 HottelingsT2-Test 使用的测试统计量的分布相匹配，则失败。这意味着模拟实验的平均值不是 0，但显然它们有。所以上下文应该有问题。有没有一些错误？

require(mvtnorm)
require(ICSNP)
subject<-50
treatment<-4
V<-matrix(c(644.03100226056, 184.319025225855, 572.5312199559, 143.106678641056, 184.319025225855, 73.5310268006399, 230.838267981476, 130.977532385651, 572.5312199559, 230.838267981476, 736.378779002912, 429.445506266528, 143.106678641056, 130.977532385651, 429.445506266528, 435.124191935888),treatment,treatment)

experiment<-list()
R<-3000
seed<-split(1:(R*subject),1:R)
for(i in 1:R){
  e<-c()
  for(j in 1:subject){
    set.seed(seed[[i]][j]) 
    e<-c(e,rmvnorm(mean=rep(0,treatment),sigma=V,n=1,method="chol"))
   }
   experiment<-c(experiment,list(matrix(e,subject,treatment,byrow=T)))
 }

 p.values<-c()
 for(e in experiment){
   fit<-lm(e~1)
   p.values<-c(p.values,HotellingsT2(e, mu=rep(0,treatment))[["p.value"]])
 }

 ks.test(p.values, punif,alternative = "two.sided")
4

2 回答 2

4

Hong Ooi 关于这是一个问题是正确的set.seed。我在发布您的代码时运行了它并得到以下结果：

> ks.test(p.values, punif,alternative = "two.sided")

    One-sample Kolmogorov-Smirnov test

data:  p.values
D = 0.0615, p-value = 2.729e-10
alternative hypothesis: two-sided

但是如果你改变你的代码，这样：

... everything the same before here ...
experiment <- list()
R <- 3000 # experiment
set.seed(42) # set new seed
for (i in 1:R) { # for each of 3000 experiments
  e <- c() # empty vector
  for (j in 1:subject){ # for each of 50 subjects
    e <- c(e,rmvnorm(mean=rep(0,treatment),sigma=V,n=1,method="chol"))
  }
  experiment <- c(experiment,list(matrix(e,subject,treatment,byrow=T)))
}
... everything the same after here ...

然后你会得到以下信息：

> ks.test(p.values, punif,alternative = "two.sided")

    One-sample Kolmogorov-Smirnov test

data:  p.values
D = 0.0122, p-value = 0.7613
alternative hypothesis: two-sided

本质上，通过在每次迭代中重新设置随机种子，即使您小心选择不同的值，您仍然会消除连续抽签的独立性。

于 2013-06-29T18:18:45.353 回答
2

我没有检查代码，但如果这与 Klaus 的另一篇文章中描述的问题相同，我不会感到惊讶：Using Kolmogorov Smirnov Test in R。基本上，不要放在set.seed循环的中间：设置一次，在代码的顶部，然后不理会它。

于 2013-06-29T18:03:56.897 回答