MAINT: stats.Mixture: fix inverse functions when mean is undefined (scipy#22337)

* MAINT: stats.Mixture: fix inverse when mean is undefined

* Apply suggestions from code review

* Update scipy/stats/_distribution_infrastructure.py

Author: mdhaber

scipy/stats/_distribution_infrastructure.py

@@ -1212,6 +1212,25 @@ def _logexpxmexpy(x, y):
     return res
 
 
+def _guess_bracket(xmin, xmax):
+    a = np.full_like(xmin, -1.0)
+    b = np.ones_like(xmax)
+
+    i = np.isfinite(xmin) & np.isfinite(xmax)
+    a[i] = xmin[i]
+    b[i] = xmax[i]
+
+    i = np.isfinite(xmin) & ~np.isfinite(xmax)
+    a[i] = xmin[i]
+    b[i] = xmin[i] + 1
+
+    i = np.isfinite(xmax) & ~np.isfinite(xmin)
+    a[i] = xmax[i] - 1
+    b[i] = xmax[i]
+
+    return a, b
+
+
 def _log_real_standardize(x):
     """Standardizes the (complex) logarithm of a real number.
 
@@ -2011,27 +2030,13 @@ def f2(x, _p, **kwargs):  # named `_p` to avoid conflict with shape `p`
         shape = xmin.shape
         xmin, xmax = np.atleast_1d(xmin, xmax)
 
-        a = -np.ones_like(xmin)
-        b = np.ones_like(xmax)
-
-        i = np.isfinite(xmin) & np.isfinite(xmax)
-        a[i] = xmin[i]
-        b[i] = xmax[i]
-
-        i = np.isfinite(xmin) & ~np.isfinite(xmax)
-        a[i] = xmin[i]
-        b[i] = xmin[i] + 1
-
-        i = np.isfinite(xmax) & ~np.isfinite(xmin)
-        a[i] = xmax[i] - 1
-        b[i] = xmax[i]
-
+        xl0, xr0 = _guess_bracket(xmin, xmax)
         xmin = xmin.reshape(shape)
         xmax = xmax.reshape(shape)
-        a = a.reshape(shape)
-        b = b.reshape(shape)
+        xl0 = xl0.reshape(shape)
+        xr0 = xr0.reshape(shape)
 
-        res = _bracket_root(f3, xl0=a, xr0=b, xmin=xmin, xmax=xmax, args=args)
+        res = _bracket_root(f3, xl0=xl0, xr0=xr0, xmin=xmin, xmax=xmax, args=args)
         # For now, we ignore the status, but I want to use the bracket width
         # as an error estimate - see question 5 at the top.
         xrtol = None if _isnull(self.tol) else self.tol
@@ -4636,7 +4641,8 @@ def _invert(self, fun, p):
         xmin, xmax = self.support()
         fun = getattr(self, fun)
         f = lambda x, p: fun(x) - p  # noqa: E731 is silly
-        res = _bracket_root(f, xl0=self.mean(), xmin=xmin, xmax=xmax, args=(p,))
+        xl0, xr0 = _guess_bracket(xmin, xmax)
+        res = _bracket_root(f, xl0=xl0, xr0=xr0, xmin=xmin, xmax=xmax, args=(p,))
         return _chandrupatla(f, a=res.xl, b=res.xr, args=(p,)).x
 
     def icdf(self, p, /, *, method=None):

scipy/stats/tests/test_continuous.py

@@ -1863,3 +1863,17 @@ def test_properties(self):
         assert X.components[0] == components[0]
         X.weights[0] = weights[1]
         assert X.weights[0] == weights[0]
+
+    def test_inverse(self):
+        # Originally, inverse relied on the mean to start the bracket search.
+        # This didn't work for distributions with non-finite mean. Check that
+        # this is resolved.
+        rng = np.random.default_rng(24358934657854237863456)
+        Cauchy = stats.make_distribution(stats.cauchy)
+        X0 = Cauchy()
+        X = stats.Mixture([X0, X0])
+        p = rng.random(size=10)
+        np.testing.assert_allclose(X.icdf(p), X0.icdf(p))
+        np.testing.assert_allclose(X.iccdf(p), X0.iccdf(p))
+        np.testing.assert_allclose(X.ilogcdf(p), X0.ilogcdf(p))
+        np.testing.assert_allclose(X.ilogccdf(p), X0.ilogccdf(p))