Repository: scipy/scipy
Pull request: MAINT: stats.shapiro: override p-value when len(x)==3
#### Reference issue
Closes gh-18322

#### What does this implement/fix?
gh-18322 reported that `shapiro` can produce negative p-values when `len(x)==3`. This PR overrides the p-value reported by `swilk` in this case.

*Lint failure unrelated. That will be fixed when gh-18321 is reverted (i.e. gh-16996 comes back)*

Problem statement:
BUG: scipy.stats.shapiro gives a negative pvalue
### Describe your issue.

p-values are supposed to be in [0,1], but I get a negative p-value using the scipy shapiro.
See the example code below.

### Reproducing Code Example

```python
from scipy.stats import shapiro
x = [0.6931471805599453, 0.0, 0.0]
res = shapiro(x)
res.pvalue # gives -6.830188112871838e-07, which is negative...
```


### Error message

```shell
# no error message, just a negative p-value
```


### SciPy/NumPy/Python version and system information

```shell
1.10.1 1.23.5 sys.version_info(major=3, minor=9, micro=15, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: cmake
    found: true
    include directory: unknown
    lib directory: unknown
    name: OpenBLAS
    openblas configuration: unknown
    pc file directory: unknown
    version: 0.3.18
  lapack:
    detection method: cmake
    found: true
    include directory: unknown
    lib directory: unknown
    name: OpenBLAS
    openblas configuration: unknown
    pc file directory: unknown
    version: 0.3.18
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 12.0.0
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 12.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 0.29.33
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 4.9.0
  pythran:
    include directory: /private/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T/pip-build-env-x9bi86_d/overlay/lib/python3.9/site-packages/pythran
    version: 0.12.1
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: darwin
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: darwin
Python Information:
  path: /private/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T/cibw-run-1cu21evt/cp39-macosx_x86_64/build/venv/bin/python
  version: '3.9'
```

