# What is the CDF of the Shapiro-Wilk $W$ statistic?

I am trying to understand how Shapiro-Wilk test functions. So far, I have come across the following links:

How are the values in the Shapiro Wilk weight table calculated?

which explains how to calculate the coefficients in the shapiro-wilk test. And then:

which explains the routine to calculate the test-statistic and presents a script for calculating the p-value. In my mind, I am comparing the test to a standard t-test where the test statistics has a certain CDF and it helps us to calculate the p-value. I cannot figure this out in the case of Shapiro-Wilk test.

• Finally got around to looking at Royston's papers on P-values of Shapiro-Wilk tests. Please see Addendum to my Answer. – BruceET Aug 24 '18 at 9:02
• BIM welcome to CV. Are you specifically looking for the CDF of $W$? – Alexis Aug 24 '18 at 16:09
• @Alexis This is exactly what I am looking for. I am mostly comparing the test to the classic t-test in my head and I was wondering if Shapiro-Wilk acts the same way or not. – BIM Aug 24 '18 at 16:56
• I recommend editing the title to something like "What is the CDF of the Shapiro-Wilk $W$ statistic?" – Alexis Aug 24 '18 at 17:08

There are many tests of normality. Some of them are based on Q-Q plots (normal probability plots) of the data, measuring according to various criteria how "nearly linear" the Q-Q plot is.

Intuitive view of Shapiro-Wilk test. In visual inspection of a Q-Q plot, the extreme values of a truly normal sample can appear to give an undue departure from a straight line. Here are Q-Q plots of three normal samples of size $n = 500.$ For each random sample from $\mathsf{Norm}(0,1),$ the central portion of the plot seems very close to linear, but the tails look a "little wobbly." (These are not cherry-picked examples, they are the first three normal samples of size 500 resulting from a set.seed statement in R based on today's date.) One feature of the Shapiro-Wilk test is that it tends to 'down-weight' observations in the two tails. The three P-values from shapiro.test are 0.3218225, 0.7221126, and 0.8429852, respectively (all above 0.05); so all three samples are consistent with sampling from a normal population.

set.seed(822)
shapiro.test(rnorm(500))$p.value  0.3218225 shapiro.test(rnorm(500))$p.value
 0.7221126
shapiro.test(rnorm(500))$p.value  0.8429852  You can check on this site and more generally online for discussions of the technical details of computing the Shapiro-Wilk test statistic. But for your purposes, perhaps this intuitive description will be a useful start. Power of Shapiro-Wilk test. As mentioned in the NIST Handbook (and its references) the Shapiro-Wilk test is known for its high power against various alternatives, compared with other tests of normality. In order to make sense of power, you must have a specific significance level and non-normal distribution in mind. Student's t distribution with$\nu = 10$degrees of freedom (symmetrical with heavy tails,$\mu = 0,\, \sigma\approx 1.12)$and the distribution$\mathsf{Gamma}(\text{shape}=10, \text{scale}=10)$(right-skewed,$\mu=1,\, \sigma\approx 0.316)$are somewhat close to normal in shape, as illustrated in the figure below. My impression is that the power in specific situations is usually obtained by simulation. The following simulations in R show approximate powers of the Shapiro-Wilk normality test at level 5% against each of these two alternative distributions for samples of size$n=500.$The power is the probability of rejection when data are from the alternative distribution. Respective power values are about 65% for$\mathsf{T}(10)$and above 99% for$\mathsf{Beta}(10,10).$[With$m = 100,000$iterations one can expect about two-place accuracy.] set.seed(818); m = 10^5; n = 500 p.val = replicate(m, shapiro.test(rt(n, 10))$p.val)
mean(p.val < 0.05)
 0.65474

set.seed(2018);  m = 10^5;  n = 500
p.val = replicate(m, shapiro.test(rgamma(n, 10, 10))$p.val) mean(p.val < 0.05)  0.99978  Addendum: P-values. Traditionally, significance of a Shapiro-Wilk test was determined using tabled values of the test statistic$W.$The P-values used in modern software programs seem to be mainly due to the work of Patrick Royston published in Applied Statistics: (1982) Vol. 31, 115-124 and 176-180, and (1995) Vol. 44, 547-551. [See References in the R documentation for shapiro.test.] The second 1982 paper extended the algorithm for P-values to accommodate$n \le 2000.$The general method is to find a serviceable transformation to make the test statistic$W$approximately normal. Currently, R accepts datasets with sizes$3 \le n \le 5000.$However, the best accuracy is not guaranteed for P-values above$0.1.$Simulation can provide an intuitive idea of the null distribution of$W$for a particular$n.$In R, the code set.seed(1234); rnorm(500) produces a standard normal sample of size$n = 500$for which R gives$W = 0.99623$with p-value$= 0.2848.$To simulate the distribution of$W$for$n = 500,$we use the program below: set.seed(1234); x = rnorm(500); w.obs=shapiro.test(x)$stat
set.seed(2018);  m = 10^5;  n = 500
w = replicate(m, shapiro.test(rnorm(n))$stat) mean(w < w.obs)  0.28396  The histogram below shows the simulated distribution of$W$along with the observed value of$W = 0.9962$for for our sample with P-value about 0.2848. The null distribution of exact P-values from a test with a continuous test statistic is$\mathsf{Unif}(0,1).$If we run a simulation similar to the one above, but capturing P-values (instead of tests statistics), we can see how close the Shapiro-Wilk P-values come to this uniform distribution. Because our P-values are not based on a continuous test statistic and are not exactly correct, we do not quite have a perfect fit to uniform. [The left-most bar for P-values$< 0.05$(rejection) has area approximately$0.05.\$] tl;dr: There is no known distribution. i.e. it doesn't have a name. You use simulations instead.

I found this to give some insight:

Basically, Shapiro and Wilk calculated the distribution of their statistic only for $$n=3$$, to be a truncated $$Beta(\frac{1}{2}, \frac{1}{2})$$ distribution for $$\frac{3}{4} \le w \le 1$$ and zero elsewhere.

For $$3 < n \le 20$$, the coefficients values in the statistics (expected values of ordered statistic, times inverse of co-variance matrix, normed) were calculated precisely, but there's no name for the distribution. For $$20 < n \le 50$$ these coefficients were only approximated, and (monte carlo) simulations were used.

Some researchers after tried to find the asymptotic distribution, which in itself has no name.

Leslie, Stephens, and Fotopoulos showed that $$n(w - \mathbb{E}(w)) \sim -\sum_{k=3}^{\infty}\frac{Z_k^2-1}{k}$$ where $$Z_k \sim N(0,1)$$ and were i.i.d.

But in any case that too was shown to converge

"at a painfully slow rate. ... For such reasons, the limit distribution seems not to be of much practical use, and for $$n > 50$$, Monte Carlo simulation seems still to be needed."