критическое значение для теста Шапиро Wilk

Question

Я пытаюсь получить критическое W значение для Wilk теста Шапиро в R.

 Shapiro-Wilk normality test 

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

С п = 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 использует аппроксимацию теста Шапиро-Wilk. Насколько я знаю, фактические таблицы критических значений S-W полностью проистекают из Shapiro и Wilk 1965 Biometrikahttp://www.jstor.org/stable/2333709 p. 605, в котором говорится только «Основанный на установленном приближении Джонсона (1949) S_B, см. Шапиро и Уилк 1965a для подробностей» и «Шапиро и Уилк 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)) 
} 

Сравнить оригинальные S & значения W (транскрибируется с бумагами) с аппроксимацией R:

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"))