Repository: scipy/scipy
Pull request: BUG/ENH: more stable recursion for 2-sample ks test exact p-values
Fixes: #12999

So this follows @mdhaber's in #12999 request to leave the code structure as is.
However, in the test sent by @sigallev (which I adapted into a test), this is about 30x slower on my laptop. I think we can recover quite a bit of that by moving into cython/pythran/ whatever (I used [numba before](https://github.com/TorchDrift/TorchDrift/blob/bda529a7a26c1a9b14d8d878a75122ece22b95d6/torchdrift/detectors/ks.py#L25), but note that I'm computing with < instead of <=).

I also have a test variant that reproduces the classic algorithm using Python's `fractions`, but I didn't include this here and instead stuck to the "should be in the same ballpark as asymp" heuristic of the issue.

I think it might be worth integrating improving speed instead of keeping the structure as close to the old one as possible.

Problem statement:
Bug in scipy.stats.ks_2samp for two-sided auto and exact modes on a few thousand samples (when n1!=n2) 
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

When you run have a few thousand samples for each group, and different a number of samples than auto goes to exact.
One sided - gets an exception and shifts to 'asymp'
Two sided - doesn't get an exception but may well finish with prob = 0 or 1, which is not the real probability.
                   the "finally" statement whould change to 'asymp' mode, if it was prob < 0 or >1, but not at 0 and 1.
                   I think this should change to prob <=0 or >=1.

This causes the code snippet below to work well for 1000, 1010 samples worse for 2000,2010 and fail completely for 3000,3010 and above, until 10,000 when the 'asymp' mode kicks in no matter what. 

#### Reproducing code example:

print("look at the p_values going to 1, for 3000 and above when in 'auto' mode. This is clearly not the right probability.")
for x in range(1000, 12000,1000):
    print("For %d vs %d auto" % (x, x))
    print(scipy.stats.ks_2samp([random.normalvariate(100,10) for i in range(x)], [random.normalvariate(95,10) for i in range(x)]))
    print("For %d vs %d asymptotic (not exact)" % (x, x))
    print(scipy.stats.ks_2samp([random.normalvariate(100, 10) for i in range(x)],
                               [random.normalvariate(95, 10) for i in range(x)], mode='asymp'))
    print()

#### Scipy/Numpy/Python version information:

1.4.1 1.18.1 sys.version_info(major=3, minor=7, micro=4, releaselevel='final', serial=0)

Hints from issue discussion:
There were some bugs of `ks_2samp` in SciPy 1.4, e.g. #11184, that were fixed in 1.5. Could you check if the problem still exists if you upgrade SciPy?
The results actually changed, *and are now even weirder*.

My code is very simple:
import random
import scipy
import scipy.stats

print(scipy.__version__)
for x in range(1000, 12000, 1000):
    vals1 = [random.normalvariate(100, 10) for i in range(x)]
    vals2 = [random.normalvariate(95, 10) for i in range(x+10)]
    print("\n%d" % x)
    print(scipy.stats.ks_2samp(vals1, vals2))
    print(scipy.stats.ks_2samp(vals1, vals2, mode='asymp'))
    # these two p-values should be in line with each other
print()

And look how it behaves. There are some major discrepancies between the two
values, but only for some of the xs (for some number of samples, not for
others):
1.5.3

1000
KstestResult(statistic=0.20679207920792078, pvalue=2.6391393024758554e-19)
KstestResult(statistic=0.2067920792079208, pvalue=2.639139302475824e-19)

2000
KstestResult(statistic=0.17853980099502487, pvalue=2.0522704781112646e-28)
KstestResult(statistic=0.17853980099502487, pvalue=2.0522704781112646e-28)

3000
KstestResult(statistic=0.22652380952380952, pvalue=3.2773195020948264e-68)
KstestResult(statistic=0.22652380952380957, pvalue=3.277319502094573e-68)

4000
KstestResult(statistic=0.22134351620947632, pvalue=2.220446049250313e-15)
KstestResult(statistic=0.2213435162094763, pvalue=1.254472004892521e-86)

5000
KstestResult(statistic=0.20126586826347306, pvalue=2.5567455151927657e-89)
KstestResult(statistic=0.20126586826347304, pvalue=2.5567455151929115e-89)

6000
KstestResult(statistic=0.2133230726566833, pvalue=4.107825191113079e-15)
KstestResult(statistic=0.21332307265668327, pvalue=2.2822895371211086e-120)

7000
KstestResult(statistic=0.18743896474424293, pvalue=3.4202940876631045e-108)
KstestResult(statistic=0.18743896474424288, pvalue=3.420294087663611e-108)

8000
KstestResult(statistic=0.19601404494382021, pvalue=2.5827118221855017e-12)
KstestResult(statistic=0.1960140449438202, pvalue=3.389282533241279e-135)

9000
KstestResult(statistic=0.2069876680231841, pvalue=0.006901878771671144)
KstestResult(statistic=0.20698766802318413, pvalue=1.2018998948969596e-169)

10000
KstestResult(statistic=0.20650299700299696, pvalue=1.5122150846105846e-187)
KstestResult(statistic=0.20650299700299696, pvalue=1.5122150846105846e-187)

11000
KstestResult(statistic=0.20663826273635538, pvalue=1.6137073644802917e-206)
KstestResult(statistic=0.20663826273635538, pvalue=1.6137073644802917e-206)


*Sigal Peled-Leviatan*


On Mon, Nov 2, 2020 at 10:19 PM Christoph Baumgarten <
notifications@github.com> wrote:

> There were some bugs of ks_2samp in SciPy 1.4, e.g. #11184
> <https://github.com/scipy/scipy/issues/11184>, that were fixed in 1.5.
> Could you check if the problem still exists if you upgrade SciPy?
>
> —
> You are receiving this because you authored the thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/scipy/scipy/issues/12999#issuecomment-720701878>, or
> unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AL3RAFGLLBDM4YTPCOU4WH3SN4H47ANCNFSM4S66VXEA>
> .
>

For the last two cases, it is expected that you get the same result since `auto` switches to `asymp` if one of the samples is larger than 10000. 

In the other cases, I am not sure how accurate the values are but for practical purposes, it looks like the pvalue is essentially zero as expected, i.e. the test detects that the samples have a different mean.

If `mode == 'exact'`, the p-value looks indeed a bit large in the case `x == 9000`, so potentially the switch to `asymp` should happen earlier? See comment in #11184:

> In general, the exact two-sided computation is best suited to situations in which one of the sample sizes is at most a few thousand, or the statistic d is "small", as the computation is O(N1 N2 d).

in your case, d is very large and so are N1 and N2.

@pvanmulbregt any view on this?
I know that atx10,000 it changes to asymp, but what's interesting is before
that.
Sometimes it computes arround 10**-15 and decides its good, and other times
arround -10**-15, realises this is less than 0 so does not make sense and
moves to asymp.

This suggests to me that any pvalue in this range is not trust worthy, and
should be replaced by the asymp calculation.

You can also see the problem by noting that 4000 points from the
distribution should be more significant than 1000 or 200, and is not. You
must agree this is a problem....

Sigal Peled-Leviatan

בתאריך יום ג׳, 3 בנוב׳ 2020, 21:37, מאת Christoph Baumgarten ‏<
notifications@github.com>:

> For the last two cases, it is expected that you get the same result since
> auto switches to asymp if one of the samples is larger than 10000.
>
> In the other cases, I am not sure how accurate the values are but for
> practical purposes, it looks like the pvalue is essentially zero as
> expected, i.e. the test detects that the samples have a different mean.
>
> If mode == 'exact', the p-value looks indeed a bit large, so potentially
> the switch to asymp should happen earlier? See comment in #11184
> <https://github.com/scipy/scipy/issues/11184>:
>
> In general, the exact two-sided computation is best suited to situations
> in which one of the sample sizes is at most a few thousand, or the
> statistic d is "small", as the computation is O(N1 N2 d).
>
> in your case, d is very large and so are N1 and N2.
>
> @pvanmulbregt <https://github.com/pvanmulbregt> any view on this?
>
> —
> You are receiving this because you authored the thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/scipy/scipy/issues/12999#issuecomment-721335057>, or
> unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AL3RAFAK553APB6CZWJRFWTSOBLY5ANCNFSM4S66VXEA>
> .
>

Indeed, that is linked to my statement above that the switch to `asymp` might need to happen earlier.

Note that `exact` does not mean the value is precise up to a certain number of digits. it just means the exact distribution of the test statistic is used. For any practical purpose, the numbers above should be small enough to reject the assumption that the distribution of the two samples is the same :)
Not necessarily, when you are dealing with multiple hypothesis testing
(GWAS, for example).
And I agree that moving to asymp earlier would have solved this, but this
needs some research with a bigger set of distributions to decide when to
move.
In general if a random partial set of samples gives a lower pvalue then
something is very wrong with an estimation.

Good luck with the code, we rely on your keeping up the great work.

Sigal Peled-Leviatan

בתאריך יום ד׳, 4 בנוב׳ 2020, 20:47, מאת Christoph Baumgarten ‏<
notifications@github.com>:

> Indeed, that is linked to my statement above that the switch to asymp
> might need to happen earlier.
>
> Note that exact does not mean the value is precise up to a certain number
> of digits. it just means the exact distribution of the test statistic is
> used. For any practical purpose, the numbers above should be small enough
> to reject the assumption that the distribution of the two samples is the
> same :)
>
> —
> You are receiving this because you authored the thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/scipy/scipy/issues/12999#issuecomment-721908316>, or
> unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AL3RAFDBWWQABUCNVR5MGNTSOGOUTANCNFSM4S66VXEA>
> .
>

