Repository: scipy/scipy
Pull request: BUG: adds warning to scipy.stats.brunnermunzel
#### Reference issue
closes #15843

#### What does this implement/fix?

emits a warning when `scipy.stats.brunnermunzel` cannot compute p-value
because both numerator and denominator are 0, and therefore p-value is nan
also suggests an alternative (using distribution="normal")

~~This change does not change any output/no code is changed. The only change
is emitting a warning in some cases where output p-value is nan~~


also add two tests:

1. `test_brunnermunzel_return_nan` tests that the correct warning is raised when using a t-distribution
2. `test_brunnermunzel_normal_dist` tests that p=0 when using a normal dist with the problematic dataset


Problem statement:
BUG: scipy.stats.brunnermunzel incorrectly returns nan for undocumented reason
### Describe your issue.

```python

x = np.array([1,2,3])
y = np.array([5,6,7,8,9])

scipy.stats.brunnermunzel(x,y)

```

returns:

```
BrunnerMunzelResult(statistic=inf, pvalue=nan)
```

### Reproducing Code Example

```python
x = np.array([1,2,3])
y = np.array([5,6,7,8,9])

scipy.stats.brunnermunzel(x,y)
```


### Error message

```shell
None.
```


### SciPy/NumPy/Python version information

1.8.0 1.21.5 sys.version_info(major=3, minor=8, micro=12, releaselevel='final', serial=0)

Hints from issue discussion:
I haven't investigated carefully, but it does not have to do with the sample sizes being unequal:
```python3
x1 = [1,2,3,4,5]
x2 = [5,6]
w, p_value = scipy.stats.brunnermunzel(x1, x2)  # OK
```
or simply that the total sample size is too small. It looks to me like it depends on the particular values. This should be assessed for correctness and, if it is correct, it might help to document it in the notes section. 

@sg-s, I'd suggest that you take a look at the definition of the `brunnermunzel` statistic and determining whether `inf` is indeed the appropriate value. This is not necessarily an implementation bug but an issue with the test statistic. (If the p-value is to be calculated with a normal or t approximation, we probably could return a p-value of `0` or `1`, whichever is appropriate. But the p-value of an exact permutation test would certainly not be `0` or `1`, so an argument might be made for NaN rather than suggesting an extreme p-value like this.) 

hi @mdhaber your counterexample is interesting. it's an issue that happens when you have small sizes and unequal sample sizes.

i should mention that it also throws this cryptic warning:

```
/Users/srinivas/opt/anaconda3/envs/longo/lib/python3.8/site-packages/scipy/stats/_stats_py.py:8055: RuntimeWarning:

divide by zero encountered in double_scalars

/Users/srinivas/opt/anaconda3/envs/longo/lib/python3.8/site-packages/scipy/stats/_stats_py.py:8061: RuntimeWarning:

invalid value encountered in double_scalars
```

meanwhile, this works:

```python
x = np.arange(0,12)
y = np.arange(5,20)

scipy.stats.brunnermunzel(x,y)
```

i suspect that when effect sizes are very large, (when statistic is inf) the pvalue should be 0, but it's being reported as `nan`

and so does this:


```python
x = np.arange(0,3)
y = x + 1

scipy.stats.brunnermunzel(x,y)
```

but this doesn't:


```python
x = np.arange(0,3)
y = x + 3

scipy.stats.brunnermunzel(x,y)
```





Well it seems to be happening when all elements of a group are greater than all elements of the other. I suspect `inf` might just be the value of the statistic you get in this case.
yup, i get that the value of the statistic can be (should be?) `inf`. the bug is that `p` is `nan`, when it should be 0
In your example, `df` is 0/0 in the t-approximation.
https://github.com/scipy/scipy/blob/b5d8bab88af61d61de09641243848df63380a67f/scipy/stats/_stats_py.py#L8057-L8062

If `df` is undefined, the distribution is undefined, and the p-value is undefined. 
One possibility would be to emit a warning suggesting the use of a normal approximation (which works) or permutation test.

If you submit a PR, I can review and merge!
@mdhaber thanks -- this is relatively straightforward, i can send you a PR

Commit message: BUG: adds warning to scipy.stats.brunnermunzel

scipy.stats.brunnermunzel returns a p-value of nan
in some cases, because the numerator and denominator
of the t-distribution were 0. It is not clear
why nan is returned, especially when it appears
that p should be 0 (when effect sizes are large).
Now, a warning is emitted when p is nan, and an
alternative is suggested

Fixes #15843

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index 577be6732ff7..9d0770b36de6 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -8232,7 +8232,16 @@ def brunnermunzel(x, y, alternative="two-sided", distribution="t",
         df_denom = np.power(nx * Sx, 2.0) / (nx - 1)
         df_denom += np.power(ny * Sy, 2.0) / (ny - 1)
         df = df_numer / df_denom
+
+        if (df_numer == 0) and (df_denom == 0):
+            warnings.warn("p-value cannot be estimated  "
+                          "because denominator and numerator of "
+                          "t-distribution are 0. Try using "
+                          "distribution='normal' ", RuntimeWarning)
+
         p = distributions.t.cdf(wbfn, df)
+
+
     elif distribution == "normal":
         p = distributions.norm.cdf(wbfn)
     else:


Code review comments:
diff_hunk: @@ -8232,7 +8232,16 @@ def brunnermunzel(x, y, alternative="two-sided", distribution="t",
         df_denom = np.power(nx * Sx, 2.0) / (nx - 1)
         df_denom += np.power(ny * Sy, 2.0) / (ny - 1)
         df = df_numer / df_denom
+
+        if (df_numer == 0) and (df_denom == 0):
+            warnings.warn("p-value cannot be estimated  "
+                          "because denominator and numerator of "
+                          "t-distribution are 0. Try using "
+                          "distribution='normal' ", RuntimeWarning)
path: scipy/stats/_stats_py.py
text: @user1:
```suggestion
        if (df_numer == 0) and (df_denom == 0):
            message = ("p-value cannot be estimated with `distribution='t' " 
                       "because degrees of freedom parameter is undefined "
                       "(0/0). Try using `distribution='normal'")
            warnings.warn(message, RuntimeWarning)
```
Slight rewording, because "denominator and numerator of t-distribution are 0" is less precise, and this saves a line.

diff_hunk: @@ -8232,7 +8232,16 @@ def brunnermunzel(x, y, alternative="two-sided", distribution="t",
         df_denom = np.power(nx * Sx, 2.0) / (nx - 1)
         df_denom += np.power(ny * Sy, 2.0) / (ny - 1)
         df = df_numer / df_denom
+
+        if (df_numer == 0) and (df_denom == 0):
+            warnings.warn("p-value cannot be estimated  "
+                          "because denominator and numerator of "
+                          "t-distribution are 0. Try using "
+                          "distribution='normal' ", RuntimeWarning)
+
         p = distributions.t.cdf(wbfn, df)
+
+
path: scipy/stats/_stats_py.py
text: @user1:
```suggestion

```
To fix `scipy/stats/_stats_py.py:8245:5: E303 too many blank lines (2)`