Repository: scipy/scipy
Pull request: TST: stats.ortho_group: improve determinant distribution test
#### Reference issue
Closes gh-18272

#### What does this implement/fix?
gh-18272 reported an occasional failure of `test_multivariate.TestOrthoGroup.test_det_and_ortho`. This improves the test and makes it deterministic.

Problem statement:
BUG: stats: occasional failure of `test_multivariate.TestOrthoGroup.test_det_and_ortho`
I see an occasional failure of the test  `scipy.stats.tests.test_multivariate.TestOrthoGroup.test_det_and_ortho`.

Here's a typical example of the failure:
```
=========================================================== FAILURES ===========================================================
______________________________________________ TestOrthoGroup.test_det_and_ortho _______________________________________________
[gw1] linux -- Python 3.10.8 /home/warren/py3.10.8/bin/python3
scipy/stats/tests/test_multivariate.py:1867: in test_det_and_ortho
    assert_array_less([np.nonzero(d < 0)[0].shape[0] for d in dets], [10]*10)
        dets       = array([[-1., -1., -1.,  1., -1., -1.,  1., -1., -1., -1.],
       [-1., -1., -1.,  1., -1.,  1.,  1., -1., -1.,  1.],
...
       [-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.],
       [ 1.,  1., -1., -1.,  1., -1., -1.,  1.,  1., -1.]])
        self       = <scipy.stats.tests.test_multivariate.TestOrthoGroup object at 0x7ff076214df0>
        xs         = [[array([[-0.29334315, -0.95600722],
       [-0.95600722,  0.29334315]]), array([[ 0.42281174, -0.90621754],
       [-...      [-0.76983009,  0.03994489,  0.00614306,  0.29541576, -0.19363484,
         0.26732913, -0.45771008]]), ...], ...]
/home/warren/py3.10.8/lib/python3.10/contextlib.py:79: in inner
    return func(*args, **kwds)
E   AssertionError: 
E   Arrays are not less-ordered
E   
E   Mismatched elements: 1 / 10 (10%)
E   Max absolute difference: 7
E   Max relative difference: 0.7
E    x: array([ 8,  6,  6,  3,  4,  6,  5,  7, 10,  5])
E    y: array([10, 10, 10, 10, 10, 10, 10, 10, 10, 10])
        args       = (<built-in function lt>, [8, 6, 6, 3, 4, 6, ...], [10, 10, 10, 10, 10, 10, ...])
        func       = <function assert_array_compare at 0x7ff079006200>
        kwds       = {'equal_inf': False, 'err_msg': '', 'header': 'Arrays are not less-ordered', 'verbose': True}
        self       = <contextlib._GeneratorContextManager object at 0x7ff079052350>
```
The test checks the determinants of a set of 10 random orthogonal matrices.  It expects at least one to have a negative determinant.  But the probability of having all determinants positive is 1/1024, which is not that small when you consider how often the SciPy tests are run.  It is surprising that we don't see this failure more often.


### SciPy/NumPy/Python version and system information