The computation in stats._compute_prob_inside_method() returns the complementary probability, and this number may be very close to 1 if the samples do not not come for the same distribution. For large-sized mismatched samples, the accuracy of the computed probability p may not be sufficient for use when needing 1-p.  Certainly it won't be if the actual, probability is less than 1e-16, as in the examples above.
gh-13168 looks for insufficient accuracy and backs off to the asymptotic calculation in that situation.

Cross posting from the PR: I wrote up a more stable recursion for the inner method / 2-sided 2-sample in https://arxiv.org/abs/2102.08037

The basic trick that you can do instead of doing 1-A(m,n)/Binom(m+n,n) at the end, you can do the entire recursion on values of this form 1-A(i,j)/Binom(i+j,j). It seems to be about 4x slower than the current implementation, which I think it is due to the extra compute, but then both the current implementation and the stable recursion can be moved to numba / other ways to compile to speed them up. I can also make a PR of it if it helps.

Commenting to bump this issue, 

Here's another example of when things go wrong:
```
import numpy as np
print(np.__version__)
import scipy
print(scipy.__version__)
from scipy.stats import ks_2samp

np.random.seed(1)
x = np.random.normal(10, 1, 5_000_000)
y = np.random.normal(12, 1, 5_000)
print(ks_2samp(x, y, mode='auto', alternative = 'two-sided'))
print(ks_2samp(x, y, mode='exact', alternative = 'two-sided'))

```
Gives:
```
1.18.5
1.5.2
KstestResult(statistic=0.6859724, pvalue=0.0)
KstestResult(statistic=0.6859724, pvalue=1.0)

```
Commenting to further bump.  I was using ks_2samp to determine whether or not a set of experimental measurements were different.  For reference, here is a histogram of the first two measurement sets which have sizes 5971 and 6000, respectively:

