Repository: scipy/scipy
Pull request: Bugfix for somersd where an integer overflow could occur
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes scipy/scipy/issues/18132

#### What does this implement/fix?
<!--Please explain your changes.-->
Bugfix for somersd where an integer overflow could occur for cases that should be computable (Issue #18132).

#### Additional information
<!--Any additional information you think is important.-->
Modified the somersd function and added unit tests where the result is compared to an independent calculation using sklearn.


Problem statement:
BUG: invalid output and behavior of scipy.stats.somersd
### Describe your issue.

The Somers' D statistic is a measure of ordinal association between two possibly dependent random variables X and Y and is implemented in scipy.stats as the function somersd(x,...). 

This implementation has some unexpected and incorrect behavior, at least for large input samples (see minimal reproducing example below):
1.. There is a RuntimeWarning
2. Somers' D is outside of the range  [-1, +1], which is the only interval it shoudl give outputs
3. The p-value is not always defined, even when the Somers' D statistic is calculated, which it should be

### Reproducing Code Example

```python
from scipy.stats import somersd
import random

# generate lists of random classes 1-5
classes = [1, 2, 3, 4, 5]
n_samples = 10**6
x = random.choices(classes, k=n_samples)
y = random.choices(classes, k=n_samples)

# Calculate ordinal association with Somers' D:
somersd(x, y)
```


### Error message

```shell
RuntimeWarning: invalid value encountered in double_scalars
  Z = (PA - QA)/(4*(S))**0.5
Out[1]: 
SomersDResult(statistic=-3.6429990341561904, pvalue=nan, table=array([[39960, 39980, 39995, 39800, 40122],
       [40011, 40053, 40119, 39874, 40259],
       [40068, 39859, 39834, 39612, 40176],
       [40288, 39964, 39905, 40075, 40190],
       [39881, 39828, 39970, 40171, 40006]]))

#---- Comments and further experiments ----#
We would expect statistic to be around 0 and definitely not outside of the interval [-1,1] as it happens here. Also the pvalue should not be nan.

Setting n_samples=10**4 or n_samples=10**5, one does _sometimes_ get a RuntimeWarning and a pvalue that is nan, but a correct statistic.

Setting n_samples=10**3 one get correct results for the statistic an pvalue
```


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

```shell
1.10.0 1.23.5 sys.version_info(major=3, minor=8, micro=16, releaselevel='final', serial=0)
lapack_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\include']
lapack_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\include']
blas_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\include']
blas_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/denni/.conda/envs/inv\\Library\\include']
```

Hints from issue discussion:
I'll take a look. Could be an integer overflow.

Do you ever get incorrect output without the RuntimeWarning or a RuntimeWarning with perfectly correct output? Otherwise, I would not consider the RuntimeWarning behavior itself to be a bug. It is signaling that there is a problem and the output cannot be trusted, which is valuable.
It does look like an integer overflow.
Simply converting `table` in this line:
https://github.com/scipy/scipy/blob/c1ed5ece8ffbf05356a22a8106affcd11bd3aee0/scipy/stats/_hypotests.py#L874

with `.astype(float)` gives, for your example:
```
SomersDResult(statistic=0.00020258604036181154, pvalue=0.800093700827777, table=array([[40046, 39902, 40027, 40291, 39595],
       [39547, 40182, 40140, 39961, 39926],
       [40095, 40097, 40065, 39964, 39998],
       [39986, 39887, 40300, 39955, 39886],
       [39963, 40182, 39835, 40157, 40013]]))
```
with no warnings. 

You mentioned cases in which you get "a RuntimeWarning and a pvalue that is nan, but a correct statistic". Are you certain that the statistic is correct, or is it just reasonable-looking? (I can see how either could happen; just checking.) Can you provide other problematic cases (e.g. fix your random seed so that it is reproducible) that you think might indicate a different problem?

The simplest fix is probably to just always treat the table with a floating point type. The slightly better solution (in that it would preserve the last few decimal places of accuracy for slightly large tables/counts) would be to predict when floats will be needed and convert the type only when necessary.
> I'll take a look. Could be an integer overflow.
> 
> Do you ever get incorrect output without the RuntimeWarning or a RuntimeWarning with perfectly correct output? Otherwise, I would not consider the RuntimeWarning behavior itself to be a bug. It is signaling that there is a problem and the output cannot be trusted, which is valuable.

Hi Matt, thanks for taking a look at this! 
No, I never got the incorrect output without the warning; I have checked the output independently in several cases by comparing with other packages (sk-learn's roc-auc function) and they agree.
Ok. Can you compute reference values for this case?

```python3
from scipy.stats import somersd
import random

# generate lists of random classes 1-5
classes = [1, 2, 3, 4, 5]
np.random.seed(6272161)

n_samples = 10**6
x = random.choices(classes, k=n_samples)
y = random.choices(classes, k=n_samples)

# Calculate ordinal association with Somers' D:
res = somersd(x, y)
```
I just fixed the random seed to `6272161`.
With the change described above (https://github.com/scipy/scipy/issues/18132#issuecomment-1464961026), I get `statistic=0.0002617496287395554, pvalue=0.7434399693577043`
> It does look like an integer overflow. Simply converting `table` in this line:
> 
> https://github.com/scipy/scipy/blob/c1ed5ece8ffbf05356a22a8106affcd11bd3aee0/scipy/stats/_hypotests.py#L874
> 
> with `.astype(float)` gives, for your example:
> 
> ```
> SomersDResult(statistic=0.00020258604036181154, pvalue=0.800093700827777, table=array([[40046, 39902, 40027, 40291, 39595],
>        [39547, 40182, 40140, 39961, 39926],
>        [40095, 40097, 40065, 39964, 39998],
>        [39986, 39887, 40300, 39955, 39886],
>        [39963, 40182, 39835, 40157, 40013]]))
> ```
> 
> with no warnings.
> 
> You mentioned cases in which you get "a RuntimeWarning and a pvalue that is nan, but a correct statistic". Are you certain that the statistic is correct, or is it just reasonable-looking? (I can see how either could happen; just checking.) Can you provide other problematic cases (e.g. fix your random seed so that it is reproducible) that you think might indicate a different problem?
> 
> The simplest fix is probably to just always treat the table with a floating point type. The slightly better solution (in that it would preserve the last few decimal places of accuracy for slightly large tables/counts) would be to predict when floats will be needed and convert the type only when necessary.

In this particular example I only checked that statistic has a value that is around what one would expect (since they are random, there should not be any ordinal association, and thus the statistic around 0.0).
Your investigation looks promising!
> Ok. Can you compute reference values for this case?
> 
> ```python
> from scipy.stats import somersd
> import random
> 
> # generate lists of random classes 1-5
> classes = [1, 2, 3, 4, 5]
> np.random.seed(6272161)
> 
> n_samples = 10**6
> x = random.choices(classes, k=n_samples)
> y = random.choices(classes, k=n_samples)
> 
> # Calculate ordinal association with Somers' D:
> res = somersd(x, y)
> ```
> 
> I just fixed the random seed to `6272161`. With the change described above ([#18132 (comment)](https://github.com/scipy/scipy/issues/18132#issuecomment-1464961026)), I get `statistic=0.0002617496287395554, pvalue=0.7434399693577043`

Unfortunately I cannot verify this one directly, as I would have to go to e.g. R, and then I am not sure that I would be able to generate the exact same random numbers. But the output looks very reasonable.

I have below here another example where I get a RuntimeWarning and wrong outputs from somersd, which I can check by other means:

```python
from scipy.stats import somersd
import random
import numpy as np
from sklearn import metrics

# generate lists of random classes 1-2 (binary)
classes = [1, 2]
n_samples = 10**6
random.seed(6272161)
x = random.choices(classes, k=n_samples)
y = random.choices(classes, k=n_samples)

# Calculate ordinal association with sklearn ROC-AUC method and Somers' D:
val1_sklearn = 2*metrics.roc_auc_score(x, y)-1  # Convert to the Gini coefficient (Gini = (AUC*2)-1) (which should be equal to the result of Somers' D)
val2_scipy = somersd(x, y).statistic
```
Here one gets the correct result from sklearn and the wrong result from somersd (and a nan pvalue):
  sklearn.metrics.roc_auc_score: -0.001528138777036947
  scipy.stats.somersd: -0.4283415175162947

These numbers you can compre with your proposed fix
Can you set the random seed (e.g. `np.random.seed(6272161)`) so that I can reproduce this data?

> Can you set the random seed (e.g. `np.random.seed(6272161)`) so that I can reproduce this data?

Yes, sorry I forgot about the seed - I have updated the post above (use random.seed and not np.random.seed, since I use the random module)
Ok, with the change, SciPy agrees to ~machine precision. Would you like to submit a PR? If so, I can review and merge. I can submit the PR if you prefer, though.
> Ok, with the change, SciPy agrees to ~machine precision. Would you like to submit a PR? If so, I can review and merge. I can submit the PR if you prefer, though.

Perfect, I will submit a PR.

BTW, there is a closely related overflow problem in scipy.stats.kendalltau; I will create a seperate issue for that one and tag you

Commit message: Bugfix for somersd where an integer overflow could occur for cases that should be computable (Issue #18132). Modified the somersd function and added unit tests.

Patch to review:
diff --git a/scipy/stats/_hypotests.py b/scipy/stats/_hypotests.py
index da98054a4fff..c68c41a2d6c9 100644
--- a/scipy/stats/_hypotests.py
+++ b/scipy/stats/_hypotests.py
@@ -871,7 +871,8 @@ def somersd(x, y=None, alternative='two-sided'):
         table = x
     else:
         raise ValueError("x must be either a 1D or 2D array")
-    d, p = _somers_d(table, alternative)
+    # The table type is converted to a float to avoid an integer overflows
+    d, p = _somers_d(table.astype(float), alternative)
 
     # add alias for consistency with other correlation functions
     res = SomersDResult(d, p, table)
diff --git a/scipy/stats/tests/test_hypotests.py b/scipy/stats/tests/test_hypotests.py
index 50834c8cac03..b2bf29d4204d 100644
--- a/scipy/stats/tests/test_hypotests.py
+++ b/scipy/stats/tests/test_hypotests.py
@@ -1,6 +1,7 @@
 from itertools import product
 
 import numpy as np
+import random
 import functools
 import pytest
 from numpy.testing import (assert_, assert_equal, assert_allclose,
@@ -994,6 +995,29 @@ def test_somersd_perfect_correlation(self, positive_correlation):
         assert res.statistic == expected_statistic
         assert res.pvalue == (0 if positive_correlation else 1)
 
+    def test_somersd_large_inputs(self):
+        # Test that large inputs where potential overflows could occur give
+        # the expected output. This is tested in the case of binary inputs.
+
+        # generate lists of random classes 1-2 (binary)
+        classes = [1, 2]
+        n_samples = 10 ** 6
+        random.seed(6272161)
+        x = random.choices(classes, k=n_samples)
+        y = random.choices(classes, k=n_samples)
+
+        # get value to compare with: sklearn output
+        # from sklearn import metrics
+        # val_auc_sklearn = metrics.roc_auc_score(x, y)
+        # # convert to the Gini coefficient (Gini = (AUC*2)-1)
+        # val_sklearn = 2 * val_auc_sklearn - 1
+        val_sklearn = -0.001528138777036947
+
+        # calculate the Somers' D statistic, which should be equal to the result
+        # of val_sklearn until approximately machine precision
+        val_scipy = stats.somersd(x, y).statistic
+        assert_allclose(val_sklearn, val_scipy, atol=1e-15)
+
 
 class TestBarnardExact:
     """Some tests to show that barnard_exact() works correctly."""


Code review comments:
diff_hunk: @@ -871,7 +871,8 @@ def somersd(x, y=None, alternative='two-sided'):
         table = x
     else:
         raise ValueError("x must be either a 1D or 2D array")
-    d, p = _somers_d(table, alternative)
+    # The table type is converted to a float to avoid an integer overflows
path: scipy/stats/_hypotests.py
text: @user1:
```suggestion
    # The table type is converted to a float to avoid an integer overflow
```