Repository: scipy/scipy
Pull request: Fix ttest permutations 
<!-- 
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 gh-16459

#### What does this implement/fix?
<!--Please explain your changes.-->
As stated in the issue, the p-values should not be zero.
If all permutations are used, the p-value remains the same as before. If only a few random permutations are used, the pvalue estimate is replaced by a biased estimate. This enables correct type-I error control for any value of the parameter `permutations`.
In practice, the estimated pvalues change as `pvalue_new = (p_value_old * permutations + 1) / (permutations + 1).


#### Additional information
<!--Any additional information you think is important.-->
I added some additional test that shows the issue.
I had to update the expected pvalues for some of the existent test according to above rule.

Problem statement:
BUG: stats.ttest_ind returns wrong p-values with permutations
### Describe your issue.

When running the two-sample t-test with permutations we do not make any assumptions about the distributions.
A user would reject the two-sample null hypothesis when obtaining a p-value smaller or equal to their significance level and we should **guarantee** that under the null hypothesis this does not happen with probability larger than the significance level.
Hence, the p-value should never be 0!

But with the current implementation it does.

Note that if a large number of permutations is used, this does not lead to a problem. But an inexperienced user might wrongly reject their null-hypothesis, because they used few permutations and received pvalue equal 0.


### Reproducing Code Example

```python
from scipy import stats
import numpy as np

rng = np.random.default_rng()

for i in range(200):
  rvs1 = stats.norm.rvs(loc=5, scale=10, size=500, random_state=rng)
  rvs2 = stats.norm.rvs(loc=5, scale=10, size=500, random_state=rng)
  print(stats.ttest_ind(rvs1, rvs2, permutations=1)[1])
