Repository: scipy/scipy
Pull request: BUG: stats: fix multivariate_hypergeom.rvs method
#### Reference issue

fixes #16171 

#### What does this implement/fix?

Fixes multivariate_hypergeom.rvs method when the last two or more elements of the population are 0.


Problem statement:
BUG: scipy.stats.multivariate_hypergeom.rvs raises ValueError when at least the last two populations are 0
### Describe your issue.

I ran into a strange issue when sampling using Scipy's multivariate hypergeometric distribution. Calling `scipy.stats.multivariate_hypergeom.rvs` on a list of populations raises the error message `ValueError: ngood + nbad < nsample`, while calling the same function on the sorted array as input does not.
Specfically, the error seems to be raised when at least the last two values of the population are 0. 

### Reproducing Code Example

```python
import scipy.stats
r = scipy.stats.multivariate_hypergeom.rvs
res1 = r([0, 20], 5) # Returns [0, 5]
res2 = r([0, 0, 20], 5) # Returns [0, 0, 5]
res3 = r([0, 0, 20, 0], 5) # # Returns [0, 0, 5, 0]
res4 = r([0, 0, 20, 0, 0], 5) # # Raises error
```


### Error message

```shell
Traceback (most recent call last):
  File "<string>", line 6, in <module>
  File "/home/alexander/miniconda3/envs/py37/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 4726, in rvs
    size=size))
  File "mtrand.pyx", line 3839, in numpy.random.mtrand.RandomState.hypergeometric
ValueError: ngood + nbad < nsample
```


### SciPy/NumPy/Python version information

1.7.3 1.21.5 sys.version_info(major=3, minor=10, micro=4, releaselevel='final', serial=0)

Hints from issue discussion:
Hi @Alexamk, thank you for reporting. Indeed I can reproduce and agree that this behaviour is confusing.
The problem is here:

https://github.com/scipy/scipy/blob/a6ee381ec61f89e8b43631bb7f3e5ee30e6bda47/scipy/stats/_multivariate.py#L4875-L4878

When n is 0, the result on the right is zero (because of the multiplication by `(n != 0)`), but the function `random_state.hypergeometric()` is still being evaluated, and in this case the arguments to `random_state.hypergeometric()` are not valid.
@WarrenWeckesser do you mean we just need to add a check to see if the right element is 0?
My comment about the arguments being invalid needs some qualification.

If we could rely on `random_state` being a `Generator` instance (i.e. if we could ensure that we were using the new NumPy random API instead of the legacy API), we could replace that line with:

```
            rvs[..., c] = random_state.hypergeometric(m[..., c], rem, n,
                                                      size=size)
```
`n` might be 0, but `Generator.hypergeometric` will do the right thing.  For example,
```
In [8]: rng = np.random.default_rng()

In [9]: rng.hypergeometric(0, 1, 0)
Out[9]: 0

In [10]: rng.hypergeometric(1, 0, 0)
Out[10]: 0

In [11]: rng.hypergeometric(0, 0, 0)
Out[11]: 0
```
The reason for the shenanigans with `n != 0` and `n == 0` is that we have to support the legacy random API (and we have to support broadcasting, so we can't use a simple `if` statement), and `np.random.hypergeometric` doesn't allow `n == 0`:
```
In [16]: np.random.hypergeometric(0, 1, 0)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [16], in <module>
----> 1 np.random.hypergeometric(0, 1, 0)

File mtrand.pyx:3871, in numpy.random.mtrand.RandomState.hypergeometric()

File _common.pyx:919, in numpy.random._common.disc()

File _common.pyx:436, in numpy.random._common.check_constraint()

ValueError: nsample < 1 or nsample is NaN
```
The code as written fixes that problem by multiplying the result of `random_state.hypergeometric()` by `(n != 0)`; this ensures that wherever `n` is 0, the result is 0.  It passes the argument `n + (n == 0)` to `random_state.hypergeometric()` to prevent the `ValueError` from being raised in that case.  However, we missed the case where *both* of the first two arguments to `hypergeometric()` are 0.  In that case, our code does, in effect:

```
In [18]: np.random.hypergeometric(0, 0, 1)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [18], in <module>
----> 1 np.random.hypergeometric(0, 0, 1)

File mtrand.pyx:3870, in numpy.random.mtrand.RandomState.hypergeometric()

ValueError: ngood + nbad < nsample

```
One way to fix this is to add `(n == 0)` to `rem` in addition to `n`:
```
            n0mask = n == 0
            rvs[..., c] = (~n0mask *
                           random_state.hypergeometric(m[..., c],
                                                       rem + n0mask,
                                                       n + n0mask,
                                                       size=size))
```
The proposed fix sounds good to me @WarrenWeckesser. Do you want to propose the PR @WarrenWeckesser or @Alexamk? If not, I can take this up.
@tirthasheshpatel, I probably wouldn't get to a PR right away, so if you have time now, go for it!

Commit message: BUG: stats: fix multivariate_hypergeom.rvs method

Patch to review:
diff --git a/scipy/stats/_multivariate.py b/scipy/stats/_multivariate.py
index 204cb2e754e4..1cdfc36d0369 100644
--- a/scipy/stats/_multivariate.py
+++ b/scipy/stats/_multivariate.py
@@ -4872,9 +4872,11 @@ def rvs(self, m, n, size=None, random_state=None):
         # https://github.com/numpy/numpy/pull/13794
         for c in range(m.shape[-1] - 1):
             rem = rem - m[..., c]
-            rvs[..., c] = ((n != 0) *
-                           random_state.hypergeometric(m[..., c], rem,
-                                                       n + (n == 0),
+            n0mask = n == 0
+            rvs[..., c] = (~n0mask *
+                           random_state.hypergeometric(m[..., c],
+                                                       rem + n0mask,
+                                                       n + n0mask,
                                                        size=size))
             n = n - rvs[..., c]
         rvs[..., m.shape[-1] - 1] = n
diff --git a/scipy/stats/tests/test_multivariate.py b/scipy/stats/tests/test_multivariate.py
index ce6ee7301c09..f86051913411 100644
--- a/scipy/stats/tests/test_multivariate.py
+++ b/scipy/stats/tests/test_multivariate.py
@@ -1976,6 +1976,17 @@ def test_rvs_broadcasting(self):
         rvs = rv.rvs(size=(1000, 2), random_state=123)
         assert_allclose(rvs.mean(0), rv.mean(), rtol=1e-2)
 
+    @pytest.mark.parametrize('m, n', (
+        ([0, 0, 20, 0, 0], 5), ([0, 0, 0, 0, 0], 0), ([0, 0], 0),
+        ([0, 0], 0), ([0], 0)
+    ))
+    def test_rvs_gh16171(self, m, n):
+        res = multivariate_hypergeom.rvs(m, n)
+        m = np.asarray(m)
+        res_ex = m.copy()
+        res_ex[m != 0] = n
+        assert_equal(res, res_ex)
+
     @pytest.mark.parametrize(
         "x, m, n, expected",
         [


Code review comments:
diff_hunk: @@ -1976,6 +1976,17 @@ def test_rvs_broadcasting(self):
         rvs = rv.rvs(size=(1000, 2), random_state=123)
         assert_allclose(rvs.mean(0), rv.mean(), rtol=1e-2)
 
+    @pytest.mark.parametrize('m, n', (
+        ([0, 0, 20, 0, 0], 5), ([0, 0, 0, 0, 0], 0), ([0, 0], 0),
path: scipy/stats/tests/test_multivariate.py
text: @user1:
You already have this case in your test matrix. Or you had another case in mind maybe?
```suggestion
        ([0, 0, 20, 0, 0], 5), ([0, 0, 0, 0, 0], 0),
```

@author:
> Or you had another case in mind maybe?

I just wanted to make sure that small arrays with all 0s don't error out.

@author:
Oh, I see that this test case is just a duplicate of the one below. Will remove it as you have suggested here. Thanks for noticing!