![image](https://user-images.githubusercontent.com/36451772/115933279-689fcf80-a443-11eb-9afc-c445c33b577f.png)

Running `kstest(sample1, sample2)` returns the following:
```KstestResult(statistic=0.8563704013844694, pvalue=1.0)```

If I switch `alternative='greater'`, I get 
```KstestResult(statistic=0.8563704013844694, pvalue=0.0)```
which does not surprise me since these are very large and very different distributions

Switching to `mode='asymp'` makes the `pvalue=0.0`.

I also ran @srcoulombe 's code above and got the same result with a newer version of scipy
```python
import numpy as np
print(np.__version__)
import scipy
print(scipy.__version__)
from scipy.stats import ks_2samp

np.random.seed(1)
x = np.random.normal(10, 1, 5_000_000)
y = np.random.normal(12, 1, 5_000)
print(ks_2samp(x, y, mode='auto', alternative = 'two-sided'))
print(ks_2samp(x, y, mode='exact', alternative = 'two-sided'))
```
```
1.19.2
1.6.2
KstestResult(statistic=0.6859724, pvalue=0.0)
KstestResult(statistic=0.6859724, pvalue=1.0)
```

Commit message: more stable recursion for 2-sample ks test

Fixes: #12999

Patch to review:
diff --git a/scipy/stats/_hypotests_pythran.py b/scipy/stats/_hypotests_pythran.py
index f965bbe60937..d83b497738e9 100644
--- a/scipy/stats/_hypotests_pythran.py
+++ b/scipy/stats/_hypotests_pythran.py
@@ -1,3 +1,5 @@
+import numpy as np
+
 #pythran export _Aij(float[:,:], int, int)
 #pythran export _Aij(int[:,:], int, int)
 def _Aij(A, i, j):
@@ -47,3 +49,100 @@ def _a_ij_Aij_Dij2(A):
         for j in range(n):
             count += A[i, j]*(_Aij(A, i, j) - _Dij(A, i, j))**2
     return count
+
+
+#pythran export _compute_outer_prob_inside_method(int64, int64, int64, int64)
+def _compute_outer_prob_inside_method(m, n, g, h):
+    """
+    Count the proportion of paths that do not stay strictly inside two
+    diagonal lines.
+
+    Parameters
+    ----------
+    m : integer
+        m > 0
+    n : integer
+        n > 0
+    g : integer
+        g is greatest common divisor of m and n
+    h : integer
+        0 <= h <= lcm(m,n)
+
+    Returns
+    -------
+    p : float
+        The proportion of paths that do not stay inside the two lines.
+
+    The classical algorithm counts the integer lattice paths from (0, 0)
+    to (m, n) which satisfy |x/m - y/n| < h / lcm(m, n).
+    The paths make steps of size +1 in either positive x or positive y
+    directions.
+    We are, however, interested in 1 - proportion to computs p-values,
+    so we change the recursion to compute 1 - p directly while staying
+    within the "inside method" a described by Hodges.
+
+    We generally follow Hodges' treatment of Drion/Gnedenko/Korolyuk.
+    Hodges, J.L. Jr.,
+    "The Significance Probability of the Smirnov Two-Sample Test,"
+    Arkiv fiur Matematik, 3, No. 43 (1958), 469-86.
+
+    For the recursion for 1-p see
+    Viehmann, T.: "Numerically more stable computation of the p-values
+    for the two-sample Kolmogorov-Smirnov test," arXiv: 2102.08037
+
+    """
+    # Probability is symmetrical in m, n.  Computation below uses m >= n.
+    if m < n:
+        m, n = n, m
+    mg = m // g
+    ng = n // g
+
+    # Count the integer lattice paths from (0, 0) to (m, n) which satisfy
+    # |nx/g - my/g| < h.
+    # Compute matrix A such that:
+    #  A(x, 0) = A(0, y) = 1
+    #  A(x, y) = A(x, y-1) + A(x-1, y), for x,y>=1, except that
+    #  A(x, y) = 0 if |x/m - y/n|>= h
+    # Probability is A(m, n)/binom(m+n, n)
+    # Optimizations exist for m==n, m==n*p.
+    # Only need to preserve a single column of A, and only a
+    # sliding window of it.
+    # minj keeps track of the slide.
+    minj, maxj = 0, min(int(np.ceil(h / mg)), n + 1)
+    curlen = maxj - minj
+    # Make a vector long enough to hold maximum window needed.
+    lenA = min(2 * maxj + 2, n + 1)
+    # This is an integer calculation, but the entries are essentially
+    # binomial coefficients, hence grow quickly.
+    # Scaling after each column is computed avoids dividing by a
+    # large binomial coefficient at the end, but is not sufficient to avoid
+    # the large dyanamic range which appears during the calculation.
+    # Instead we rescale based on the magnitude of the right most term in
+    # the column and keep track of an exponent separately and apply
+    # it at the end of the calculation.  Similarly when multiplying by
+    # the binomial coefficint
+    dtype = np.float64
+    A = np.ones(lenA, dtype=dtype)
+    # Initialize the first column
+    A[minj:maxj] = 0.0
+    for i in range(1, m + 1):
+        # Generate the next column.
+        # First calculate the sliding window
+        lastminj, lastlen = minj, curlen
+        minj = max(int(np.floor((ng * i - h) / mg)) + 1, 0)
+        minj = min(minj, n)
+        maxj = min(int(np.ceil((ng * i + h) / mg)), n + 1)
+        if maxj <= minj:
+            return 1.0
+        # Now fill in the values. We cannot use cumsum, unfortunately.
+        val = 0.0 if minj == 0 else 1.0
+        for jj in range(maxj - minj):
+            j = jj + minj
+            val = (A[jj + minj - lastminj] * i + val * j) / (i + j)
+            A[jj] = val
+        curlen = maxj - minj
+        if lastlen > curlen:
+            # Set some carried-over elements to 1
+            A[maxj - minj:maxj - minj + (lastlen - curlen)] = 1
+
+    return A[maxj - minj - 1]
diff --git a/scipy/stats/stats.py b/scipy/stats/stats.py
index 8eadb5e4549e..c5f07a071651 100644
--- a/scipy/stats/stats.py
+++ b/scipy/stats/stats.py
@@ -48,6 +48,7 @@
                      _local_correlations)
 from dataclasses import make_dataclass
 from ._hypotests import _all_partitions
