kritische Wert für Shapiro Wilk Test

Question

Ich versuche, den kritischen W-Wert für einen Shapiro Wilk-Test in R.

 Shapiro-Wilk normality test 

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

mit n = 50 und Alpha zu erhalten = .05, weiß ich, dass die kritische Wert W = .947, indem die Tabelle für kritische Werte durchgeführt wird. Wie bekomme ich diesen kritischen Wert, wenn ich R verwende?

Answer 1

Rechnen kritische Werte direkt ist nicht einfach (siehe CrossValidated answer); Was ich hier habe, ist im Wesentlichen dasselbe wie das, was in dieser Antwort enthalten ist (obwohl ich es selbständig gefunden habe und diese Antwort etwas verbessert, indem ich Ordnungsstatistiken anstelle von Stichproben verwende). Die Idee ist, dass wir eine Probe progressiv mehr nicht-normal machen können, bis sie genau den gewünschten p-Wert erreicht (in diesem Fall 0,05), dann sehen Sie, welche W-Statistik dieser Probe entspricht.

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

Die Antwort (0,9540175) ist nicht ganz die gleiche wie die Antwort, die Sie erhalten, weil R eine Annäherung an den Shapiro-Wilk-Test verwendet. Soweit ich weiß, stammten die tatsächlichen S-W-Tabellen für kritische Werte vollständig von Shapiro und Wilk 1965 Biometrikahttp://www.jstor.org/stable/2333709 p. 605, die nur sagt "Basierend auf angepasster Johnson (1949) S_B Approximation, siehe Shapiro und Wilk 1965a für Details" - und "Shapiro und Wilk 1965a" bezieht sich auf ein unveröffentlichtes Manuskript! (S & W hat im Wesentlichen viele Normalwerte abgetastet, die SW-Statistik berechnet, glatte Näherungen der SW-Statistik über einen Bereich von Werten konstruiert und die kritischen Werte aus dieser Verteilung genommen).

Ich habe auch versucht, dies mit brutaler Gewalt zu tun, aber (siehe unten), wenn wir naiv sein wollen und nicht Kurvenanpassung als SW tat, werden wir viel größere Proben brauchen ...

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

original S & W Werte vergleichen (transkribiert von den Papieren) mit der R Näherung:

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