Hints from issue discussion:
The null distribution of the Shapiro-Wilk test is usually estimated using a Monte Carlo approach, but for $n=3$ there is an [exact solution](https://math.mit.edu/~rmd/465/shapiro.pdf). 
![image](https://user-images.githubusercontent.com/6570539/233172167-449fb54b-e458-4743-a948-941c2e58efa3.png)

I think we want to override whatever `shapiro` is doing in this case.

<details>

```python3
import numpy as np
from scipy import stats
x = [1, 2, 5]
res = stats.shapiro(x)

dist = stats.beta(a=0.5, b=0.5)
l = dist.cdf(0.75)
r = 1
p = dist.cdf(res.statistic)
pvalue = (p - l)/(r - l)

np.testing.assert_allclose(res.pvalue, pvalue, rtol=1e-5)

# could also integrate the pdf numerically for statistic close to 0.75
from scipy import integrate
def f(x):
  return 3/np.pi*(1-x)**(-0.5)*x**(-0.5)
pvalue2 = integrate.quad(f, 0.75, res.statistic)[0]
np.testing.assert_allclose(pvalue2, pvalue, rtol=1e-5)

# And there's an analytical solution
# https://www.wolframalpha.com/input?i=integrate+%281-x%29%5E%28-1%2F2%29+*+x%5E%28-1%2F2%29+from+3%2F4+to+p
pvalue3 = 1 - 6/np.pi*np.arccos(np.sqrt(res.statistic))
np.testing.assert_allclose(pvalue3, pvalue2, rtol=1e-5)
```

There is an [asymptotic expansion](https://www.wolframalpha.com/input?i=2*arccos%28sqrt%28x%29%29+expansion+near+x%3D3%2F4) that could be used to avoid catastropic cancellation near the left end.

</details>

Looks like the true p-value here is close to `1.31447057e-7`. 
That would fix the case n=3, but are we sure that negative values are only appearing for n=3? Since the underlying Fortran is more than 20 years old, I guess we would have heard by now, but who knows.
Not certain, but $n=3$ is indeed treated differently from other cases (confirmed in `swilk.f`), so I'm content with fixing the reported bug to close this.

Commit message: MAINT: stats.shapiro: override p-value when len(x)==3

Patch to review:
diff --git a/scipy/stats/_morestats.py b/scipy/stats/_morestats.py
index 65afcdc20308..34b16b4bbf4c 100644
--- a/scipy/stats/_morestats.py
+++ b/scipy/stats/_morestats.py
@@ -1881,6 +1881,10 @@ def shapiro(x):
     if N > 5000:
         warnings.warn("p-value may not be accurate for N > 5000.")
 
+    # `swilk` can return negative p-values for N==3; see gh-18322.
+    if N == 3:
+        # Potential improvement: precision for small p-values
+        pw = 1 - 6/np.pi*np.arccos(np.sqrt(w))
     return ShapiroResult(w, pw)
 
 
diff --git a/scipy/stats/tests/test_morestats.py b/scipy/stats/tests/test_morestats.py
index 0ebef783f5d6..98e740332d23 100644
--- a/scipy/stats/tests/test_morestats.py
+++ b/scipy/stats/tests/test_morestats.py
@@ -225,6 +225,22 @@ def test_gh14462(self):
 
         assert_allclose(res, ref, rtol=1e-5)
 
+    def test_length_3_gh18322(self):
+        # gh-18322 reported that the p-value could be negative for input of
+        # length 3. Check that this is resolved.
+        res = stats.shapiro([0.6931471805599453, 0.0, 0.0])
+        assert res.pvalue > 0
+
+        # R `shapiro.test` doesn't produce an accurate p-value in the case
+        # above. Check that the formula used in `stats.shapiro` is not wrong.
+        # options(digits=16)
+        # x = c(-0.7746653110021126, -0.4344432067942129, 1.8157053280290931)
+        # shapiro.test(x)
+        rng = np.random.default_rng(376593478292346598752763)
+        res = stats.shapiro(rng.normal(size=3))
+        assert_allclose(res.statistic, 0.84658770645509)
+        assert_allclose(res.pvalue, 0.2313666489882, rtol=1e-6)
+
 
 class TestAnderson:
     def test_normal(self):


Code review comments:
diff_hunk: @@ -225,6 +225,22 @@ def test_gh14462(self):
 
         assert_allclose(res, ref, rtol=1e-5)
 
+    def test_length_3_gh18322(self):
+        # gh-18322 reported that the p-value could be negative for input of
+        # length 3. Check that this is resolved.
+        res = stats.shapiro([0.6931471805599453, 0.0, 0.0])
+        assert res.pvalue > 0
path: scipy/stats/tests/test_morestats.py
text: @author:
```suggestion
        assert res.pvalue >= 0
```

diff_hunk: @@ -225,6 +225,22 @@ def test_gh14462(self):
 
         assert_allclose(res, ref, rtol=1e-5)
 
+    def test_length_3_gh18322(self):
+        # gh-18322 reported that the p-value could be negative for input of
+        # length 3. Check that this is resolved.
+        res = stats.shapiro([0.6931471805599453, 0.0, 0.0])
+        assert res.pvalue > 0
+
+        # R `shapiro.test` doesn't produce an accurate p-value in the case
+        # above. Check that the formula used in `stats.shapiro` is not wrong.
+        # options(digits=16)
+        # x = c(-0.7746653110021126, -0.4344432067942129, 1.8157053280290931)
+        # shapiro.test(x)
+        rng = np.random.default_rng(376593478292346598752763)
+        res = stats.shapiro(rng.normal(size=3))
path: scipy/stats/tests/test_morestats.py
text: @user1:
Why the use of the RNG here? I would suggest putting this in comment if you want but explicitly having the values there already is better to me.

@author:
To show that the numbers are random.

@author:
If you prefer:
```suggestion
        x = [-0.7746653110021126, -0.4344432067942129, 1.8157053280290931]
        res = stats.shapiro(x)
```

@user1:
I do yes :)