Repository: JuliaStats/HypothesisTests.jl
Pull request #124 (2018-01-04 18:51:59)

Shapiro-Wilk normality test
implements ShapiroWilkTest following
> PATRICK ROYSTON
> Approximating the Shapiro-Wilk W-test for non-normality
> *Statistics and Computing* (1992) **2**, 117-119
> DOI: [10.1007/BF01891203](https://doi.org/10.1007/BF01891203)


This is work in progress, all comments are welcome! 
I tried to follow the paper closely, but **copied**  e.g. computed constants from the original [swilk.f](https://github.com/scipy/scipy/blob/master/scipy/stats/statlib/swilk.f). 
Please let me know if You are ok (license-wise) with this.

These polynomial interpolations are computed at loadtime for  speed.
Currently
```
julia> using BenchmarkTools
julia> srand(1);
julia> k = 5000;
julia> X = sort(randn(k));
julia> @btime ShapiroWilkTest(X);
  112.523 μs (18 allocations: 137.45 KiB)
julia> swc = SWCoeffs(k);
julia> @btime ShapiroWilkTest(X, swc);
  62.476 μs (11 allocations: 78.50 KiB)
``` 
whereas calling `swilkfort` directly
```
julia> srand(1);
julia> X = sort(randn(5000));
julia> A = zeros(2500);
julia> @btime swilkfort!(X, A);
  205.763 μs (10 allocations: 224 bytes)
julia> @btime swilkfort!(X, A, false);
  75.692 μs (10 allocations: 224 bytes)
```

Still missing:

- [x] documentation
- [ ] censored data
- [ ] polynomial interpolation with better precision (?),

and probably more.

I tried to compute exact values of `SWCoeffs` (via MonteCarlo simulation), but the results I'm getting are off the reported ones in Table 1 *op. cit.* Woule be glad if anyone could help.

```julia
> function approximate_A(N, samps)
    s = [sort(randn(N)) for i in 1:samps]
    m = sum(s[i][1:N] for i in 1:samps)/samps
    ss = vcat(s...)
    Vinv = inv(cov(ss, ss))
    A = (-m'*Vinv/sqrt(m'*Vinv*Vinv*m))'
    return ((A - reverse(A))/2)[1:div(N,2)]
end;
```
`approximate_A(10,1000000)` results in
```
[0.5469493640066541, 0.3559315327467397, 0.23322491989287789, 0.1336531031109471, 0.043609401791162974]
```
Compare `swilk.f`'s, and `SWCoeffs(10).A`:
```
[0.5737146934337275, 0.3289700603561813, 0.21434902886647542, 0.12279063037794058, 0.04008871216291867]
[0.5737147066903874, 0.3289700464878161, 0.21434901803439887, 0.1227906248657577,  0.04008871105102476]
```
## Rahman & Govindarajulu
A different implementation of SWCoeffs following 
> M. Mahibbur Rahman & Z. Govindarajulu,  
> A modification of the test of Shapiro and Wilk for normality  
> *Journal of Applied Statistics* 24:2, 219-236, 1997  
> DOI: [10.1080/02664769723828](http://dx.doi.org/10.1080/02664769723828)

```julia
function HypothesisTests.SWCoeffs(N::Int, ::Type{Val{:Rahman}})
    M = [norminvcdf(i/(N+1)) for i in 1:N];
    normM = normpdf.(M)

    c = (N+1)*(N+2)
    d = 2c.*normM.^2
    dl= [-c*normM[i]*normM[i+1] for i in 1:N-1]
    Vinv = SymTridiagonal(d, dl)
    C = sqrt(M'*(Vinv*Vinv)*M)

    Astar = -(M'*Vinv)'/C
    return SWCoeffs(N, Astar[1:div(N,2)])
end
```
But they don't provide critical values for their statistic (or I couldn't find it), so that would require further simulation. Also their `SWCoeffs(10, Val{:Rahman})` are quite different:
```
[0.6010309510135708, 0.2975089851946684, 0.19212508233094944, 0.10979215067183574, 0.03583065785524651]
```
On the other hand `ALGORITHM AS R94` is a **standard**.
  
  

kalmarek: >The main issue is that the swilk.f is written for Float32 and not Float64 so you'll probably need to compute new polynomial coefficient for the Float64 implementation.

Yep, I tried to do this, but failed obtaining the "true" `SWCoeffs` by semi-exact means

Even obtaining the **"exact"** expected values of order statistics `E(r:n)` (our vector `m`) seems to be a problem. Convention: we obtain the first half of `-order statistics N(0,1)`.

A simple MC produces highly variable results:
```julia
function mMC(N::Int,samps=10^5)
    X = Array{Float64}(N)
    s = zeros(div(N,2))
    for i in 1:samps
        for j in 1:N
            X[j] = randn()
        end
        sort!(X)
        for j in 1:div(N,2)
            s[j] += X[j]
        end
    end
    return -s/samps
end
d = mMC(101, 10^7) - mMC(101, 10^7);
norm(d,1)
0.0019123699000122493
```
A more or less exact order statistics can be obtained by integration directly, eg. following 
> `Algorithm AS 177`, 
> Royston, J.P. 
> [Expected Normal Order Statistics (Exact and Approximate)](http://www.jstor.org/stable/2347982)
But I don't know how exact are normcdfs, etc., quadgk, i.e. how exact is the following method?

```julia
function Royston(n,r, logs=[log(i) for i in 1:n], o=7*sqrt(n))
    C = sum(logs) - sum(view(logs, 2:r-1)) - sum(view(logs, 2:n-r))
    I(x) = C + (r-1)*log(normccdf(x)) + (n-r)*log(normcdf(x)) - 1/2*(log(2π) + x^2)
    return quadgk(x->x*exp(I(x)), -9.0-log(n),9.0+log(n), 
            order=floor(Int, o), 
            maxevals=floor(Int, log(n)*10^7))
end
function mRoyston(n, o=7*sqrt(n))
    logs=[log(i) for i in 1:n]
    m = [Royston(n,r, logs, o) for r in 1:div(n,2)]
    err = norm([i[2] for i in m],1)
    @show err
    return [i[1] for i in m]
end
```
Royston's `Algoritm AS R94` uses Blom approximation:
```julia
Blom(n,r) = -norminvcdf((r - 3/8)/(n + 1/4))
mBlom(n) = [Blom(n,r) for r in 1:div(n,2)]
```
and Rahman (cited above) used a simple
```julia
Rahman(n,r) = -norminvcdf(r/(n + 1))
mRahman(n) = [Rahman(n,r) for r in 1:div(n,2)]
```
# Comparing
I believe the direct integration to be the most accurate, so I'd lean towards regarding it as "true". Comparing other methods against it gives
```julia
const N = 5000
m_Ro = mRoyston(N);
m_mc = mMC(N, 10^6);
m_Bl = mBlom(N);
m_Rh = mRahman(N);

p = plot(yaxis=:log10, title="Relative to Royston",size=(800,600))
for (X, l) in [(m_mc, "MonteCarlo"), 
        (m_Bl, "Blom approximation"), 
        (m_Rh, "Rahman approximation"), 
   ]
    plot!(abs.((X-m_Ro)./m_Ro), label=l)
end
p
```
![newplot 2](https://user-images.githubusercontent.com/1666118/34645436-aa939930-f34d-11e7-83d2-60568bb80801.png) 

Question:
 1. Is storing/caching values of `m` (and later of `SWCoeffs`) an option (for sizes `N in 4:1:5000)
 2. Should we follow `AS R94` for the ShapiroWilk test, or should we fit our own curves for distributions of `W`?
andreasnoack: > But I don't know how exact are `normcdf`s, etc., `quadgk`, i.e. how exact is the following method?

Both should be high quality and `quadgk` gives can provide an estimate of the error. Notice also that there are `normlogcdf` which might give you a little extra accuracy over the explicit `log`.

So all-in-all I think this should be fine for computing the mean order statistics. Have you tried to compute the covariances yet? They are needed for fitting the polynomials, right?

> 1. Is storing/caching values of m (and later of SWCoeffs) an option (for sizes `N in 4:1:5000)

I guess we could do that. We could also try to do this as part of a generated function that would compute everything on demand but it might make the first call pretty slow. It would be pretty cool on the other hand.

> Should we follow AS R94 for the ShapiroWilk test, or should we fit our own curves for distributions of W?

Given how inaccuratly they computed the order statistics back then, I guess it would be good if we could recompute them.

Another thing, do you have experience with Weisberg-Bingham version? It is so much easier to compute so do you know why people seem to prefer SW? Path dependency? I tried to look for power comparisons I couldn't find anything that suggested it to be significantly worse.
kalmarek: For covariances I plan to do the following: [DavisStephens1977] goves an algorithm to compute covariance matrix of normal order statistics given
  1. $V_11$ - variance of the largers order statistics → Approximated by [SheaScallon1988]
  2. $EX_1$, $EX_2$ - the expected value of the largest and second largest order statistics → ~exact by [Royston1982]
  3. $\Sigma^2$ - sum of squared of order statistics → Approximated (with bounds on the error) by [Balakrishnan1984]

- [DavidStephens1977]: C. S. Davis & M. A. Stephens The covariance matrix of normal order statistics *Communications in Statistics - Simulation and Computation* Vol. **6** , Iss. 1, 1977 [DOI:10.1080/03610917708812028](https://doi.org/10.1080/03610917708812028)  
- [SheaScallon1988]: B. L. Shea and A. J. Scallon Remark AS R72: A Remark on Algorithm AS 128. Approximating the Covariance Matrix of Normal Order Statistics  *Journal of the Royal Statistical Society. Series C (Applied Statistics)* Vol. **37**, No. 1 (1988), pp. 151-155 [DOI:10.2307/2347514](https://doi.org/10.2307/2347514)  
- [Royston1982]:  J. P. Royston Algorithm AS 177: Expected Normal Order Statistics (Exact and Approximate) *Journal of the Royal Statistical Society. Series C (Applied Statistics)* Vol. **31**, No. 2 (1982), pp. 161-165 [DOI:10.2307/2347982](https://doi.org/10.2307/2347982)  
- [Balakrishnan1984]: N. Balakrishnan Algorithm AS 200: Approximating the Sum of Squares of Normal Scores *Journal of the Royal Statistical Society. Series C (Applied Statistics)* Vol. **33**, No. 2 (1984), pp. 242-245 [DOI:10.2307/2347458](https://doi.org/10.2307/2347458)

But I may not have enough time to this during this week.

> We could also try to do this as part of a generated function that would compute everything on demand but it might make the first call pretty slow. It would be pretty cool on the other hand.

We may have a look at speeds and decide when I finish the plan above; but I guess cached (or `@memoized`?) version of `SWCoeff(N::Int)` will be the way to go.

> Another thing, do you have experience with Weisberg-Bingham version? It is so much easier to compute so do you know why people seem to prefer SW? Path dependency? I tried to look for power comparisons I couldn't find anything that suggested it to be significantly worse.

Roystons `AS R94` is essentially Weisberg-Bingham + corrections + fitting for the distribution of `W`-statistic. `AS R94` uses [Bloms expectation of normal order statistics](https://github.com/JuliaStats/HypothesisTests.jl/pull/124/commits/8774c6f7c186f1588b61b610df1336e05099c572#diff-afb30590cfaf785037d4c7a73b375177R70) `m` and then corrects just the [first two coefficients](https://github.com/JuliaStats/HypothesisTests.jl/pull/124/commits/8774c6f7c186f1588b61b610df1336e05099c572#diff-afb30590cfaf785037d4c7a73b375177R85)
This is an approximation of (the perfect) Shapiro-Wilk, if You like. If we can (on modern hardware) do it the "exact" way (and make it fast) maybe we should?
 
kalmarek: > Another thing, do you have experience with Weisberg-Bingham version?

As a matter of fact, I'm not a statistician, doing pure math, coding in julia "after hours" ;-). But for most of my "empirical science" friends SW was the one to go to, and seems to have the best power for small samples.
kalmarek: Ok, variances/covariances turned out harder than I thought:
I implemented the exact formulas from [Godwin Some Low Moments of Order Statistics](https://projecteuclid.org/euclid.aoms/1177730036) in [normordstats.jl](https://gist.github.com/kalmarek/67b6459589823dba0cacbb1aec17be8d), please have a look. These become numerically unstable for `n ~ 25`, but for `n=20` agree to `~10` digits with [Table 1 of Parrish](https://doi.org/10.1080/03610919208813009). He used `128` bits of precision to compute tables of products of order statistics up to `n=50` (and I guess he didn't progress any further because of loss of precision). So maybe my cheap implementation isn't that bad ;-) 

(BTW. Do You have, by a chance, access to a paper copy of Parrishs article (to scan and send me)? The digital one I have is without text layer and of such poor quality that OCRing those tables produces errors in every row. Which for an article to BORROW for 24h for $50 is **OUTRAGEOUS**. )

With these I can tell that Roystons approximation of `SWCoeffs` is surprisingly good, whereas Rahmans is (surprisingly) bad? Maybe there is an error in my implementation of the latter:
```julia
blom(n,i) = norminvcdf((i - 3/8)/(n + 1/4))
rahman(n,i) = norminvcdf(i/(n + 1))

struct Blom
    n::Int
end

struct Rahman
    n::Int
end

OrderStatistics.expectation(B::Blom) = [blom(B.n,i) for i in 1:B.n]
OrderStatistics.expectation(R::Rahman) = [rahman(R.n,i) for i in 1:R.n]

# M. Mahibbur Rahman & Z. Govindarajulu,  
# A modification of the test of Shapiro and Wilk for normality  
# *Journal of Applied Statistics* 24:2, 219-236, 1997  
# DOI: [10.1080/02664769723828](http://dx.doi.org/10.1080/02664769723828)

function invV(N, M)
    pdfM = normpdf.(M)
    d = pdfM.^2
    dl= [-pdfM[i]*pdfM[i+1] for i in 1:N-1]
    
    return (N+1)*(N+2)*SymTridiagonal(d, dl)
end

function HypothesisTests.SWCoeffs(R::Rahman)    
    M = expectation(R)
    A = invV(R.n, M)
    C = sqrt(M'*(A*A)*M)

    Astar = (M'*A)'/C
    return SWCoeffs(R.n, Astar[1:div(R.n,2)])
end

function HypothesisTests.SWCoeffs(OS::NormOrderStatistic)
    m = expectation(OS)
    iV = inv(cov(OS))
    A = m'*iV./sqrt(m'*iV*iV*m)
    return SWCoeffs(OS.n, A[1,1:div(OS.n,2)])
end
```
Then:
```julia
N = 25
@time ExactA = SWCoeffs(NormOrderStatistic(N));
plot(SWCoeffs(N).A, label="Blom + 0-correlation OrdStats")
plot!(SWCoeffs(Rahman(N)).A, label="Rahman + limit inverse correlation")
plot!(-ExactA.A, label="Exact formulas, BigFloat")
```
![newplot 4](https://user-images.githubusercontent.com/1666118/35403570-67887108-0200-11e8-9fd6-d4417d935e85.png)
(You can see instability kicking in)

Right now I can not think of anything better than the mentioned above approximation by Stephen&David (using Taylor expansion of the expectation of product of order statistics).

Using BigFloats, ArbFloats, etc will not help much, as the precision of e.g. `erf`, hence `normcdf` is [limited to `Float64`](https://github.com/JuliaMath/SpecialFunctions.jl/issues/75), at least for now.
kalmarek: It took more time than I had predicted, but turned out a fun project ;-)
I'm looking for a piece of advise what to do with the code, based on the results below. At the moment I think that numerically Shapiro-Wilk Test should be left "as is", until someone could point what did Royston fit his polynomials against. We should of course use the exact values of Parrish for `N≤20` (or even `N≤50` if I will be able to fix the loss of precision mentioned below).

# Normal order statistics
For implementation see https://github.com/kalmarek/HypothesisTests.jl/tree/normordstats
Using this code You can confirm the covariances computed by Parrish (the 25dp values) for `N≤20`, but beyond that point there seems to be precision loss, which I was not able to trace. 
# David-Johnson formula: 
```julia
function diffG(x)
    s = normpdf(x)
    s² = s^2
    s⁴ = s²*s²
    x² = x^2
    x⁴ = x²*x²
    
    return (
        1/s, 
        x/s², 
        (1 + 2x²)/(s*s²), 
        x*(7 + 6x²)/s⁴, 
        (7 + 46x² + 24*x⁴)/(s*s⁴),
        x*(127 + 326x² + 120x⁴)/(s²*s⁴)
        )
end
function covDavidJohnson(i::Int,j::Int,n::Int, order=3)
    if j < i
        return covDavidJohnson(j,i,n)
    end
        
    pᵢ, pⱼ = i/(n+1), j/(n+1)
    qᵢ, qⱼ = 1 - pᵢ, 1 - pⱼ
    Gᵢ = diffG(norminvcdf(pᵢ))
    Gⱼ = diffG(norminvcdf(pⱼ))
    X = 1/(n+2)
    T = Vector{eltype(Gᵢ)}()
    if order >= 1
        T1 = X  *pᵢ*qⱼ*Gᵢ[1]*Gⱼ[1]
        push!(T, T1)
    end
    if order >= 2
        T2 = X^2*pᵢ*qⱼ*(
                (qᵢ-pᵢ)*Gᵢ[2]*Gⱼ[1] + 
                (qⱼ-pⱼ)*Gᵢ[1]*Gⱼ[2] + 
                0.5(
                    pᵢ*qᵢ*Gᵢ[3]*Gⱼ[1] +
                    pⱼ*qⱼ*Gᵢ[1]*Gⱼ[3] +
                    pᵢ*pⱼ*Gᵢ[2]*Gⱼ[2] 
                    )
                )
        push!(T, T2)
    end
    if order >= 3
        T3 = X^3*pᵢ*qⱼ*(
                -(qᵢ-pᵢ)*Gᵢ[2]*Gⱼ[1] - (qⱼ-pⱼ)*Gᵢ[1]*Gⱼ[2] + 
                ((qᵢ-pᵢ)^2 - pᵢ*qᵢ)*Gᵢ[3]*Gⱼ[1] + 
                ((qⱼ-pⱼ)^2 - pⱼ*qⱼ)*Gᵢ[1]*Gⱼ[3] +
                1/2*(3(qᵢ-pᵢ)*(qⱼ-pⱼ) + pⱼ*qᵢ - 4pᵢ*qⱼ)*Gᵢ[2]*Gⱼ[2] + 
                5/6*pᵢ*qᵢ*(qᵢ-pᵢ)*Gᵢ[4]*Gⱼ[1] + 
                5/6*pⱼ*qⱼ*(qⱼ-pⱼ)*Gᵢ[1]*Gⱼ[4] + 
                (pᵢ*qⱼ*(qᵢ-pᵢ) + 0.5pᵢ*qᵢ*(qⱼ-pⱼ))*Gᵢ[3]*Gⱼ[2] + 
                (pᵢ*qⱼ*(qⱼ-pⱼ) + 0.5pⱼ*qⱼ*(qᵢ-pᵢ))*Gᵢ[2]*Gⱼ[3] + 
                1/8*(pᵢ^2*qᵢ^2*Gᵢ[5]*Gⱼ[1] + pⱼ^2*qⱼ^2*Gᵢ[1]*Gⱼ[5]) + 
                1/4*(pᵢ^2*qᵢ*qⱼ*Gᵢ[4]*Gⱼ[2] + pᵢ*pⱼ*qⱼ^2*Gᵢ[2]*Gⱼ[4]) + 
                1/12*(2pᵢ^2*qⱼ^2 + 3pᵢ*pⱼ*qᵢ*qⱼ)*Gᵢ[3]*Gⱼ[3]
                )
        push!(T, T3)
    end
    return sum(T)
end
```
Comparing to the high-precision (25d) values obtained by Parrish (`N = 20; log10.(abs.(differences))`):
![newplot 5](https://user-images.githubusercontent.com/1666118/36852148-4c4d1d72-1d6b-11e8-8624-a1b5f4451746.png)
# DavisStephens _Covariance Matrix of Normal order Statistics_
This method uses the David-Johnson formula + relations between correlations, etc:
```julia
function normalizerel!(v::AbstractVector, relative::Union{Vector{Int}, UnitRange})
    Sr = sum(v[relative])
    c = (1-Sr)/(sum(v) - Sr)
    for i in eachindex(v)
        if !(i in relative)
            v[i] *= c
        end
    end
    return v
end
function fillcov!(V, j::Int)
    # copy reversed j-th column to (end-j)-th column
    @views V[j:end-j+1, end-j+1] .= reverse(V[j:end-j+1, j])
    # copy j-th column to j-th row
    @views V[j, j+1:end-j] .= V[j+1:end-j, j]
    # copy (end-j)-th column to (end-j)-th row 
    @views V[end-j+1, j+1:end-j] .= V[j+1:end-j, end-j+1] 
    return V
end
function covDavisStephens(OS::NormOrderStatistic{T}, V11::T=var(OS, 1)) where T
    V = zeros(T, OS.n, OS.n)
    
    V[1,1] = V11
    V[2,1] = moment(OS, 1, 2) - expectation(OS,1)*expectation(OS,2) - T(1)
    
    for i in 3:OS.n
        V[i, 1] = covDavidJohnson(i, 1, OS.n)
    end
    normalizerel!(view(V, :, 1), [1,2])
    fillcov!(V, 1)
    j=1

    for j in 2:div(OS.n, 2)
        V[j,j] = var(OS, j)
        for i in j+1:OS.n-j+1
            V[i,j] = covDavidJohnson(i, j, OS.n)
        end
        normalizerel!(view(V, :, j), [collect(1:j); collect(OS.n-j+1:OS.n)])
        fillcov!(V, j)
    end
    
    if isodd(OS.n) # the middle point
        j = div(OS.n, 2)+1 
        V[j, j] = var(OS, j)
        normalizerel!(view(V, :, j), [j])
        fillcov!(V, j)
    end
    return Symmetric(V, :U)
end

covDavisStephens(N::Int) = covDavisStephens(NormOrderStatistic(N))
```
Again comparing against high-precision values of Parrish (`N = 20; log10.(abs.(differences))`):
![newplot 6](https://user-images.githubusercontent.com/1666118/36852589-525621cc-1d6c-11e8-8c69-549356794e47.png)

The formula is pretty quick:
```julia
N = 2000
@time M = covDavisStephens(N)
@time invM = inv(M);
heatmap(abs.(invM./maximum(invM)), size=(800,700), yflip=true)
# 2.853262 seconds (12.26 M allocations: 365.602 MiB, 2.97% gc time)
# 0.552725 seconds (22 allocations: 62.043 MiB, 1.02% gc time)
```
![newplot 7](https://user-images.githubusercontent.com/1666118/36852774-f4912572-1d6c-11e8-914a-5ede84f6bb78.png)

# Shapiro-Wilk coefficients
Finally the coefficients:
```julia
swcoeff(invcov, m) = ((m'*invcov)'/sqrt(m'*(invcov*invcov)*m))[1:div(length(m),2)]
```

Relative to the exact (=using high-precision covariance by Parrish) value of SW coefficients:
```julia
N=20
m = Float64.(expectation(NormOrderStatistic{BigFloat}(N)))
exact = -swcoeff(inv(Float64.(covd[N])), m)
p = plot(yaxis=:log10)
# plot!(exact, label="Exact (Parrish)")
plot!(abs.(SWCoeffs(N).A)./exact, label="Royston")
plot!(abs.(swcoeff(eye(N), expectation(Blom(N))))./exact, label="Weisberg-Bingham")
plot!(abs.(swcoeff(inv(covDavisStephens(N)), m))./exact, label="DavisStephens")
```
![newplot 8](https://user-images.githubusercontent.com/1666118/36855527-01a1ec90-1d74-11e8-91fd-b70b098aa822.png)

I still have no idea what Royston fitted his coefficients against, but that must have been pretty accurate!

for `N=2000`, the third coefficient of Davis-Stephens is still larger than it should by ~25% (relative to Royston):
![newplot 9](https://user-images.githubusercontent.com/1666118/36858072-6ad957c4-1d7a-11e8-9d57-0f2ed1caa95d.png)

andreasnoack: @kalmarek Sorry for dropping the ball here. What is the current status?
kalmarek: I also had too much to do, couldn't finish this cleanly;)

long story short we either:
0. use Royston approximation (already in place), forgetting that he fit using 32-bits of precision (as everybody do)
1. use Parrish exact values (`N≤25`) of covariances to compute exact SW-coeffs and Royston otherwise (as everybody should do)
2. try to get the order statistics going (and pre-compute the SW-coeffs, e.g. for `N ≤ 50`, or `100`), and use Royston method above;
  - I tried my own solution to compute exact covariances of normal order statistics here https://github.com/kalmarek/HypothesisTests.jl/tree/normordstats, but as noted above, there is some problem with precision loss which I was not able to track down;
  - use Nemo's high quality acb_integrate https://github.com/Nemocas/Nemo.jl/blob/master/src/arb/acb_calc.jl (was the next thing on my list, but never happened)

Happy to finish this, starting the next week, actually (conferences, etc.)
Ph0non: This dedicated work on Shapiro-Wilk-Test is exceptional fantastic 🥇 

Are there any news to complete the work on the SWT?
kalmarek: @Ph0non this effort (it was fun!) grew our from my frustration on the teaching I was forced to do.
and yes, that included teaching statistics ;-)
that semester is over, so I don't have enough frustration in my store right to forge it into code ;-)

but if You are willing to provide documentation, I may look into finishing this, but only in December (conferences abound).
kalmarek: There are three things to be finished here:
1. fix problems with accuracy
2. produce test suite
3. write documentation

I volunteer to finish 1. 
If I succeed then 2 should be rather easy, but not zero-work task, In the documentation we should carefully explain why do we differ from standard packages, and why everybody should be using exact order statistics in SW. 
Ph0non: Some small fixes https://github.com/kalmarek/HypothesisTests.jl/pull/1
kalmarek: @albz could You point me to a python implementation? the one from scipy https://github.com/scipy/scipy/blob/v1.1.0/scipy/stats/morestats.py#L1304 uses standard fortran code, the approximating polynomials were not fit in double precision there. 

right now copyright is not a problem, since I wrote everything from scratch following the paper only. 
kalmarek: I see that people are generally interested so here is some details. In 1992 paper Royston 
1. obtained (by 10000-MC samples) what he calls "exact" coefficients (`A`),
2. took Weisberg-Bingham approximation (`C`) and 
3. fit two quintic polynomials (in `n^{1/2}`) to differences `A_n - C_n` and  `A_{n-1} - C_{n-1}`. Coefficients of those polynomials are visible in the fortran code as `C1` and `C2`;
4. His coefficients for SW consist of corrected `C_n` and `C_{n-1}` and re-normalized rest of `C`. 

I'm happy with both options here:
* providing a pre-computed (but this time exact!) coefficients `A` for sizes up to e.g. 200 and then
* compare the exact coefficients (again) with Weisberg-Bingham and use the Royston approximation method above

Note that my experiments with MC show that it can not provide better precision than `1e-5` or `1e-6`;
Then additional fitting must be provided for the distribution of W-statistic.
Royston claims it's close to ~lognormal for two different set of parameters which needs to be refitted again (these are C3-C7 in the fortran code).

albz: You are right, the python-fortran code is in single precision. That is of little help.
andreasnoack: @kalmarek It's really a shame that this never was completed given all the work you put into this. So hereby my encouragement to resume the work.
kalmarek: @andreasnoack Thanks for the encouragement ;) However unlikely it is, SW test is still in the back of my head.
Hopefully we will be soon bringing [Arblib.jl](https://github.com/kalmarek/Arblib.jl) to a usable state and then I could use its certified integration to re-compute `SW` coefficients.
The plan is just as it was before:
* use exact SW coeffs for small number of samples, say `N≤200` (using `Arblib`)
* use Royston approximation afterwards (this is implemented) (but SW test becomes pretty useless for large sample sizes :)
* re-fit polynomials (logpolynomials?) to distribution(s) of `W`-statistics and find critical thresholds (TODO) 


kalmarek: @andreasnoack in https://github.com/kalmarek/ShapiroWilk/pull/1 `SWCoeffs` are implemented using Arblib;
one can do (an embarrassingly parallel) computation e.g. for `(n=50, prec=128)`, which gives in a couple of minutes ~ 7 digits of precision; one gets more than Float64 has for `prec=256`.

To compute critical/p-values I'd need to
1. approximate a normalizing transform `(1-Wₙ)^t ~ N(µ, σ)` (`µ` could be computed analytically, I'm not sure about `σ`) via MC
2. fit a (exponential? biexponential?, inverse log?) function to `t = t(n)` and use it for all sample sizes
3. use z-scores for the transformed `W`-statistic

I used of the shelf Optim.jl methods to do so:

```julia
using Statistics
using Random

using Distributions
using LinearAlgebra
using Optim


using Revise
using ShapiroWilk


function sample_Wstatistic(sw::ShapiroWilk.SWCoeffs{T}, sample_size) where T
    Wstats = zeros(T, sample_size)
    tmp = similar(Wstats, length(sw))
    for i in 1:sample_size
        randn!(tmp)
        sort!(tmp)
        Wstats[i] = ShapiroWilk.Wstatistic(tmp, sw)
    end
    Wstats
end;

function normal_power_fit(
    t,
    sample,
    qs = 0.005:0.005:0.995,
    normal_qs=quantile.(Normal(0,1), qs),
)

    sample_t = similar(sample)

    for idx in 1:length(sample_t)
        sample_t[idx] = (one(eltype(sample)) - sample[idx])^t
    end

    µ, σ = mean(sample_t), std(sample_t)
    for idx in 1:length(sample_t)
        sample_t[idx] = (sample_t[idx] - μ)/σ
    end

    #@assert issorted(sample_t)

    sample_t_qs = quantile(sample_t, qs, sorted=true)

    for idx in 1:length(sample_t_qs)
        sample_t_qs[idx] = sample_t_qs[idx] - normal_qs[idx]
    end

    return norm(sample_t_qs, 2)
end
```
optionally e.g.
```julia
@time let n = 30, prec=128
    OS = ShapiroWilk.OrderStatisticsArblib.NormOrderStatistic(n, prec=prec)
    @time ShapiroWilk.OrderStatisticsArblib._precompute_ψ(n, prec=prec, R=18.0)
end
```
then
```julia
exponents = map(5:30) do n
    prec=128
    qs=0.005:0.005:0.995

    OS = ShapiroWilk.OrderStatisticsArblib.NormOrderStatistic(n, prec=prec)
    @info "Computing coefficients..."
    @time swfl = ShapiroWilk.SWCoeffs(n, Float64.(ShapiroWilk.SWCoeffs(OS)))

    @info "Sampling W-statistic..."
    @time W_sample = sort!(sample_Wstatistic(swfl, 10_000_000), rev=true)

    @info "Fitting power-transform for normality..."
    @time res = let qs=qs
        normal_qs=quantile.(Normal(0,1), qs)

        optimize(
            x->normal_power_fit(first(x), W_sample, qs, normal_qs),
            0.01,
            0.5,
            Brent(),
        )
    end
    @info "For n=$n optimization yielded" res

    Optim.minimizer(res)
end
```
yields the normalizing exponents, i.e. `τ` such that `(1-W)^τ ~ N(µ, σ)`. fitting of `τ=τ(n)` is straightforward, e.g.
```julia
┌ Info: For n=30 optimization yielded
│   res =
│    Results of Optimization Algorithm
│     * Algorithm: Brent's Method
│     * Search Interval: [0.010000, 0.500000]
│     * Minimizer: 5.386236e-02
│     * Minimum: 6.233191e-02
│     * Iterations: 12
│     * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true
└     * Objective Function Calls: 13
```

The whole table of these `τ`s look as follows
 
![tau_fit_5:30](https://user-images.githubusercontent.com/1666118/104092204-a7d65600-5282-11eb-936a-1d02e0abaa00.png)




BioTurboNick: Was just looking into Shapiro-Wilk, am I reading it right that it's pretty close?
kalmarek: @BioTurboNick You can already use code from this branch directly (well, the only thing that doesn't work is censored samples) for an implementation based on Royston approximation of S-W test; (well updating to modern HypothesisTests will require a bit coding, but not much).

What is (relatively) close is a new approximation that I'm doing in my spare time ;) Whether it will make a difference, I don't know, but it will be much closer to "mathematical correctness" (whatever that means:).
kalmarek: @andreasnoack or anyone else in this thread ;) I'd need your advice.

I have the SW-coeffs computed up to `eps(Float64)` precision for `n=3:189` (I could go further, but I wasn't patient enough ;).

what is left is the critical levels for the `W`-statistic, i.e. `wstatcdf`. Royston/Parrish use different normalizing transforms (functions that take W to Z ~N(0, 1))

* for `n∈4:11` approximate/interpolate pdf/cdf of `W=W_n` using MC. Royston in 1993 suggests `-ln(γ-ln(1-W))` as normalizing transform, but this seems very crude, given that pdf seems  non-differentiable. I'm not even sure this is taken into the account in the existing Shapiro-Wilk implementations.
*  find a normalizing transform for `n∈12:2000` and estimate its parameters based on MC for `n∈12:189`. Here there are two candidates: `W → (ln(1-W) - µ)/σ` (Royston) or `W → ((1-W)^t - µ)/σ` (Parrish):

After looking closely at MC samples I arrived at this:
 - `n = 3` → analytic solution
 - `n = 4:6` → MC + Piecewise fit ?? 
 - `n = 7:28` → Parrishes power transform `(1-W)^t ~ N(µₙ,σₙ)`: we need to store `(tₙ, µₙ, σₙ)`
 - `n >= 29` → Roystons lognormal `log(1-W) → N(µ′ₙ, σ′ₙ)`: fit biexps/someother models to µ′, and σ′ on `n = 29:189` and use it for all n.

(The transition Parrish → Royston is determined based on MSE on ~200 uniformly distributed quantiles)

What do you think?

I asked about the shape of the W-distributions on discourse: https://discourse.julialang.org/t/can-anyone-tell-me-what-is-this-distribution/53726 but got no input ;)
kalmarek: ok, it seems that I was fitting Normal distribution wrong;) Box-Cox transform `W → (1-W)^λ/λ` for `n in 7:189` seems superior to shifted logarithmic, when optimizing log likelihood (currently using `BoxCoxTrans.jl`).

using `1_000_000` MC samples of `W`-statistic taken from normal (for each of `n∈7:189`) I computed `λ` (the exponent of Box-Cox transform), `µ` and `σ` (the mean and std of the transformed sample of `W`) as functions of `n`
* `inverse(t, a₀, a₁, γ) = a₀ + a₁/(t - γ)`
* `logarithmic(t, a₀, a₁, γ) = a₀ + a₁*log(t - γ)`.
The models (green) were fit using least squares for `n=7:140` (blue) and you can see how well they fit the "test" data (in red). The bottom plot is normalized residuals with their `mean` and `std`.

![λ_7:140](https://user-images.githubusercontent.com/1666118/109663643-a29fd400-7b6c-11eb-8101-3d69b9036239.png)
![µ_7:140](https://user-images.githubusercontent.com/1666118/109663654-a5022e00-7b6c-11eb-905b-f4a6ad712ff0.png)
![σ_7:140](https://user-images.githubusercontent.com/1666118/109663672-a7fd1e80-7b6c-11eb-8858-85cae76f6876.png)



[λ_7:140.pdf](https://github.com/JuliaStats/HypothesisTests.jl/files/6069308/_7.140.pdf)
[µ_7:140.pdf](https://github.com/JuliaStats/HypothesisTests.jl/files/6069310/_7.140.pdf)
[σ_7:140.pdf](https://github.com/JuliaStats/HypothesisTests.jl/files/6069313/_7.140.pdf)

@andreasnoack As I really have very little knowledge about statistics/fitting models/parameters estimation I need someone knowledgeable to have a look at this and provide a stamp of approval ;)
jarredclloyd: @kalmarek @andreasnoack just curious if this is still being looked into/implemented
kalmarek: @jarredclloyd this kind of stalled because I was not able to produce a meaningful fits (except eyeballing plots and saying: this looks seriously gooood...). If you're willing to help (and have more experience with this) be my guest
jarredclloyd: Definitely willing to help, but I think the specific problem you are trying to overcome might be beyond my level of knowledge currently.
Correct me if I'm wrong, but the implementation here has a working version consistent with other languages?
I think it would be great to have this implemented into HypothesisTest so that Julia has a functional Shapiro-Wilk test, especially given all the fantastic work you've done on this. From there we could attempt to resolve the fitting issue (I can ask some of my colleagues who might be able to help). 
So, what needs to be completed to have a functional version of Shapiro-Wilk that can be merged into the main branch?
kalmarek: @jarredclloyd an old todo is here https://github.com/JuliaStats/HypothesisTests.jl/pull/124#issuecomment-439651770;

Minimal set:
1) produce more tests (optional)
2) write documentation
(what is implemented follows the paper of Royston precisely)
3) Resolve the conflict (trivial) + update to the newest version of the HT.

Extended set:
1) perform computations for accurate SWCoeffs (I computed those for n≤400, done)
2) extrapolate these for larger n by finding corrections to Weisberg-Bingham method (for a few first + normalizing the rest)
3) fit a model of the normalizing transform for all n based on the first n=400
4) produce documentation for this.

Can you try 3) ?
codecov-commenter: # [Codecov](https://codecov.io/gh/JuliaStats/HypothesisTests.jl/pull/124?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats) Report
Base: **93.67**% // Head: **92.39**% // Decreases project coverage by **`-1.27%`** :warning:
> Coverage data is based on head [(`c269e66`)](https://codecov.io/gh/JuliaStats/HypothesisTests.jl/pull/124?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats) compared to base [(`a891451`)](https://codecov.io/gh/JuliaStats/HypothesisTests.jl/commit/a89145162648df8b37303894211c0ce7f197f065?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats).
> Patch coverage: 65.00% of modified lines in pull request are covered.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master     #124      +/-   ##
==========================================
- Coverage   93.67%   92.39%   -1.28%     
==========================================
  Files          28       29       +1     
  Lines        1723     1802      +79     
==========================================
+ Hits         1614     1665      +51     
- Misses        109      137      +28     
```


| [Impacted Files](https://codecov.io/gh/JuliaStats/HypothesisTests.jl/pull/124?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats) | Coverage Δ | |
|---|---|---|
| [src/HypothesisTests.jl](https://codecov.io/gh/JuliaStats/HypothesisTests.jl/pull/124/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats#diff-c3JjL0h5cG90aGVzaXNUZXN0cy5qbA==) | `83.09% <ø> (-0.24%)` | :arrow_down: |
| [src/shapirowilk.jl](https://codecov.io/gh/JuliaStats/HypothesisTests.jl/pull/124/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats#diff-c3JjL3NoYXBpcm93aWxrLmps) | `65.00% <65.00%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats)

</details>

[:umbrella: View full report at Codecov](https://codecov.io/gh/JuliaStats/HypothesisTests.jl/pull/124?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaStats).

kalmarek: @andreasnoack, @nalimilan or other maintainers ;) This seems ready for merging.

# What it does:
* Implements an approximate Shapiro-Wilk normality test following Royston 1992
* From `swilk.f` I copied constants which
  - define polynomial approximations of sw-coefficients
  - define normalizing transforms for the W-statistic.
* everything else is a clean-room implementation from the article directly.

This is **not** exactly binary compatible with e.g. numpy which uses `swilk.f`, since there all coeffs are computed as `Float32`s. here even for `Float32` data we use `Float64` precision.

## What it does not:
- [x] re-compute swcoeffs exactly (done for N < 400)
- [ ] update Royston to modern computing standards (partially done)
- [ ] fit new normal transforms

# Speed/benchmarks:
```julia
using Random
Random.seed!(1)
k = 5000
X = sort(randn(k));
A = zeros(k÷2);
```

## Fortran `swilk.f`:
```julia
julia> W, pval, _, A = @btime swilkfort($X)
  77.942 μs (2 allocations: 19.61 KiB)
julia> W, pval, _, A = @btime swilkfort!($X, $A, false)
  75.149 μs (0 allocations: 0 bytes)
(0.9997552011806234, 0.8593080439083717, 0, [...])

```
## This pull:
```
julia> @btime ShapiroWilkTest($X)
  40.647 μs (2 allocations: 19.61 KiB)
Shapiro-Wilk normality test
---------------------------
Population details:
    parameter of interest:   Squared correlation of data and SWCoeffs (W)
    value under h_0:         1.0
    point estimate:          0.999755

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.8593

Details:
    number of observations: 5000
    censored ratio:         0.0
    W-statistic:            0.999755

julia> @btime ShapiroWilkTest($X, SWc=$swc, is_sorted=true);
  5.727 μs (0 allocations: 0 bytes)
```

This has been waiting long enough ;) Many thanks to @jarredclloyd for the latest motivation!

kalmarek: @andreasnoack a friendly ping ;)
kalmarek: @andreasnoack a friendly ping ;)
kalmarek: @andreasnoack a friendly ping ;)
kalmarek: the new benchmarks:
```julia
julia> let N = 5000
           Random.seed!(1234)
           swc = HypothesisTests.SWCoeffs(N)
           sample = randn(N)
           @info "not sorted"
           @btime ShapiroWilkTest($sample, $swc)
           @btime ShapiroWilkTest($sample, $swc, sample_sorted=false)
           sort!(sample)
           @info "sorted"
           @btime ShapiroWilkTest($sample, $swc, sample_sorted=false)
           @btime ShapiroWilkTest($sample, $swc)
           @btime ShapiroWilkTest($sample, $swc, sample_sorted=true)
           A = similar(sample, N÷2)
           @info "swilkfort"
           @btime swilkfort($sample)
           @btime swilkfort!($sample, $A, true)
           @btime swilkfort!($sample, $(A), false)
       end;
```
result in:
```
[ Info: not sorted
  171.111 μs (2 allocations: 39.11 KiB)
  177.118 μs (2 allocations: 39.11 KiB)
[ Info: sorted
  31.777 μs (2 allocations: 39.11 KiB)
  12.292 μs (0 allocations: 0 bytes)
  5.797 μs (0 allocations: 0 bytes)
[ Info: swilkfort
  84.158 μs (2 allocations: 19.61 KiB)
  82.762 μs (0 allocations: 0 bytes)
  82.343 μs (0 allocations: 0 bytes)
```
(swilkfort needs sorted data)


kalmarek: @nalimilan it seems that somebody did not follow semver (`StatsAPI`?) exporting a previously unexported symbol (`pvalue`) requires major version bump; Should I merge `master` here?;
With old manifest tests continue to be green locally. 
nalimilan: > @nalimilan it seems that somebody did not follow semver (`StatsAPI`?) exporting a previously unexported symbol (`pvalue`) requires major version bump; Should I merge `master` here?;
> With old manifest tests continue to be green locally.

Actually, exporting new symbols isn't considered as breaking (otherwise packages would have to tag new major versions all the time). It's actually HypothesisTests which should have used `import StatsAPI: ...` instead of `using StatsAPI`, but most packages do that for historical reasons. Anyway, yes it's better to rebase/merge against master to fix that.

Patch:
diff --git a/Project.toml b/Project.toml
index 8501522..fcdef39 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "HypothesisTests"
 uuid = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
-version = "0.11.1"
+version = "0.11.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/docs/src/nonparametric.md b/docs/src/nonparametric.md
index 048a0ba..889b976 100644
--- a/docs/src/nonparametric.md
+++ b/docs/src/nonparametric.md
@@ -81,3 +81,8 @@ ApproximatePermutationTest
 ```@docs
 FlignerKilleenTest
 ```
+
+## Shapiro-Wilk test
+```@docs
+ShapiroWilkTest
+```
diff --git a/src/kolmogorov_smirnov.jl b/src/kolmogorov_smirnov.jl
index 7be7d05..2413c38 100644
--- a/src/kolmogorov_smirnov.jl
+++ b/src/kolmogorov_smirnov.jl
@@ -63,8 +63,8 @@ sample is not drawn from `d`.
 
 Implements: [`pvalue`](@ref)
 """
-function ExactOneSampleKSTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)
-    if !allunique(x)
+function ExactOneSampleKSTest(x::AbstractVector{T}, d::UnivariateDistribution) where T<:Real
+    if length(x) > length(unique(x))
         @warn("This test is inaccurate with ties")
     end
 
@@ -107,8 +107,8 @@ that the sample is not drawn from `d`.
 
 Implements: [`pvalue`](@ref)
 """
-function ApproximateOneSampleKSTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)
-    if !allunique(x)
+function ApproximateOneSampleKSTest(x::AbstractVector{T}, d::UnivariateDistribution) where T<:Real
+    if length(x) > length(unique(x))
         @warn("This test is inaccurate with ties")
     end
 
diff --git a/src/shapiro_wilk.jl b/src/shapiro_wilk.jl
new file mode 100644
index 0000000..8ccb8f7
--- /dev/null
+++ b/src/shapiro_wilk.jl
@@ -0,0 +1,252 @@
+export ShapiroWilkTest
+
+#=
+From:
+PATRICK ROYSTON
+Approximating the Shapiro-Wilk W-test for non-normality
+*Statistics and Computing* (1992) **2**, 117-119
+DOI: [10.1007/BF01891203](https://doi.org/10.1007/BF01891203)
+=#
+
+# Coefficients from Royston (1992)
+for (s, c) in [(:C1, [0.0, 0.221157, -0.147981, -2.07119, 4.434685, -2.706056]),
+               (:C2, [0.0, 0.042981, -0.293762, -1.752461, 5.682633, -3.582633]),
+               (:C3, [0.5440, -0.39978, 0.025054, -0.0006714]),
+               (:C4, [1.3822, -0.77857, 0.062767, -0.0020322]),
+               (:C5, [-1.5861, -0.31082, -0.083751, 0.0038915]),
+               (:C6, [-0.4803, -0.082676, 0.0030302]),
+               (:C7, [0.164, 0.533]),
+               (:C8, [0.1736, 0.315]),
+               (:C9, [0.256, -0.00635]),
+               (:G, [-2.273, 0.459])]
+    @eval $(Symbol(:__RS92_, s))(x) = Base.Math.@horner(x, $(c...))
+end
+
+"""
+    HypothesisTests.shapiro_wilk_coefs(N::Integer)
+
+Construct a vector of de-correlated expected order statistics for Shapiro-Wilk
+test for a sample of size `N`.
+
+If multiple tests on samples of size `N` are performed, it is beneficial to
+construct and pass a single vector of coefficients to [`ShapiroWilkTest`](@ref).
+"""
+function shapiro_wilk_coefs(N::Integer)
+    if N < 3
+        throw(ArgumentError("N must be greater than or equal to 3, got $N"))
+    elseif N == 3 # exact
+        w = sqrt(2.0) / 2.0
+        return [w,zero(w),-w]
+    else
+        n = div(N, 2)
+        swc = Vector{Float64}(undef, N)
+        m = @view swc[1:n] # store only positive half of swc:
+        # it is (anti-)symmetric; hence '2' factor below
+        for i in eachindex(m)
+            # Weisberg&Bingham 1975 statistic
+            m[i] = -quantile(Normal(), (i - 3 / 8) / (N + 1 / 4))
+        end
+        mᵀm = 2sum(abs2, m)
+        x = 1 / sqrt(N)
+        a₁ = m[1] / sqrt(mᵀm) + __RS92_C1(x) # aₙ = cₙ + (...)
+        if N ≤ 5
+            ϕ = (mᵀm - 2m[1]^2) / (1 - 2a₁^2)
+            m .= m ./ sqrt(ϕ) # A, but reusing m to save allocs
+            m[1] = a₁
+        else
+            a₂ = m[2] / sqrt(mᵀm) + __RS92_C2(x) # aₙ₋₁ = cₙ₋₁ + (...)
+            ϕ = (mᵀm - 2m[1]^2 - 2m[2]^2) / (1 - 2a₁^2 - 2a₂^2)
+            m .= m ./ sqrt(ϕ) # A, but reusing m to save allocs
+            m[1], m[2] = a₁, a₂
+        end
+
+        for i in 1:n
+            swc[N-i+1] = -swc[i]
+        end
+        if isodd(N)
+            swc[n+1] = zero(eltype(swc))
+        end
+        return swc
+    end
+end
+
+function unsafe_swstat(X::AbstractVector{<:Real}, A::AbstractVector{<:Real})
+    _begin = firstindex(A)
+    AX = @inbounds dot(view(A, _begin:(_begin + length(X) - 1)), X)
+    m = mean(X)
+    S² = sum(x -> abs2(x - m), X)
+    W = AX^2 / S²
+    return clamp(W, 0, 1) # to guard against numeric errors
+end
+
+struct ShapiroWilkTest <: HypothesisTest
+    coefs::Vector{Float64}  # expectation of order statistics for Shapiro-Wilk test
+    W::Float64              # test statistic
+    censored::Int           # only the smallest N₁ samples were used
+end
+
+testname(::ShapiroWilkTest) = "Shapiro-Wilk normality test"
+function population_param_of_interest(t::ShapiroWilkTest)
+    return ("Squared correlation of data and expected order statistics of N(0,1) (W)",
+            1.0, t.W)
+end
+default_tail(::ShapiroWilkTest) = :left
+censored_ratio(t::ShapiroWilkTest) = t.censored / length(t.coefs)
+
+function show_params(io::IO, t::ShapiroWilkTest, indent)
+    l = 24
+    println(io, indent, rpad("number of observations:", l), length(t.coefs))
+    println(io, indent, rpad("censored ratio:", l), censored_ratio(t))
+    return println(io, indent, rpad("W-statistic:", l), t.W)
+end
+
+function StatsAPI.pvalue(t::ShapiroWilkTest)
+    n = length(t.coefs)
+    W = t.W
+
+    if iszero(censored_ratio(t))
+        if n == 3 # exact by Shapiro&Wilk 1965
+            # equivalent to 6/π * (asin(sqrt(W)) - asin(sqrt(3/4)))
+            return 1 - 6acos(sqrt(W)) / π
+        elseif n ≤ 11 # Royston 1992
+            γ = __RS92_G(n)
+            w = -log(γ - log1p(-W))
+            μ = __RS92_C3(n)
+            σ = exp(__RS92_C4(n))
+        elseif 12 ≤ n # Royston 1992
+            w = log1p(-W)
+            μ = __RS92_C5(log(n))
+            σ = exp(__RS92_C6(log(n)))
+        end
+        return ccdf(Normal(μ, σ), w)
+    else
+        throw(ArgumentError("censored samples not implemented yet"))
+        # to implement censored samples follow Royston 1993 Section 3.3
+    end
+end
+
+"""
+    ShapiroWilkTest(X::AbstractVector{<:Real},
+                    swc::AbstractVector{<:Real}=shapiro_wilk_coefs(length(X));
+                    sorted::Bool=issorted(X),
+                    censored::Integer=0)
+
+Perform a Shapiro-Wilk test of the null hypothesis that the data in vector `X` come from a
+normal distribution.
+
+This implementation is based on the method by Royston (1992).
+The calculation of the p-value is exact for sample size `N = 3`, and for ranges
+`4 ≤ N ≤ 11` and `12 ≤ N ≤ 5000` (Royston 1992) two separate approximations
+for p-values are used.
+
+# Keyword arguments
+The following keyword arguments may be passed.
+* `sorted::Bool=issorted(X)`: to indicate that sample `X` is already sorted.
+* `censored::Integer=0`: to censor the largest samples from `X`
+  (so called upper-tail censoring)
+
+Implements: [`pvalue`](@ref)
+
+!!! warning
+    As noted by Royston (1993), (approximated) W-statistic will be accurate
+    but returned p-values may not be reliable if either of these apply:
+    * Sample size is large  (`N > 2000`) or small (`N < 20`)
+    * Too much data is censored (`censored / N > 0.8`)
+
+# Implementation notes
+* The current implementation DOES NOT implement p-values for censored data.
+* If multiple Shapiro-Wilk tests are to be performed on samples of same
+  size, it is faster to construct `swc = shapiro_wilk_coefs(length(X))` once
+  and pass it to the test via `ShapiroWilkTest(X, swc)` for re-use.
+* For maximal performance sorted `X` should be passed and indicated with
+  `sorted=true` keyword argument.
+
+# References
+Shapiro, S. S., & Wilk, M. B. (1965). An Analysis of Variance Test for Normality
+(Complete Samples). *Biometrika*, 52, 591–611.
+[doi:10.1093/BIOMET/52.3-4.591](https://doi.org/10.1093/BIOMET/52.3-4.591).
+
+Royston, P. (1992). Approximating the Shapiro-Wilk W-test for non-normality.
+*Statistics and Computing*, 2(3), 117–119.
+[doi:10.1007/BF01891203](https://doi.org/10.1007/BF01891203)
+
+Royston, P. (1993). A Toolkit for Testing for Non-Normality in Complete and
+Censored Samples. Journal of the Royal Statistical Society Series D
+(The Statistician), 42(1), 37–43.
+[doi:10.2307/2348109](https://doi.org/10.2307/2348109)
+
+Royston, P. (1995). Remark AS R94: A Remark on Algorithm AS 181: The W-test for
+Normality. *Journal of the Royal Statistical Society Series C
+(Applied Statistics)*, 44(4), 547–551.
+[doi:10.2307/2986146](https://doi.org/10.2307/2986146).
+"""
+function ShapiroWilkTest(sample::AbstractVector{<:Real},
+                         swcoefs::AbstractVector{<:Real}=shapiro_wilk_coefs(length(sample));
+                         sorted::Bool=issorted(sample),
+                         censored::Integer=0)
+    N = length(sample)
+    if N < 3
+        throw(ArgumentError("at least 3 samples are required, got $N"))
+    elseif censored ≥ N
+        throw(ArgumentError("`censored` must be less than `length(sample)`, " *
+                            "got `censored = $censored > $N = length(sample)`"))
+    elseif length(swcoefs) ≠ length(sample)
+        throw(DimensionMismatch("length of sample and Shapiro-Wilk coefficients " *
+                                "differ, got $N and $(length(swcoefs))"))
+    end
+
+    W = if !sorted
+        X = sort!(sample[1:(end - censored)])
+        if abs(last(X) - first(X)) < length(X) * eps()
+            throw(ArgumentError("sample is constant (up to numerical accuracy)"))
+        end
+        unsafe_swstat(X, swcoefs)
+    else
+        X = @view sample[1:(end - censored)]
+        if last(X) - first(X) < length(X) * eps()
+            throw(ArgumentError("sample doesn't seem to be sorted or " *
+                                "is constant (up to numerical accuracy)"))
+        end
+        unsafe_swstat(X, swcoefs)
+    end
+
+    return ShapiroWilkTest(swcoefs, W, censored)
+end
+
+#=
+# To compare with the standard ALGORITHM AS R94 fortran subroutine
+#  * grab scipys (swilk.f)[https://raw.githubusercontent.com/scipy/scipy/main/scipy/stats/statlib/swilk.f];
+#  * compile
+#  ```
+#  gfortran -shared -fPIC -o swilk.so swilk.f
+#  gfortran -fdefault-integer-8 -fdefault-real-8 -shared -fPIC swilk.f -o swilk64.so
+#  ```
+
+for (lib, I, F) in (("./swilk64.so", Int64, Float64),
+                    ("./swilk.so"  , Int32, Float32))
+    @eval begin
+        function swilkfort!(X::AbstractVector{$F}, A::AbstractVector{$F}, computeA=true)
+            X = issorted(X) ? X : sort(X)
+            w, pval = Ref{$F}(0.0), Ref{$F}(0.0)
+            ifault = Ref{$I}(0)
+
+            ccall((:swilk_, $lib),
+                Cvoid,
+                (
+                Ref{Bool},  # INIT if false compute SWCoeffs in A, else use A
+                Ref{$F},    # X    sample
+                Ref{$I},    # N    samples length
+                Ref{$I},    # N1   (upper) uncensored data length
+                Ref{$I},    # N2   length of A
+                Ref{$F},    # A    A
+                Ref{$F},    # W    W-statistic
+                Ref{$F},    # PW   p-value
+                Ref{$I},    # IFAULT error code (see swilk.f for meaning)
+                ),
+                !computeA, X, length(X), length(X), div(N,2), A, w, pval, ifault)
+            return (w[], pval[], ifault[], A)
+        end
+        swilkfort(X::Vector{$F}) = swilkfort!(X, zeros($F, div(length(X),2)))
+    end
+end
+=#


Test patch:
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index e3d12e2..cc33ed8 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -148,5 +148,6 @@ include("diebold_mariano.jl")
 include("clark_west.jl")
 include("white.jl")
 include("var_equality.jl")
+include("shapiro_wilk.jl")
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 18fff3d..5cbbfd7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -30,6 +30,7 @@ include("mann_whitney.jl")
 include("correlation.jl")
 include("permutation.jl")
 include("power_divergence.jl")
+include("shapiro_wilk.jl")
 include("show.jl")
 include("t.jl")
 include("wald_wolfowitz.jl")
diff --git a/test/shapiro_wilk.jl b/test/shapiro_wilk.jl
new file mode 100644
index 0000000..acfba78
--- /dev/null
+++ b/test/shapiro_wilk.jl
@@ -0,0 +1,130 @@
+using HypothesisTests, LinearAlgebra, Test, Random
+using StableRNGs
+
+@testset "Shapiro-Wilk" begin
+    @testset "shapiro_wilk_coefs" begin
+        @test HypothesisTests.shapiro_wilk_coefs(3) == [sqrt(2.0) / 2.0, 0.0, -sqrt(2.0) / 2.0]
+
+        swc = HypothesisTests.shapiro_wilk_coefs(10)
+        @test swc[4] == -swc[7]
+        @test swc[2] == -swc[9]
+
+        swc = HypothesisTests.shapiro_wilk_coefs(11)
+
+        @test swc[6] == 0.0
+        @test swc[5] == -swc[7]
+        @test swc[3] == -swc[9]
+        @test swc[1] == -swc[11]
+
+        #anti-symmetry
+        swc = HypothesisTests.shapiro_wilk_coefs(20)
+        @test all([swc[i] == -swc[end - i + 1] for i in eachindex(swc)])
+
+        # Values obtained by calling `_swilkfort` fortran subroutine directly.
+
+        swc10 = HypothesisTests.shapiro_wilk_coefs(10)
+        res = swc10[1:5] .- [0.573715, 0.32897, 0.214349, 0.122791, 0.0400887]
+        @test norm(res, 1) ≈ 0.0 atol = length(swc10) * eps(Float32)
+
+        swc20 = HypothesisTests.shapiro_wilk_coefs(20)
+        res = swc20[1:10] .-
+              [0.473371, 0.32174, 0.255663, 0.208297, 0.16864, 0.133584, 0.101474,
+               0.0712893, 0.0423232, 0.0140351]
+        @test norm(res, 1) ≈ 0.0 atol = length(swc20) * eps(Float32)
+
+        rng = StableRNG(0x5bca7c69b794f8ce)
+        X = sort(rand(rng, Float64, 10))
+        # W, pval, _ = swilkfort(X)
+        W = 0.9088434774710951
+        # pval = 0.2731410626084226
+        @test HypothesisTests.unsafe_swstat(X, swc10) ≈ W atol = eps(Float32)
+    end
+
+    @testset "Shapiro-Wilk" begin
+        # syntactic tests
+
+        @test_throws ArgumentError ShapiroWilkTest([1, 2])
+        @test_throws ArgumentError("at least 3 samples are required, got 2") ShapiroWilkTest([1, 2], HypothesisTests.shapiro_wilk_coefs(3))
+        @test_throws ArgumentError ShapiroWilkTest([1, 2, 3], censored=4)
+        @test_throws DimensionMismatch ShapiroWilkTest([1, 2, 3],
+                                                       HypothesisTests.shapiro_wilk_coefs(4))
+
+        @test_throws ArgumentError("sample doesn't seem to be sorted or is constant (up to numerical accuracy)") ShapiroWilkTest([1,1,1])
+        @test_throws ArgumentError("sample is constant (up to numerical accuracy)") ShapiroWilkTest([1,1,1], sorted=false)
+
+        t = ShapiroWilkTest([1, 2, 3])
+        @test t.W == 1.0
+        @test pvalue(t) == 1.0
+
+        @test_throws ArgumentError("censored samples not implemented yet") pvalue(ShapiroWilkTest(1:4, censored=1))
+
+        str = sprint(show, t)
+        @test str ==
+              """Shapiro-Wilk normality test
+              ---------------------------
+              Population details:
+                  parameter of interest:   Squared correlation of data and expected order statistics of N(0,1) (W)
+                  value under h_0:         1.0
+                  point estimate:          1.0
+
+              Test summary:
+                  outcome with 95% confidence: fail to reject h_0
+                  one-sided p-value:           1.0000
+
+              Details:
+                  number of observations: 3
+                  censored ratio:         0.0
+                  W-statistic:            1.0
+              """
+
+        # testing different cases of N
+        for N in (3, 5, 11, 12)
+            rng = StableRNG(0x5bca7c69b794f8ce)
+            X = sort(randn(rng, N))
+            t = ShapiroWilkTest(X; sorted=true)
+
+            # analytic properties from Shapiro-Wilk 1965:
+            # Lemma 1: Scale and origin invariance:
+            scale, shift = rand(rng, 2)
+            @test t.W ≈ ShapiroWilkTest(X .+ shift).W
+            @test t.W ≈ ShapiroWilkTest(scale .* X).W
+            @test t.W ≈ ShapiroWilkTest(scale .* X .+ shift).W
+            # Lemma 2, 3: upper and lower bounds
+            @test N * t.coefs[1]^2 / (N - 1) ≤ t.W ≤ 1.0
+
+            # test the computation of pvalue in those cases
+            @test pvalue(t) ≥ 0.05
+        end
+
+        @testset "Worked Example" begin
+            # **Worked Example** (Section 4) from
+            # PATRICK ROYSTON Approximating the Shapiro-Wilk W-test for non-normality
+            # *Statistics and Computing* (1992) **2**, 117-119
+
+            X = [48.4, 49.0, 59.5, 59.6, 60.7, 88.8, 98.2, 109.4, 169.1, 227.1]
+            swc = HypothesisTests.shapiro_wilk_coefs(length(X))
+            @test norm(swc[1:5] .- [0.5737, 0.3290, 0.2143, 0.1228, 0.0401], Inf) < 5.0e-5
+            W = HypothesisTests.unsafe_swstat(X, swc)
+            @test W ≈ 0.8078 atol = 2.9e-5
+
+            t = ShapiroWilkTest(X)
+            @test t.W == W
+            @test pvalue(t) ≈ 0.018 atol = 4.7e-5
+            @test pvalue(t) ≈ pvalue(ShapiroWilkTest(X, sorted=true))
+            @test iszero(HypothesisTests.censored_ratio(t))
+            @test length(t.coefs) == length(X)
+
+            # test for un-sorted sample
+            X2 = X[[9,8,2,3,4,5,7,10,1,6]]
+            t2 = ShapiroWilkTest(X2)
+            @test_throws ArgumentError ShapiroWilkTest(X2, sorted=true)
+            @test t2.W ≈ t.W
+            @test pvalue(t2) ≈ pvalue(t)
+            X3 = X[[2,8,9,3,4,5,7,10,1,6]]
+            t3 = ShapiroWilkTest(X3)
+            @test t3.W ≈ t.W
+            @test pvalue(t3) ≈ pvalue(t)
+            @test pvalue(t3) ≠ pvalue(ShapiroWilkTest(X3, sorted=true))
+        end
+    end
+end