为夏皮罗Wilk检验临界值

Question

我试图得到一个夏皮罗Wilk检验临界W值在R.

 Shapiro-Wilk normality test 

data: samplematrix[, 1] 
W = 0.69661, p-value = 7.198e-09 

其中n = 50和α= 0.05，我知道临界值W = .947，通过执行临界值表。但是，如何使用R来获得这个临界值？

Answer 1

直接计算临界值并不容易（请参见CrossValidated answer）;我在这里得到的结果基本上与答案中的相同（尽管我独立提出了它，并且通过使用顺序统计量而不是随机样本稍微改进了该答案）。我们的想法是，我们可以使样本逐渐变得更加非正态，直到获得完全所需的p值（在这种情况下为0.05），然后查看该样本对应的W统计量。

## compute S-W for a given Gamma shape parameter and sample size 
tmpf <- function(gshape=20,n=50) { 
    shapiro.test(qgamma((1:n)/(n+1),scale=1,shape=gshape)) 
} 
## find shape parameter that corresponds to a particular p-value 
find.shape <- function(n,alpha) { 
    uniroot(function(x) tmpf(x,n)$p.value-alpha, 
      interval=c(0.01,100))$root 
} 
find.W <- function(n,alpha) { 
    s <- find.shape(n,alpha) 
    tmpf(s,n=n)$statistic 
} 
find.W(50,0.05) 

答案（0.9540175）是不太一样的，你得到的答案，因为R使用的近似值夏皮罗 - 威尔克测试。据我所知，实际的S-W临界值表完全来自Shapiro和Wilk 1965 Biometrikahttp://www.jstor.org/stable/2333709 p。 605，它只说“基于拟合约翰逊（1949）S_B近似，详见Shapiro and Wilk 1965a” - 和“Shapiro and Wilk 1965a”是指未发表的手稿！ （S & W基本上取样了许多正常偏差，计算了SW统计量，在一系列值上构造了SW统计量的平滑近似值，并从该分布中取出了临界值）。

我也试图通过强力做到这一点，但（见下文），如果我们想成为幼稚，而不是做曲线拟合为SW一样，我们需要更大的样本...

find.W.stoch <- function(n=50,alpha=0.05,N=200000,.progress="none") { 
    d <- plyr::raply(N,.Call(stats:::C_SWilk,sort(rnorm(n))), 
        .progress=.progress) 
    return(quantile(d[1,],1-alpha)) 
} 

的R近似比较原件S &的W值（从文件转录）：

SW1965 <- c(0.767,0.748,0.762,0.788,0.803,0.818,0.829,0.842, 
    0.850,0.859,0.866,0.874,0.881,0.887,0.892,0.897,0.901,0.905, 
    0.908,0.911,0.914,0.916,0.918,0.920,0.923,0.924,0.926,0.927, 
    0.929,0.930,0.931,0.933,0.934,0.935,0.936,0.938,0.939,0.940, 
    0.941,0.942,0.943,0.944,0.945,0.945,0.946,0.947,0.947,0.947) 
    Rapprox <- sapply(3:50,find.W,alpha=0.05) 
    Rapprox.stoch <- sapply(3:50,find.W.stoch,alpha=0.05,.progress="text") 
    par(bty="l",las=1) 
    matplot(3:50,cbind(SW1965,Rapprox,Rapprox.stoch),col=c(1,2,4), 
      type="l", 
      xlab="n",ylab=~W[crit]) 
    legend("bottomright",col=c(1,2,4),lty=1:3, 
     c("SW orig","R approx","stoch"))