BUG: address the overflow bug of inverse sf of gamma when input proba…

…bility is very small

scipy/stats/_continuous_distns.py

@@ -3578,7 +3578,8 @@ def f_small_q(q, mu, k):
             """select the start values for the newton using the quantiles of
             the gamma (if right tail is required) or normal distribution"""
             if upper:
-                return gamma.isf(q, 1 / mu, scale=mu * mu)
+                mu3 = mu ** 3
+                return sc.gammainccinv(1 / mu3, q) / mu3
             return mu / (_norm_ppf(q) ** 2),
 
         k = 1.5 * mu

scipy/stats/tests/test_distributions.py

@@ -1894,14 +1894,12 @@ def test_ppf_sf(self):
 
         # tests if algorithm does not diverge for small probabilities.
         assert np.isclose(stats.invgauss.ppf(0.00013, mu=1 / 3) * 3,
-                          0.15039762631802803, atol=1e-08)
+                          0.15039762631802803, atol=1e-17)
         # test if it returns right tail values accurately
-        assert np.isclose(stats.invgauss.isf(1e-16, mu=1.05) / 0.7,
-                          98.47905825943425, atol=1e-08)
-        # test if it overflows for input 1e-17 or less
-        with pytest.warns(RuntimeWarning,
-            match="invalid value encountered in multiply"):
-             assert_equal(stats.invgauss.isf(1e-17, mu=1.05) / 0.7, np.nan)
+        assert np.isclose(stats.invgauss.isf(1e-17, mu=1.05) / 0.7,
+                          105.42079072389444, atol=1e-14)
+        assert np.isclose(stats.invgauss.isf(1e-50, mu=1.05) / 0.7,
+                          339.3560284841935, atol=1e-13)
         # test if correct  out is returned for boundary values
         with np.errstate(invalid='ignore'):
             # because probabilities must be in the interval [0, 1] a