MAINT: stats.rv_discrete: avoid infinite loops

Author: mdhaber

scipy/stats/_distn_infrastructure.py

@@ -3048,34 +3048,47 @@ def _drv2_ppfsingle(self, q, *args):  # Use basic bisection algorithm
     _a, _b = self._get_support(*args)
     b = _b
     a = _a
+
+    step = 10
     if isinf(b):            # Be sure ending point is > q
-        b = int(max(100*q, 10))
+        b = float(max(100*q, 10))
         while 1:
             if b >= _b:
                 qb = 1.0
                 break
             qb = self._cdf(b, *args)
             if (qb < q):
-                b += 10
+                b += step
+                step *= 2
             else:
                 break
     else:
         qb = 1.0
+
+    step = 10
     if isinf(a):    # be sure starting point < q
-        a = int(min(-100*q, -10))
+        a = float(min(-100*q, -10))
         while 1:
             if a <= _a:
                 qb = 0.0
                 break
             qa = self._cdf(a, *args)
             if (qa > q):
-                a -= 10
+                a -= step
+                step *= 2
             else:
                 break
     else:
         qa = self._cdf(a, *args)
 
-    while 1:
+    if np.isinf(a) or np.isinf(b):
+        message = "Arguments that bracket the requested quantile could not be found."
+        raise RuntimeError(message)
+
+    # maximum number of bisections within the normal numbers
+    # maxiter = int(np.log2(finfo.max) - np.log2(finfo.smallest_normal))
+    maxiter = 2046
+    for i in range(maxiter):
         if (qa == q):
             return a
         if (qb == q):

scipy/stats/tests/test_discrete_distns.py

@@ -1,6 +1,7 @@
 import pytest
 import itertools
 
+from scipy import stats
 from scipy.stats import (betabinom, betanbinom, hypergeom, nhypergeom,
                          bernoulli, boltzmann, skellam, zipf, zipfian, binom,
                          nbinom, nchypergeom_fisher, nchypergeom_wallenius,
@@ -646,3 +647,15 @@ def test_gh20692(self):
         pmf = dist.pmf(k)
         pmf_k_int32 = dist.pmf(k_int32)
         assert_equal(pmf, pmf_k_int32)
+
+
+def test_gh20048():
+    class test_dist_gen(stats.rv_discrete):
+        def _cdf(self, k):
+            return min(k / 100, 0.99)
+
+    test_dist = test_dist_gen(b=np.inf)
+
+    message = "Arguments that bracket..."
+    with pytest.raises(RuntimeError, match=message):
+        test_dist.ppf(0.999)

scipy/stats/tests/test_distributions.py

@@ -2180,13 +2180,6 @@ def test_mean_small_p(self):
         #   1.000000005
         assert_allclose(m, 1.000000005)
 
-    def test_gh20048(self):
-        # gh-20048 reported an infinite loop; see if we can reproduce it
-        p, q = 0.25739239425774363, 0.9999999999999999
-        res = stats.logser(p).ppf(q)
-        ref = 99
-        assert_equal(res, ref)
-
 
 class TestGumbel_r_l:
     @pytest.fixture(scope='function')