爲夏皮羅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"))