+from ._hypotests_pythran import _compute_outer_prob_inside_method
 from ._axis_nan_policy import (_axis_nan_policy_factory,
                                _broadcast_concatenate)
 
@@ -7085,109 +7086,6 @@ def ks_1samp(x, cdf, args=(), alternative='two-sided', mode='auto'):
 Ks_2sampResult = KstestResult
 
 
-def _compute_prob_inside_method(m, n, g, h):
-    """
-    Count the proportion of paths that stay strictly inside two diagonal lines.
-
-    Parameters
-    ----------
-    m : integer
-        m > 0
-    n : integer
-        n > 0
-    g : integer
-        g is greatest common divisor of m and n
-    h : integer
-        0 <= h <= lcm(m,n)
-
-    Returns
-    -------
-    p : float
-        The proportion of paths that stay inside the two lines.
-
-    Count the integer lattice paths from (0, 0) to (m, n) which satisfy
-    |x/m - y/n| < h / lcm(m, n).
-    The paths make steps of size +1 in either positive x or positive y
-    directions.
-
-    We generally follow Hodges' treatment of Drion/Gnedenko/Korolyuk.
-    Hodges, J.L. Jr.,
-    "The Significance Probability of the Smirnov Two-Sample Test,"
-    Arkiv fiur Matematik, 3, No. 43 (1958), 469-86.
-
-    """
-    # Probability is symmetrical in m, n.  Computation below uses m >= n.
-    if m < n:
-        m, n = n, m
-    mg = m // g
-    ng = n // g
-
-    # Count the integer lattice paths from (0, 0) to (m, n) which satisfy
-    # |nx/g - my/g| < h.
-    # Compute matrix A such that:
-    #  A(x, 0) = A(0, y) = 1
-    #  A(x, y) = A(x, y-1) + A(x-1, y), for x,y>=1, except that
-    #  A(x, y) = 0 if |x/m - y/n|>= h
-    # Probability is A(m, n)/binom(m+n, n)
-    # Optimizations exist for m==n, m==n*p.
-    # Only need to preserve a single column of A, and only a
-    # sliding window of it.
-    # minj keeps track of the slide.
-    minj, maxj = 0, min(int(np.ceil(h / mg)), n + 1)
-    curlen = maxj - minj
-    # Make a vector long enough to hold maximum window needed.
-    lenA = min(2 * maxj + 2, n + 1)
-    # This is an integer calculation, but the entries are essentially
-    # binomial coefficients, hence grow quickly.
-    # Scaling after each column is computed avoids dividing by a
-    # large binomial coefficient at the end, but is not sufficient to avoid
-    # the large dyanamic range which appears during the calculation.
-    # Instead we rescale based on the magnitude of the right most term in
-    # the column and keep track of an exponent separately and apply
-    # it at the end of the calculation.  Similarly when multiplying by
-    # the binomial coefficint
-    dtype = np.float64
-    A = np.zeros(lenA, dtype=dtype)
-    # Initialize the first column
-    A[minj:maxj] = 1
-    expnt = 0
-    for i in range(1, m + 1):
-        # Generate the next column.
-        # First calculate the sliding window
-        lastminj, lastlen = minj, curlen
-        minj = max(int(np.floor((ng * i - h) / mg)) + 1, 0)
-        minj = min(minj, n)
-        maxj = min(int(np.ceil((ng * i + h) / mg)), n + 1)
-        if maxj <= minj:
-            return 0
-        # Now fill in the values
-        A[0:maxj - minj] = np.cumsum(A[minj - lastminj:maxj - lastminj])
-        curlen = maxj - minj
-        if lastlen > curlen:
-            # Set some carried-over elements to 0
-            A[maxj - minj:maxj - minj + (lastlen - curlen)] = 0
-        # Rescale if the right most value is over 2**900
-        val = A[maxj - minj - 1]
-        _, valexpt = math.frexp(val)
-        if valexpt > 900:
-            # Scaling to bring down to about 2**800 appears
-            # sufficient for sizes under 10000.
-            valexpt -= 800
-            A = np.ldexp(A, -valexpt)
-            expnt += valexpt
-
-    val = A[maxj - minj - 1]
-    # Now divide by the binomial (m+n)!/m!/n!
-    for i in range(1, n + 1):
-        val = (val * i) / (m + i)
-        _, valexpt = math.frexp(val)
-        if valexpt < -128:
-            val = np.ldexp(val, -valexpt)
-            expnt += valexpt
-    # Finally scale if needed.
-    return np.ldexp(val, expnt)
-
-
 def _compute_prob_outside_square(n, h):
     """
     Compute the proportion of paths that pass outside the two diagonal lines.
@@ -7333,7 +7231,7 @@ def _attempt_exact_2kssamp(n1, n2, g, d, alternative):
             if n1 == n2:
                 prob = _compute_prob_outside_square(n1, h)
             else:
-                prob = 1 - _compute_prob_inside_method(n1, n2, g, h)
+                prob = _compute_outer_prob_inside_method(n1, n2, g, h)
         else:
             if n1 == n2:
                 # prob = binom(2n, n-h) / binom(2n, n)
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 9cfd2a1439c1..b9ed682e86d5 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -16,7 +16,8 @@
                            assert_almost_equal, assert_array_almost_equal,
                            assert_array_equal, assert_approx_equal,
                            assert_allclose, assert_warns, suppress_warnings,
-                           assert_string_equal)
+                           assert_string_equal,
+                           assert_array_less)
 import pytest
 from pytest import raises as assert_raises
 import numpy.ma.testutils as mat
@@ -3885,6 +3886,17 @@ def test_gh11184_bigger(self):
         self._testOne(x, y, 'greater', 0.10597913208679133, 2.7947433906389253e-41, mode='asymp')
         self._testOne(x, y, 'less', 0.09658002199780022, 2.7947433906389253e-41, mode='asymp')
 
+    @pytest.mark.slow
+    def test_gh12999(self):
+        for x in range(1000, 12000, 1000):
+            vals1 = np.random.normal(size=(x))
+            vals2 = np.random.normal(size=(x + 10), loc=0.5)
+            exact = stats.ks_2samp(vals1, vals2, mode='exact').pvalue
+            asymp = stats.ks_2samp(vals1, vals2, mode='asymp').pvalue
+            # these two p-values should be in line with each other
+            assert_array_less(exact, 3 * asymp)
+            assert_array_less(asymp, 3 * exact)
+
     @pytest.mark.slow
     def testLargeBoth(self):
         # 10000, 11000
@@ -3912,9 +3924,12 @@ def testNamedAttributes(self):
     @pytest.mark.slow
     def test_some_code_paths(self):
         # Check that some code paths are executed
-        from scipy.stats.stats import _count_paths_outside_method, _compute_prob_inside_method
+        from scipy.stats.stats import (
+            _count_paths_outside_method,
+            _compute_outer_prob_inside_method
+        )
 
-        _compute_prob_inside_method(1, 1, 1, 1)
+        _compute_outer_prob_inside_method(1, 1, 1, 1)
         _count_paths_outside_method(1000, 1, 1, 1001)
 
         assert_raises(FloatingPointError, _count_paths_outside_method, 1100, 1099, 1, 1)


Code review comments:
diff_hunk: @@ -3885,6 +3886,17 @@ def test_gh11184_bigger(self):
         self._testOne(x, y, 'greater', 0.10597913208679133, 2.7947433906389253e-41, mode='asymp')
         self._testOne(x, y, 'less', 0.09658002199780022, 2.7947433906389253e-41, mode='asymp')
 
+    @pytest.mark.slow
path: scipy/stats/tests/test_stats.py
text: @user1:
This is quite fast on my machine. Is this marked slow only because Pythran is still optional and it's slow in pure-Python mode?

@author:
So it is rather slow on my notebook, but it is quite fast for my larger machine. I think the main reason I added it is because of the CI timeouts I saw. Can you advise whether you want to keep the slow tag? I'm quite happy to remove it if you don't need it.

@user1:
Ah okay, a large difference is a little surprising, since it's single-threaded code. Let me time it properly on my machine.

@user1:
1.2 sec on my machine. Probably 1.5-2x slower in CI, but I think that's okay. So I'd say it's fine to remove `@user2` unless it's way more than that in your notebook.

@user3:
@user1 What should we do about this one? It's the slowest test on a few builds, e.g. on Windows Python310-64bit-full recently:
```
============================ slowest 10 durations =============================
275.92s call     stats/tests/test_stats.py::TestKSTwoSamples::test_gh12999
```

_Update: marked as xslow in gh-15234, but perhaps this deserves a closer look, as it's unusual for there to be a huge difference on different machines._

@user1:
Yes, it may indeed warrant a closer look. That test has been problematic before, +1 for moving to `xslow`.

diff_hunk: @@ -47,3 +49,100 @@ def _a_ij_Aij_Dij2(A):
         for j in range(n):
             count += A[i, j]*(_Aij(A, i, j) - _Dij(A, i, j))**2
     return count
+
+
+#pythran export _compute_outer_prob_inside_method(int64, int64, int64, int64)
+def _compute_outer_prob_inside_method(m, n, g, h):
+    """
+    Count the proportion of paths that do not stay strictly inside two
+    diagonal lines.
+
+    Parameters
+    ----------
+    m : integer
+        m > 0
+    n : integer
+        n > 0
+    g : integer
+        g is greatest common divisor of m and n
+    h : integer
+        0 <= h <= lcm(m,n)
+
+    Returns
+    -------
+    p : float
+        The proportion of paths that do not stay inside the two lines.
+
+    The classical algorithm counts the integer lattice paths from (0, 0)
+    to (m, n) which satisfy |x/m - y/n| < h / lcm(m, n).
+    The paths make steps of size +1 in either positive x or positive y
+    directions.
+    We are, however, interested in 1 - proportion to computs p-values,
+    so we change the recursion to compute 1 - p directly while staying
+    within the "inside method" a described by Hodges.
+
+    We generally follow Hodges' treatment of Drion/Gnedenko/Korolyuk.
+    Hodges, J.L. Jr.,
+    "The Significance Probability of the Smirnov Two-Sample Test,"
+    Arkiv fiur Matematik, 3, No. 43 (1958), 469-86.
+
+    For the recursion for 1-p see
+    Viehmann, T.: "Numerically more stable computation of the p-values
+    for the two-sample Kolmogorov-Smirnov test," arXiv: 2102.08037
+
+    """
+    # Probability is symmetrical in m, n.  Computation below uses m >= n.
+    if m < n:
+        m, n = n, m
+    mg = m // g
+    ng = n // g
+
+    # Count the integer lattice paths from (0, 0) to (m, n) which satisfy
+    # |nx/g - my/g| < h.
+    # Compute matrix A such that:
+    #  A(x, 0) = A(0, y) = 1
+    #  A(x, y) = A(x, y-1) + A(x-1, y), for x,y>=1, except that
+    #  A(x, y) = 0 if |x/m - y/n|>= h
+    # Probability is A(m, n)/binom(m+n, n)
+    # Optimizations exist for m==n, m==n*p.
+    # Only need to preserve a single column of A, and only a
+    # sliding window of it.
+    # minj keeps track of the slide.
+    minj, maxj = 0, min(int(np.ceil(h / mg)), n + 1)
+    curlen = maxj - minj
+    # Make a vector long enough to hold maximum window needed.
+    lenA = min(2 * maxj + 2, n + 1)
+    # This is an integer calculation, but the entries are essentially
+    # binomial coefficients, hence grow quickly.
+    # Scaling after each column is computed avoids dividing by a
+    # large binomial coefficient at the end, but is not sufficient to avoid
+    # the large dyanamic range which appears during the calculation.
+    # Instead we rescale based on the magnitude of the right most term in
+    # the column and keep track of an exponent separately and apply
+    # it at the end of the calculation.  Similarly when multiplying by
+    # the binomial coefficint
path: scipy/stats/_hypotests_pythran.py
text: @user1:
minor: typo in "coefficient" - may as well fix it

diff_hunk: @@ -47,3 +49,100 @@ def _a_ij_Aij_Dij2(A):
         for j in range(n):
             count += A[i, j]*(_Aij(A, i, j) - _Dij(A, i, j))**2
     return count
+
+
+#pythran export _compute_outer_prob_inside_method(int64, int64, int64, int64)
+def _compute_outer_prob_inside_method(m, n, g, h):
+    """
+    Count the proportion of paths that do not stay strictly inside two
+    diagonal lines.
+
+    Parameters
+    ----------
+    m : integer
+        m > 0
+    n : integer
+        n > 0
+    g : integer
+        g is greatest common divisor of m and n
+    h : integer
+        0 <= h <= lcm(m,n)
+
+    Returns
+    -------
+    p : float
+        The proportion of paths that do not stay inside the two lines.
+
+    The classical algorithm counts the integer lattice paths from (0, 0)
+    to (m, n) which satisfy |x/m - y/n| < h / lcm(m, n).
+    The paths make steps of size +1 in either positive x or positive y
+    directions.
+    We are, however, interested in 1 - proportion to computs p-values,
path: scipy/stats/_hypotests_pythran.py
text: @user1:
minor: typo in "compute"