```


### Error message

```shell
N/A
```


### SciPy/NumPy/Python version information

1.7.3 1.21.6 sys.version_info(major=3, minor=7, micro=13, releaselevel='final', serial=0)

Hints from issue discussion:
@mortonjt maybe you can have a look here. I think you committed the original code. 
This is correct in `stats.permutation_test`. References that support the need to add +1 in the numerator and denominator for randomized tests are [2] and [3] of the [documentation](https://docs.scipy.org/doc//scipy/reference/generated/scipy.stats.permutation_test.html).

I'll have to look back at why I didn't already fix this in the permutation t-test in particular, but we should.
@mdhaber 
Maybe the reason was that if the set number of permutations is larger than the number of different permutations, the method takes a deterministic approach. 
In this case, adding +1 is not necessary and should actually not be done!
Therefore, one additionally has to check whether this was the case.
Furthermore we have to adjust the hard-coded expected p-values for the tests.
I have done that for some (not all yet) in my PR https://github.com/scipy/scipy/pull/16460

Commit message: STY: stats.ttest_ind: fix PEP8 lines too long in tests

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index 5cb0cb8e83b7..b008f282111c 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -6692,7 +6692,7 @@ def _permutation_distribution_t(data, permutations, size_a, equal_var,
 
     t_stat = np.concatenate(t_stat, axis=0)
 
-    return t_stat, permutations
+    return t_stat, permutations, n_max
 
 
 def _calc_t_stat(a, b, equal_var, axis=-1):
@@ -6761,7 +6761,7 @@ def _permutation_ttest(a, b, permutations, axis=0, equal_var=True,
     mat = _broadcast_concatenate((a, b), axis=axis)
     mat = np.moveaxis(mat, axis, -1)
 
-    t_stat, permutations = _permutation_distribution_t(
+    t_stat, permutations, n_max = _permutation_distribution_t(
         mat, permutations, size_a=na, equal_var=equal_var,
         random_state=random_state)
 
@@ -6771,7 +6771,10 @@ def _permutation_ttest(a, b, permutations, axis=0, equal_var=True,
 
     # Calculate the p-values
     cmps = compare[alternative](t_stat, t_stat_observed)
-    pvalues = cmps.sum(axis=0) / permutations
+    if n_max > permutations:
+        pvalues = (cmps.sum(axis=0) + 1) / (permutations + 1)
+    else:
+        pvalues = cmps.sum(axis=0) / permutations
 
     # nans propagate naturally in statistic calculation, but need to be
     # propagated manually into pvalues
diff --git a/scipy/stats/tests/test_resampling.py b/scipy/stats/tests/test_resampling.py
index 2e3d053f4405..465d2d21e638 100644
--- a/scipy/stats/tests/test_resampling.py
+++ b/scipy/stats/tests/test_resampling.py
@@ -994,13 +994,7 @@ def statistic(x, y, axis):
                                 random_state=0)
 
         assert_allclose(res1.statistic, res2.statistic, rtol=self.rtol)
-
-        if permutations == 30:
-            # Even one-sided p-value is defined differently in ttest_ind for
-            # randomized tests. See permutation_test references [2] and [3].
-            assert_allclose((res1.pvalue*30+1)/31, res2.pvalue, rtol=self.rtol)
-        else:
-            assert_allclose(res1.pvalue, res2.pvalue, rtol=self.rtol)
+        assert_allclose(res1.pvalue, res2.pvalue, rtol=self.rtol)
 
     # -- Independent (Unpaired) Sample Tests -- #
 
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 8961c478c3df..ca8e9ad4efb1 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -4385,9 +4385,10 @@ class Test_ttest_ind_permutations():
                           size=500).reshape(100, 5).T
     rvs2 = stats.norm.rvs(loc=8, scale=20, size=100)  # type: ignore
 
-    p_d = [0, 0.676]  # desired pvalues
-    p_d_gen = [0, 0.672]  # desired pvalues for Generator seed
-    p_d_big = [0.993, 0.685, 0.84, 0.955, 0.255]
+    p_d = [1/1001, (676+1)/1001]  # desired pvalues
+    p_d_gen = [1/1001, (672 + 1)/1001]  # desired pvalues for Generator seed
+    p_d_big = [(993+1)/1001, (685+1)/1001, (840+1)/1001,
+               (955+1)/1001, (255+1)/1001]
 
     params = [
         (a, b, {"axis": 1}, p_d),                     # basic test
@@ -4397,7 +4398,7 @@ class Test_ttest_ind_permutations():
         # different seeds
         (a, b, {'random_state': 0, "axis": 1}, p_d),
         (a, b, {'random_state': np.random.RandomState(0), "axis": 1}, p_d),
-        (a2, b2, {'equal_var': True}, 0),  # equal variances
+        (a2, b2, {'equal_var': True}, 1/1001),  # equal variances
         (rvs1, rvs2, {'axis': -1, 'random_state': 0}, p_d_big),  # bigger test
         (a3, b3, {}, 1/3)  # exact test
         ]
@@ -4519,8 +4520,10 @@ def test_ttest_ind_randperm_alternative(self):
 
         # For random permutations, the chance of ties between the observed
         # test statistic and the population is small, so:
-        assert_equal(res_g_ab.pvalue + res_l_ab.pvalue, 1)
-        assert_equal(res_g_ba.pvalue + res_l_ba.pvalue, 1)
+        assert_equal(res_g_ab.pvalue + res_l_ab.pvalue,
+                     1 + 1/(options_p['permutations'] + 1))
+        assert_equal(res_g_ba.pvalue + res_l_ba.pvalue,
+                     1 + 1/(options_p['permutations'] + 1))
 
     @pytest.mark.slow()
     def test_ttest_ind_randperm_alternative2(self):
@@ -4600,6 +4603,15 @@ def test_ttest_ind_permutation_check_inputs(self):
             stats.ttest_ind(self.a, self.b, permutations=1,
                             random_state='hello')
 
+    def test_ttest_ind_permutation_check_p_values(self):
+        # p-values should never be exactly zero
+        N = 10
+        a = np.random.rand(N, 20)
+        b = np.random.rand(N, 20)
+        p_values = stats.ttest_ind(a, b, permutations=1).pvalue
+        print(0.0 not in p_values)
+        assert 0.0 not in p_values
+
 
 class Test_ttest_ind_common:
     # for tests that are performed on variations of the t-test such as


Code review comments:
diff_hunk: @@ -6771,7 +6771,10 @@ def _permutation_ttest(a, b, permutations, axis=0, equal_var=True,
 
     # Calculate the p-values
     cmps = compare[alternative](t_stat, t_stat_observed)
-    pvalues = cmps.sum(axis=0) / permutations
+    if n_max > permutations:
+        pvalues = (cmps.sum(axis=0) + 1) / (permutations + 1)
+    else:
+        pvalues = cmps.sum(axis=0) / permutations
path: scipy/stats/_stats_py.py
text: @user1:
Considered:
```suggestion
    adjustment = 1 if n_max > permutations else 0
    pvalues = (cmps.sum(axis=0) + adjustment) / (permutations + adjustment)
```
but feel free to leave it alone, as many would prefer the explicit `if`s.

diff_hunk: @@ -6771,7 +6771,10 @@ def _permutation_ttest(a, b, permutations, axis=0, equal_var=True,
 
     # Calculate the p-values
     cmps = compare[alternative](t_stat, t_stat_observed)
-    pvalues = cmps.sum(axis=0) / permutations
+    if n_max > permutations:
path: scipy/stats/_stats_py.py
text: @user1:
```suggestion
    # Randomized test p-value calculation should use biased estimate; see e.g.
    # https://www.degruyter.com/document/doi/10.2202/1544-6115.1585/
    if n_max > permutations:
```