```shell
1.11.0.dev0+1827.26cd496 1.25.0.dev0+1080.gd79fe4364c sys.version_info(major=3, minor=10, micro=8, releaselevel='final', serial=0)
/home/warren/repos/git/forks/scipy/build-install/lib/python3.10/site-packages/scipy/__config__.py:140: UserWarning: Install `pyyaml` for better output
  warnings.warn("Install `pyyaml` for better output", stacklevel=1)
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "11.3.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "0.29.33",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "11.3.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "11.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.12.1",
      "include directory": "../../../../../py3.10.8/lib/python3.10/site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "linux"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "linux"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas",
      "found": true,
      "version": "0.3.20",
      "detection method": "pkgconfig",
      "include directory": "/usr/include/x86_64-linux-gnu/openblas-pthread/",
      "lib directory": "/usr/lib/x86_64-linux-gnu/openblas-pthread/",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER=1 NO_CBLAS= NO_LAPACK= NO_LAPACKE=1 NO_AFFINITY=1 USE_OPENMP=0 generic MAX_THREADS=64",
      "pc file directory": "/usr/lib/x86_64-linux-gnu/pkgconfig"
    },
    "lapack": {
      "name": "openblas",
      "found": true,
      "version": "0.3.20",
      "detection method": "pkgconfig",
      "include directory": "/usr/include/x86_64-linux-gnu/openblas-pthread/",
      "lib directory": "/usr/lib/x86_64-linux-gnu/openblas-pthread/",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER=1 NO_CBLAS= NO_LAPACK= NO_LAPACKE=1 NO_AFFINITY=1 USE_OPENMP=0 generic MAX_THREADS=64",
      "pc file directory": "/usr/lib/x86_64-linux-gnu/pkgconfig"
    }
  },
  "Python Information": {
    "path": "/home/warren/py3.10.8/bin/python3",
    "version": "3.10"
  }
}
```

Hints from issue discussion:
Thanks Warren. Sounds like it would have been better to manually set some matrices instead of using random values without meaning.
In this case, we're testing a statistical property of the `rvs()` method, so generating a random sample is appropriate.  But we need a smarter test.  For example: generate more samples (maybe 200 or more), and then apply the binomial test to the number of positive determinants, with the threshold of the p-value for not rejecting the null hypothesis set extremely small (e.g. 1e-16) , because we know this test will be repeated a lot.
> It is surprising that we don't see this failure more often.

Does the seed being set in the first test of the class make the behavior deterministic if the tests are run serially?

https://github.com/scipy/scipy/blob/e574cbcabf8d25955d1aafeed02794f8b5f250cd/scipy/stats/tests/test_multivariate.py#L1822-L1825

I don't know that they are always run serially in CI; just suggesting a possible explanation.

Fortunately, the determinant seems equally likely to be positive or negative.
<details>

```python3
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
rng = np.random.default_rng(246985373091345)
dist = stats.ortho_group(dim=20)  # tried a few `dim` {2, 5, 10, 20}; qualitatively same results
rvs = dist.rvs(size=10000, random_state=rng)
dets = np.linalg.det(rvs)
plt.hist(dets, bins=[-1.1, -.1, 0.1, 1.1])
```

</details>

<img width="1138" alt="image" src="https://user-images.githubusercontent.com/6570539/231372201-9691f2b8-9bd9-4350-b10e-7e566bb4cd1a.png">


> But we need a smarter test.

Definitely. 

> with the threshold of the p-value for not rejecting the null hypothesis set extremely small (e.g. 1e-16) , because we know this test will be repeated a lot.

But I'm fairly confident that the determinant is positive or negative with equal probability right now, and a test with such a small threshold would be less sensitive to small changes.

`ortho_group` `rvs` speed was substantially improved last year, so we can use a sample much larger than 200.
```python3
dist = stats.ortho_group(dim=20)
%timeit rvs = dist.rvs(size=10000, random_state=rng)
# 274 ms ± 10.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

Also, I'm a fan of stochastic (e.g. `np.random.default_rng()`, no seed) tests in principle, but not with the current CI system or for a test of this sort. So I would suggest making the test deterministic but draw as large a sample as we are willing (considering test time), and use a more moderate threshold on the p-value. Similarly, we could check that a confidence interval with a moderately _loose_ confidence contains 0.5:

```python3
k = np.sum(dets > 0)
n = len(dets)
res = stats.binomtest(k, n)
res.proportion_ci(confidence_level=0.95)
# ConfidenceInterval(low=0.4945504780403626, high=0.5142469581380023)
```

Commit message: STY: stats.tests.test_multivariate: remove blank line whitespace

[ci skip]

Patch to review:
diff --git a/scipy/stats/tests/test_multivariate.py b/scipy/stats/tests/test_multivariate.py
index 1964728cd6cf..6ffbb5780ff4 100644
--- a/scipy/stats/tests/test_multivariate.py
+++ b/scipy/stats/tests/test_multivariate.py
@@ -532,7 +532,7 @@ def test_frozen(self):
         assert_allclose(norm_frozen.cdf(x), multivariate_normal.cdf(x, mean, cov))
         assert_allclose(norm_frozen.logcdf(x),
                         multivariate_normal.logcdf(x, mean, cov))
-    
+
     @pytest.mark.parametrize(
         'covariance',
         [
@@ -1860,18 +1860,25 @@ def test_det_and_ortho(self):
         dets = np.array([[np.linalg.det(x) for x in xx] for xx in xs])
         assert_allclose(np.fabs(dets), np.ones(dets.shape), rtol=1e-13)
 
-        # Test that we get both positive and negative determinants
-        # Check that we have at least one and less than 10 negative dets in a sample of 10. The rest are positive by the previous test.
-        # Test each dimension separately
-        assert_array_less([0]*10, [np.nonzero(d < 0)[0].shape[0] for d in dets])
-        assert_array_less([np.nonzero(d < 0)[0].shape[0] for d in dets], [10]*10)
-
         # Test that these are orthogonal matrices
         for xx in xs:
             for x in xx:
                 assert_array_almost_equal(np.dot(x, x.T),
                                           np.eye(x.shape[0]))
 
+    @pytest.mark.parametrize("dim", [2, 5, 10, 20])
+    def test_det_distribution_gh18272(self, dim):
+        # Test that positive and negative determinants are equally likely.
+        rng = np.random.default_rng(6796248956179332344)
+        dist = ortho_group(dim=dim)
+        rvs = dist.rvs(size=5000, random_state=rng)
+        dets = np.linalg.det(rvs)
+        k = np.sum(dets > 0)
+        n = len(dets)
+        res = stats.binomtest(k, n)
+        low, high = res.proportion_ci(confidence_level=0.95)
+        assert low < 0.5 < high
+
     def test_haar(self):
         # Test that the distribution is constant under rotation
         # Every column should have the same distribution


Code review comments:
diff_hunk: @@ -1860,18 +1860,25 @@ def test_det_and_ortho(self):
         dets = np.array([[np.linalg.det(x) for x in xx] for xx in xs])
         assert_allclose(np.fabs(dets), np.ones(dets.shape), rtol=1e-13)
 
-        # Test that we get both positive and negative determinants
-        # Check that we have at least one and less than 10 negative dets in a sample of 10. The rest are positive by the previous test.
-        # Test each dimension separately
-        assert_array_less([0]*10, [np.nonzero(d < 0)[0].shape[0] for d in dets])
-        assert_array_less([np.nonzero(d < 0)[0].shape[0] for d in dets], [10]*10)
-
         # Test that these are orthogonal matrices
         for xx in xs:
             for x in xx:
                 assert_array_almost_equal(np.dot(x, x.T),
                                           np.eye(x.shape[0]))
 
+    @pytest.mark.parametrize("dim", [2, 5, 10, 20])
+    def test_det_distribution_gh18272(self, dim):
+        # Test that positive and negative determinants are equally likely.
+        rng = np.random.default_rng(6796248956179332344)
+        dist = ortho_group(dim=dim)
+        rvs = dist.rvs(size=5000, random_state=rng)
+        dets = np.linalg.det(rvs)
path: scipy/stats/tests/test_multivariate.py
text: @author:
```suggestion
        dets = scipy.linalg.det(rvs)
```
@user1 it should just be a drop-in replacement, right? Luckily, `scipy.linalg` is already imported at the top, so feel free to commit if this looks right.

@user1:
Yes, fingers crossed, but you need to rebase I think

@author:
IIRC it's automatically done in CI before tests run. (I seem to remember CircleCI not doing that, but we don't run regular tests there.) Anyway, I committed and tests passed. 

@user1:
Phew! 😅