Repository: scipy/scipy
Pull request: ENH: `stats.ContinuousDistribution`: improve method resolution logic
#### Reference issue
Closes gh-2509

#### What does this implement/fix?
Previously, when the `tol` property of `ContinuousDistribution` was left unspecified, the CCDF/CDF could be calculated as the numerical complement of the CDF/CCDF. This logic was unusual because it ignored the *magnitude* of `tol`, and it was suboptimal because it did not protect against catastrophic cancellation in the default case (`tol` unspecified). 

This PR assesses the loss of precision when using the complement method to compute the CCDF/CDF and uses `quadrature` (in the appropriate tail) when the precision is worse than `tol`. It makes a similar adjustment to the logic of the inverse CCDF/CDF and 2-arg CDF.

There was similar logic in place for evaluating the log-CDF/CCDF as the logarithm of the CDF/CCDF. This is inaccurate when the CDF/CCDF underflows, so this PR checks for underflow and falls back to log-quadrature. It makes similar adjustments to the logic of `logexp` and 2-arg CDF.  

This PR also makes a few other small improvements/adjustments:
- The result dtype of 2-arg `logcdf` previously depended on the values of the arguments. Now, the result dtype is complex regardless of whether `y < x`.
- Fixes a copy-paste error in the `ContinuousDistribution.__sub__` method.
- Adds a private method for computing distribution moments by integrating the inverse CDF (rather than the PDF).

#### Additional information
Other obvious improvements to the default logic:

- ~~We could assess whether it is safe to compute a log-method as the logarithm of the non-log method depending on whether the non-log method over/under flows.~~ Done here.
- When the tolerance is quite strict, it might not be safe to compute a non-log method as the exponential of the log-method. We could detect this by checking `exp(logf(x)) - exp(logf(x) + np.nextafter(logf(x), np.inf)` (which is the absolute error corresponding to 1 ULP of error in `logf(x)`).
- There are also some situations in which it might be preferable to use numerical calculations rather than using log-complements or `logexpxmexpy`.
- For improvements when `locdf` and `logccdf` are not defined and the log CDF/CCDF probability is close to 0, we could apply the ideas from gh-21301 to the new infrastructure. (Note that the user can easily perform the main trick by choosing `method='complement'`; it's just not the default behavior.)

<details>
Regarding the safety of computing a non-log method as the exponential of the log-method:

```python3
import numpy as np
import matplotlib.pyplot as plt

x0 = np.logspace(-300, 300, 1000)
x1 = np.log(x0)
x2 = np.nextafter(x1, np.inf)
err = abs((np.exp(x2) - np.exp(x1))/np.exp(x2))
plt.semilogy(x1, err)
plt.xlabel('`x`')
plt.ylabel('Relative Error')
plt.title('Error in `exp(x)` corresponding with 1 ULP error in `x`')
plt.show()
```

![image](https://github.com/user-attachments/assets/a2167bea-8350-435c-9c0b-24c3815c5e4f)

The maximum relative error is ~1e-13. Exponentiation is only used if there is a formula for the log-method, so assuming the formula for the log-method is accurate, the exponential of the log-method will still be quite accurate.

</details>

Problem statement:
distributions cdf, sf evaluation in the wrong tail
If we calculate cdf or sf from the wrong tail, then we often run into numerical problems

for example :
gh-2238 numerical issues in the wrong tail, discrete, hypergeom
gh-2494 fix for truncnorm

two problems that I think still exist for some distributions
- missing calculations for _sf : when sf is defined as 1 - cdf, then we can have numerical problems in the right tail
- no automatic dispatch to the numerically more precise tail calculations: If we can recognize if the probabilities are calculated for one of the tails, then we could use either the _cdf or the _sf calculations to get more precise answers.

gh-2494 does a similar conditional calculation based on the support of the distribution.

I haven't checked any other example, but for example the generic method in rv_continuous is

```
    def _sf(self, x, *args):
        return 1.0-self._cdf(x,*args)
```

which means that we use in some cases integrate.quad on the wrong tail.

Hints from issue discussion:
Adding generic `_sf` shouldn't be too much work, recognizing the tail may be tricky.

The approach that would generalize what I did for the truncnorm fix would be as follows.  First, integrate from both ends, regardless of whether you want cdf or sf.  One of those values will be smaller, call it eps.  Depending on whether you want cdf or sf, return either eps or 1-eps as appropriate.  Ignore the larger of the two integrals, because it is a less precise evaluation of 1-eps.  I'm not sure if this method would work in general, but it would be kind of easy to implement.  It would also be slower because it would require evaluating more things.

I agree (reply to @rgommers )

Conditional calculations is tricky, also because of vectorization (need masks to dispatch)

I just realized that gh-2494 rules out the case when a or b are not scalars. There are already problems with vectorized solutions when the support is an argument, but we managed to resolve some of those IIRC.

@argriffing I think calculating both tails would be too expensive in many cases. When we need to calculate, for example, ppf generically, we need to solve for the root of a cdf equation, which requires that we need to evaluate the cdf very often. 

Calculating both, _sf and _cdf, would roughly double the time, and punish users who use the correct tail.

It _might_ be a solution where the cdf, sf calculations are cheap and numerical problems serious.

It could be better to be accurate but twice as slow by default, and then patch the distributions on a case-by-case basis to speed them up while still passing the unit tests?  Also some optimizations could possibly be applied to whole classes of distributions (e.g. symmetric distributions).

argriffing maybe the idea works if we do it sequentially

first do the standard calculation
then, for values close to 1 (eg. > 0.99 or 0.999) redo the calculation with the other tail

This would restrict the double calculations to only the extreme cases.

(It's still quite a bit of work to implement this)

> It could be better to be accurate but twice as slow by default

not if users find a way to do it correctly IMO, and it makes "standard" usage slow. 
We are using stats.distributions quite a lot in statsmodels, and I wouldn't want even more overhead for "correct" usage. (kind of, we don't use exotic distributions much)

For example, I _never_ do pval = 1 - distr.cdf(...) in any statistical test, (even if I just read it this morning in an article, and it took me 3 minutes of thinking.)

(aside: another use case: 
calling distr.rvs(..., size=1) repeatedly in a loop, spends most of the time in overhead. But that's a problem of usage not of design.)

In that case, maybe it could raise an exception?  So in the case of the reported bug, when they try to do `1-hypergeom.cdf(30,13397950,4363,12390)` they might see something like `StatsException('bro you are doing it wrong, maybe you want to use the survival function')` when scipy.stats sees that the cdf is near 1?  I guess that your sequential evaluation idea would also work in this case.

the problem is that we don't have any precision estimates

2.0581985005208026e-09  versus 1.228906e-11 in matlab in gh-2238 might in general still be the regular precision.
for example integrate.quad is used with default precision 1e-8 and the actual precision might fluctuate below that.

We can do something to improve specific cases, and the numerical precision in stats.distribution has increased a lot in the last few years (various fixes in scipy.stats, many improvements in scipy.special), but in general there is no guarantee that the 8th decimal is correct.
(And when I write a statistical test, I don't care much about the 8th decimal.)

We work our way slowly to improve precision, but I think raising precision errors everywhere is not really useful.

when the call to the distribution is hypergeom.cdf(30,13397950,4363,12390), then we don't know if the user actually wants 1-cdf or sf, and not just the cdf close to 1.   (That it's above 1 is a bug.)

hi,

i'm working on an implementation of the generalized gaussian distribution (this one: https://en.wikipedia.org/wiki/Generalized_normal_distribution#Version_1) and plan on contributing that once i have implemented all the parts.

while testing my implementation against its special cases (gaussian and laplacian), i noticed that the sf and isf of the laplacian aren't accurate in the upper tail (where they should be accurate). i haven't followed all of the discussion above, but it seemed this is the right place to report my finding.

for any symmetric distribution, it seems easy enough to define the sf in terms of the cdf as:

```
def _sf( self, x, pars ):
    return self._cdf( -x, pars )

def _isf( self, x, pars ):
    return -self._ppf( x, pars )
```

i did so in my implementation. i don't know enough about the implementation of the base class to decide whether it would make sense to define these automatically for all symmetric distributions. but it seems, at least for the laplace distribution, this should be fixed.

thanks!

Commit message: Update scipy/stats/tests/test_continuous.py

Patch to review:
diff --git a/scipy/stats/_distribution_infrastructure.py b/scipy/stats/_distribution_infrastructure.py
index 0f70397003ba..7b59a4646b0c 100644
--- a/scipy/stats/_distribution_infrastructure.py
+++ b/scipy/stats/_distribution_infrastructure.py
@@ -1103,7 +1103,6 @@ def wrapped(self, x, y, *args, **kwargs):
         if func_name in {'_cdf2', '_ccdf2'}:
             res = np.clip(res, 0., 1.)
         else:
-            res = res.real  # exp(res) > 0
             res = np.clip(res, None, 0.)  # exp(res) < 1
 
         # Transform the result to account for swapped argument order
@@ -1810,7 +1809,7 @@ def __add__(self, loc):
         return ShiftedScaledDistribution(self, loc=loc)
 
     def __sub__(self, loc):
-        return ShiftedScaledDistribution(self, loc=-self.loc)
+        return ShiftedScaledDistribution(self, loc=-loc)
 
     def __mul__(self, scale):
         return ShiftedScaledDistribution(self, scale=scale)
@@ -2095,8 +2094,8 @@ def logentropy(self, *, method=None):
     def _logentropy_dispatch(self, method=None, **params):
         if self._overrides('_logentropy_formula'):
             method = self._logentropy_formula
-        elif _isnull(self.tol) and self._overrides('_entropy_formula'):
-            method = self._logentropy_logexp
+        elif self._overrides('_entropy_formula'):
+            method = self._logentropy_logexp_safe
         else:
             method = self._logentropy_quadrature
         return method
@@ -2108,6 +2107,15 @@ def _logentropy_logexp(self, **params):
         res = np.log(self._entropy_dispatch(**params)+0j)
         return _log_real_standardize(res)
 
+    def _logentropy_logexp_safe(self, **params):
+        out = self._logentropy_logexp(**params)
+        mask = np.isinf(out.real)
+        if np.any(mask):
+            params_mask = {key:val[mask] for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._logentropy_quadrature(**params_mask)
+        return out[()]
+
     def _logentropy_quadrature(self, **params):
         def logintegrand(x, **params):
             logpdf = self._logpdf_dispatch(x, **params)
@@ -2299,7 +2307,8 @@ def logcdf(self, x, y=None, /, *, method=None):
 
     @_cdf2_input_validation
     def _logcdf2(self, x, y, *, method):
-        return self._logcdf2_dispatch(x, y, method=method, **self._parameters)
+        out = self._logcdf2_dispatch(x, y, method=method, **self._parameters)
+        return (out + 0j) if not np.issubdtype(out.dtype, np.complexfloating) else out
 
     @_dispatch
     def _logcdf2_dispatch(self, x, y, *, method=None, **params):
@@ -2310,9 +2319,9 @@ def _logcdf2_dispatch(self, x, y, *, method=None, **params):
         elif (self._overrides('_logcdf_formula')
               or self._overrides('_logccdf_formula')):
             method = self._logcdf2_subtraction
-        elif _isnull(self.tol) and (self._overrides('_cdf_formula')
-                                    or self._overrides('_ccdf_formula')):
-            method = self._logcdf2_logexp
+        elif (self._overrides('_cdf_formula')
+              or self._overrides('_ccdf_formula')):
+            method = self._logcdf2_logexp_safe
         else:
             method = self._logcdf2_quadrature
         return method
@@ -2335,13 +2344,23 @@ def _logcdf2_subtraction(self, x, y, **params):
         log_tail = np.logaddexp(logcdf_x, logccdf_y)[case_central]
         log_mass[case_central] = _log1mexp(log_tail)
         log_mass[flip_sign] += np.pi * 1j
-        return np.real_if_close(log_mass[()])
+        return log_mass[()]
 
     def _logcdf2_logexp(self, x, y, **params):
         expres = self._cdf2_dispatch(x, y, **params)
         expres = expres + 0j if np.any(x > y) else expres
         return np.log(expres)
 
+    def _logcdf2_logexp_safe(self, x, y, **params):
+        out = self._logcdf2_logexp(x, y, **params)
+        mask = np.isinf(out.real)
+        if np.any(mask):
+            params_mask = {key: np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._logcdf2_quadrature(x[mask], y[mask], **params_mask)
+        return out[()]
+
     def _logcdf2_quadrature(self, x, y, **params):
         logres = self._quadrature(self._logpdf_dispatch, limits=(x, y),
                                   log=True, params=params)
@@ -2355,10 +2374,10 @@ def _logcdf1(self, x, *, method=None):
     def _logcdf_dispatch(self, x, *, method=None, **params):
         if self._overrides('_logcdf_formula'):
             method = self._logcdf_formula
-        elif _isnull(self.tol) and self._overrides('_cdf_formula'):
-            method = self._logcdf_logexp
         elif self._overrides('_logccdf_formula'):
             method = self._logcdf_complement
+        elif self._overrides('_cdf_formula'):
+            method = self._logcdf_logexp_safe
         else:
             method = self._logcdf_quadrature
         return method
@@ -2366,11 +2385,21 @@ def _logcdf_dispatch(self, x, *, method=None, **params):
     def _logcdf_formula(self, x, **params):
         raise NotImplementedError(self._not_implemented)
 
+    def _logcdf_complement(self, x, **params):
+        return _log1mexp(self._logccdf_dispatch(x, **params))
+
     def _logcdf_logexp(self, x, **params):
         return np.log(self._cdf_dispatch(x, **params))
 
-    def _logcdf_complement(self, x, **params):
-        return _log1mexp(self._logccdf_dispatch(x, **params))
+    def _logcdf_logexp_safe(self, x, **params):
+        out = self._logcdf_logexp(x, **params)
+        mask = np.isinf(out)
+        if np.any(mask):
+            params_mask = {key:np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._logcdf_quadrature(x[mask], **params_mask)
+        return out[()]
 
     def _logcdf_quadrature(self, x, **params):
         a, _ = self._support(**params)
@@ -2395,9 +2424,8 @@ def _cdf2_dispatch(self, x, y, *, method=None, **params):
         elif (self._overrides('_logcdf_formula')
               or self._overrides('_logccdf_formula')):
             method = self._cdf2_logexp
-        elif _isnull(self.tol) and (self._overrides('_cdf_formula')
-                                    or self._overrides('_ccdf_formula')):
-            method = self._cdf2_subtraction
+        elif self._overrides('_cdf_formula') or self._overrides('_ccdf_formula'):
+            method = self._cdf2_subtraction_safe
         else:
             method = self._cdf2_quadrature
         return method
@@ -2416,8 +2444,33 @@ def _cdf2_subtraction(self, x, y, **params):
         cdf_y = self._cdf_dispatch(y, **params)
         ccdf_x = self._ccdf_dispatch(x, **params)
         ccdf_y = self._ccdf_dispatch(y, **params)
-        i = (cdf_x < 0.5) & (cdf_y < 0.5)
-        return np.where(i, cdf_y-cdf_x, ccdf_x-ccdf_y)
+        i = (ccdf_x < 0.5) & (ccdf_y < 0.5)
+        return np.where(i, ccdf_x-ccdf_y, cdf_y-cdf_x)
+
+    def _cdf2_subtraction_safe(self, x, y, **params):
+        cdf_x = self._cdf_dispatch(x, **params)
+        cdf_y = self._cdf_dispatch(y, **params)
+        ccdf_x = self._ccdf_dispatch(x, **params)
+        ccdf_y = self._ccdf_dispatch(y, **params)
+        i = (ccdf_x < 0.5) & (ccdf_y < 0.5)
+        out = np.where(i, ccdf_x-ccdf_y, cdf_y-cdf_x)
+
+        eps = np.finfo(self._dtype).eps
+        tol = self.tol if not _isnull(self.tol) else np.sqrt(eps)
+
+        cdf_max = np.maximum(cdf_x, cdf_y)
+        ccdf_max = np.maximum(ccdf_x, ccdf_y)
+        spacing = np.spacing(np.where(i, ccdf_max, cdf_max))
+
+        with np.errstate(divide='ignore'):
+            mask = tol < abs(spacing/out)
+
+        if np.any(mask):
+            params_mask = {key: np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._cdf2_quadrature(x[mask], y[mask], *params_mask)
+        return out[()]
 
     def _cdf2_quadrature(self, x, y, **params):
         return self._quadrature(self._pdf_dispatch, limits=(x, y), params=params)
@@ -2432,8 +2485,8 @@ def _cdf_dispatch(self, x, *, method=None, **params):
             method = self._cdf_formula
         elif self._overrides('_logcdf_formula'):
             method = self._cdf_logexp
-        elif _isnull(self.tol) and self._overrides('_ccdf_formula'):
-            method = self._cdf_complement
+        elif self._overrides('_ccdf_formula'):
+            method = self._cdf_complement_safe
         else:
             method = self._cdf_quadrature
         return method
@@ -2447,6 +2500,19 @@ def _cdf_logexp(self, x, **params):
     def _cdf_complement(self, x, **params):
         return 1 - self._ccdf_dispatch(x, **params)
 
+    def _cdf_complement_safe(self, x, **params):
+        out = self._cdf_complement(x, **params)
+        eps = np.finfo(self._dtype).eps
+        tol = self.tol if not _isnull(self.tol) else np.sqrt(eps)
+        with np.errstate(divide='ignore'):
+             mask = tol < eps/out
+        if np.any(mask):
+            params_mask = {key: np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._cdf_quadrature(x[mask], *params_mask)
+        return out[()]
+
     def _cdf_quadrature(self, x, **params):
         a, _ = self._support(**params)
         return self._quadrature(self._pdf_dispatch, limits=(a, x),
@@ -2488,10 +2554,10 @@ def _logccdf1(self, x, *, method=None):
     def _logccdf_dispatch(self, x, method=None, **params):
         if self._overrides('_logccdf_formula'):
             method = self._logccdf_formula
-        elif _isnull(self.tol) and self._overrides('_ccdf_formula'):
-            method = self._logccdf_logexp
         elif self._overrides('_logcdf_formula'):
             method = self._logccdf_complement
+        elif self._overrides('_ccdf_formula'):
+            method = self._logccdf_logexp_safe
         else:
             method = self._logccdf_quadrature
         return method
@@ -2499,11 +2565,21 @@ def _logccdf_dispatch(self, x, method=None, **params):
     def _logccdf_formula(self):
         raise NotImplementedError(self._not_implemented)
 
+    def _logccdf_complement(self, x, **params):
+        return _log1mexp(self._logcdf_dispatch(x, **params))
+
     def _logccdf_logexp(self, x, **params):
         return np.log(self._ccdf_dispatch(x, **params))
 
-    def _logccdf_complement(self, x, **params):
-        return _log1mexp(self._logcdf_dispatch(x, **params))
+    def _logccdf_logexp_safe(self, x, **params):
+        out = self._logccdf_logexp(x, **params)
+        mask = np.isinf(out)
+        if np.any(mask):
+            params_mask = {key: np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._logccdf_quadrature(x[mask], **params_mask)
+        return out[()]
 
     def _logccdf_quadrature(self, x, **params):
         _, b = self._support(**params)
@@ -2547,8 +2623,8 @@ def _ccdf_dispatch(self, x, method=None, **params):
             method = self._ccdf_formula
         elif self._overrides('_logccdf_formula'):
             method = self._ccdf_logexp
-        elif _isnull(self.tol) and self._overrides('_cdf_formula'):
-            method = self._ccdf_complement
+        elif self._overrides('_cdf_formula'):
+            method = self._ccdf_complement_safe
         else:
             method = self._ccdf_quadrature
         return method
@@ -2562,6 +2638,19 @@ def _ccdf_logexp(self, x, **params):
     def _ccdf_complement(self, x, **params):
         return 1 - self._cdf_dispatch(x, **params)
 
+    def _ccdf_complement_safe(self, x, **params):
+        out = self._ccdf_complement(x, **params)
+        eps = np.finfo(self._dtype).eps
+        tol = self.tol if not _isnull(self.tol) else np.sqrt(eps)
+        with np.errstate(divide='ignore'):
+            mask = tol < eps/out
+        if np.any(mask):
+            params_mask = {key: np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._ccdf_quadrature(x[mask], *params_mask)
+        return out[()]
+
     def _ccdf_quadrature(self, x, **params):
         _, b = self._support(**params)
         return self._quadrature(self._pdf_dispatch, limits=(x, b),
@@ -2600,8 +2689,8 @@ def icdf(self, p, /, *, method=None):
     def _icdf_dispatch(self, x, method=None, **params):
         if self._overrides('_icdf_formula'):
             method = self._icdf_formula
-        elif _isnull(self.tol) and self._overrides('_iccdf_formula'):
-            method = self._icdf_complement
+        elif self._overrides('_iccdf_formula'):
+            method = self._icdf_complement_safe
         else:
             method = self._icdf_inversion
         return method
@@ -2612,6 +2701,18 @@ def _icdf_formula(self, x, **params):
     def _icdf_complement(self, x, **params):
         return self._iccdf_dispatch(1 - x, **params)
 
+    def _icdf_complement_safe(self, x, **params):
+        out = self._icdf_complement(x, **params)
+        eps = np.finfo(self._dtype).eps
+        tol = self.tol if not _isnull(self.tol) else np.sqrt(eps)
+        mask = tol < eps/x
+        if np.any(mask):
+            params_mask = {key: np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._icdf_inversion(x[mask], *params_mask)
+        return out[()]
+
     def _icdf_inversion(self, x, **params):
         return self._solve_bounded(self._cdf_dispatch, x, params=params)
 
@@ -2646,8 +2747,8 @@ def iccdf(self, p, /, *, method=None):
     def _iccdf_dispatch(self, x, method=None, **params):
         if self._overrides('_iccdf_formula'):
             method = self._iccdf_formula
-        elif _isnull(self.tol) and self._overrides('_icdf_formula'):
-            method = self._iccdf_complement
+        elif self._overrides('_icdf_formula'):
+            method = self._iccdf_complement_safe
         else:
             method = self._iccdf_inversion
         return method
@@ -2658,6 +2759,18 @@ def _iccdf_formula(self, x, **params):
     def _iccdf_complement(self, x, **params):
         return self._icdf_dispatch(1 - x, **params)
 
+    def _iccdf_complement_safe(self, x, **params):
+        out = self._iccdf_complement(x, **params)
+        eps = np.finfo(self._dtype).eps
+        tol = self.tol if not _isnull(self.tol) else np.sqrt(eps)
+        mask = tol < eps/x
+        if np.any(mask):
+            params_mask = {key: np.broadcast_to(val, mask.shape)[mask]
+                           for key, val in params.items()}
+            out = np.asarray(out)
+            out[mask] = self._iccdf_inversion(x[mask], *params_mask)
+        return out[()]
+
     def _iccdf_inversion(self, x, **params):
         return self._solve_bounded(self._ccdf_dispatch, x, params=params)
 
@@ -2810,6 +2923,9 @@ def _moment_raw_dispatch(self, order, *, methods, **params):
         if moment is None and 'quadrature' in methods:
             moment = self._moment_integrate_pdf(order, center=self._zero, **params)
 
+        if moment is None and 'quadrature_icdf' in methods:
+            moment = self._moment_integrate_icdf(order, center=self._zero, **params)
+
         if moment is not None and self.cache_policy != _NO_CACHE:
             self._moment_raw_cache[order] = moment
 
@@ -2870,6 +2986,11 @@ def _moment_central_dispatch(self, order, *, methods, **params):
                                              methods=self._moment_methods)
             moment = self._moment_integrate_pdf(order, center=mean, **params)
 
+        if moment is None and 'quadrature_icdf' in methods:
+            mean = self._moment_raw_dispatch(self._one, **params,
+                                             methods=self._moment_methods)
+            moment = self._moment_integrate_icdf(order, center=mean, **params)
+
         if moment is not None and self.cache_policy != _NO_CACHE:
             self._moment_central_cache[order] = moment
 
@@ -2960,6 +3081,13 @@ def integrand(x, order, center, **params):
             return pdf*(x-center)**order
         return self._quadrature(integrand, args=(order, center), params=params)
 
+    def _moment_integrate_icdf(self, order, center, **params):
+        def integrand(x, order, center, **params):
+            x = self._icdf_dispatch(x, **params)
+            return (x-center)**order
+        return self._quadrature(integrand, limits=(0., 1.),
+                                args=(order, center), params=params)
+
     def _moment_transform_center(self, order, moment_as, a, b):
         a, b, *moment_as = np.broadcast_arrays(a, b, *moment_as)
         n = order
@@ -3683,7 +3811,7 @@ def __init__(self, components, *, weights=None):
         self._shape = np.broadcast_shapes(*(var._shape for var in components))
         self._dtype, self._components = dtype, components
         self._weights = np.full(n, 1/n, dtype=dtype) if weights is None else weights
-        self.validation_policy = None  # needed for
+        self.validation_policy = None
 
     @property
     def components(self):
diff --git a/scipy/stats/_probability_distribution.py b/scipy/stats/_probability_distribution.py
index 77fad53736e5..2e7e7cab839c 100644
--- a/scipy/stats/_probability_distribution.py
+++ b/scipy/stats/_probability_distribution.py
@@ -1285,6 +1285,11 @@ def logcdf(self, x, y, /, *, method):
         Similarly, the term "logarithmic difference" of :math:`w` and :math:`z`
         is used here to mean :math:`\log(\exp(w)-\exp(z))`.
 
+        If ``y < x``, the CDF is negative, and therefore the log-CCDF
+        is complex with imaginary part :math:`\pi`. For
+        consistency, the result of this function always has complex dtype
+        when `y` is provided, regardless of the value of the imaginary part.
+
         References
         ----------
         .. [1] Cumulative distribution function, *Wikipedia*,
@@ -1477,7 +1482,6 @@ def logccdf(self, x, y, /, *, method):
         is used here to mean the :math:`\log(\exp(w)+\exp(z))`, AKA
         :math:`\text{LogSumExp}(w, z)`.
 
-
         References
         ----------
         .. [1] Cumulative distribution function, *Wikipedia*,
@@ -1630,7 +1634,7 @@ def logentropy(self, *, method):
         Notes
         -----
         If the entropy of a distribution is negative, then the log-entropy
-        is complex with imaginary part divisible by :math:`\pi`. For
+        is complex with imaginary part :math:`\pi`. For
         consistency, the result of this function always has complex dtype,
         regardless of the value of the imaginary part.
 
diff --git a/scipy/stats/tests/test_continuous.py b/scipy/stats/tests/test_continuous.py
index 911bd3b892a5..9195dd6fe79e 100644
--- a/scipy/stats/tests/test_continuous.py
+++ b/scipy/stats/tests/test_continuous.py
@@ -240,6 +240,105 @@ def test_plot(self):
         ax = X.plot()
         assert ax == plt.gca()
 
+    @pytest.mark.parametrize('method_name', ['cdf', 'ccdf'])
+    def test_complement_safe(self, method_name):
+        X = stats.Normal()
+        X.tol = 1e-12
+        p = np.asarray([1e-4, 1e-3])
+        func = getattr(X, method_name)
+        ifunc = getattr(X, 'i'+method_name)
+        x = ifunc(p, method='formula')
+        p1 = func(x, method='complement_safe')
+        p2 = func(x, method='complement')
+        assert_equal(p1[1], p2[1])
+        assert p1[0] != p2[0]
+        assert_allclose(p1[0], p[0], rtol=X.tol)
+
+    @pytest.mark.parametrize('method_name', ['cdf', 'ccdf'])
+    def test_icomplement_safe(self, method_name):
+        X = stats.Normal()
+        X.tol = 1e-12
+        p = np.asarray([1e-4, 1e-3])
+        func = getattr(X, method_name)
+        ifunc = getattr(X, 'i'+method_name)
+        x1 = ifunc(p, method='complement_safe')
+        x2 = ifunc(p, method='complement')
+        assert_equal(x1[1], x2[1])
+        assert x1[0] != x2[0]
+        assert_allclose(func(x1[0]), p[0], rtol=X.tol)
+
+    def test_subtraction_safe(self):
+        X = stats.Normal()
+        X.tol = 1e-12
+
+        # Regular subtraction is fine in either tail (and of course, across tails)
+        x = [-11, -10, 10, 11]
+        y = [-10, -11, 11, 10]
+        p0 = X.cdf(x, y, method='quadrature')
+        p1 = X.cdf(x, y, method='subtraction_safe')
+        p2 = X.cdf(x, y, method='subtraction')
+        assert_equal(p2, p1)
+        assert_allclose(p1, p0, rtol=X.tol)
+
+        # Safe subtraction is needed in special cases
+        x = np.asarray([-1e-20, -1e-21, 1e-20, 1e-21, -1e-20])
+        y = np.asarray([-1e-21, -1e-20, 1e-21, 1e-20, 1e-20])
+        p0 = X.pdf(0)*(y-x)
+        p1 = X.cdf(x, y, method='subtraction_safe')
+        p2 = X.cdf(x, y, method='subtraction')
+        assert_equal(p2, 0)
+        assert_allclose(p1, p0, rtol=X.tol)
+
+    def test_logentropy_safe(self):
+        # simulate an `entropy` calculation over/underflowing with extreme parameters
+        class _Normal(stats.Normal):
+            def _entropy_formula(self, **params):
+                out = np.asarray(super()._entropy_formula(**params))
+                out[0] = 0
+                out[-1] = np.inf
+                return out
+
+        X = _Normal(sigma=[1, 2, 3])
+        with np.errstate(divide='ignore'):
+            res1 = X.logentropy(method='logexp_safe')
+            res2 = X.logentropy(method='logexp')
+        ref = X.logentropy(method='quadrature')
+        i_fl = [0, -1]  # first and last
+        assert np.isinf(res2[i_fl]).all()
+        assert res1[1] == res2[1]
+        assert res1[1] != ref[1]
+        assert_equal(res1[i_fl], ref[i_fl])
+
+    def test_logcdf2_safe(self):
+        # test what happens when 2-arg `cdf` underflows
+        X = stats.Normal(sigma=[1, 2, 3])
+        x = [-301, 1, 300]
+        y = [-300, 2, 301]
+        with np.errstate(divide='ignore'):
+            res1 = X.logcdf(x, y, method='logexp_safe')
+            res2 = X.logcdf(x, y, method='logexp')
+        ref = X.logcdf(x, y, method='quadrature')
+        i_fl = [0, -1]  # first and last
+        assert np.isinf(res2[i_fl]).all()
+        assert res1[1] == res2[1]
+        # quadrature happens to be perfectly accurate on some platforms
+        # assert res1[1] != ref[1]
+        assert_equal(res1[i_fl], ref[i_fl])
+
+    @pytest.mark.parametrize('method_name', ['logcdf', 'logccdf'])
+    def test_logexp_safe(self, method_name):
+        # test what happens when `cdf`/`ccdf` underflows
+        X = stats.Normal(sigma=2)
+        x = [-301, 1] if method_name == 'logcdf' else [301, 1]
+        func = getattr(X, method_name)
+        with np.errstate(divide='ignore'):
+            res1 = func(x, method='logexp_safe')
+            res2 = func(x, method='logexp')
+        ref = func(x, method='quadrature')
+        assert res1[0] == ref[0]
+        assert res1[0] != res2[0]
+        assert res1[1] == res2[1]
+        assert res1[1] != ref[1]
 
 def check_sample_shape_NaNs(dist, fname, sample_shape, result_shape, rng):
     full_shape = sample_shape + result_shape
@@ -369,7 +468,7 @@ def check_cdf2(dist, log, x, y, result_shape, methods):
         res = (np.exp(dist.logcdf(x, y, method=method)) if log
                else dist.cdf(x, y, method=method))
         np.testing.assert_allclose(res, ref, atol=1e-14)
-        if log and np.any(x > y) and ref.size:
+        if log:
             np.testing.assert_equal(res.dtype, (ref + 0j).dtype)
         else:
             np.testing.assert_equal(res.dtype, ref.dtype)
@@ -515,7 +614,10 @@ def has_formula(order, kind):
         assert ref.shape == result_shape
         check(i, 'raw','cache', ref, success=True)  # cached now
         check(i, 'raw', 'formula', ref, success=has_formula(i, 'raw'))
-        check(i, 'raw', 'general', ref, i == 0)
+        check(i, 'raw', 'general', ref, success=(i == 0))
+        if dist.__class__ == stats.Normal:
+            check(i, 'raw', 'quadrature_icdf', ref, success=True)
+
 
     # Clearing caches to better check their behavior
     dist.reset_cache()
@@ -542,6 +644,8 @@ def has_formula(order, kind):
         check(i, 'central', 'cache', ref, success=True)
         check(i, 'central', 'formula', ref, success=has_formula(i, 'central'))
         check(i, 'central', 'general', ref, success=i <= 1)
+        if dist.__class__ == stats.Normal:
+            check(i, 'central', 'quadrature_icdf', ref, success=True)
         check(i, 'central', 'transform', ref,
               success=has_formula(i, 'raw') or (i <= 1))
         if not has_formula(i, 'raw'):
@@ -562,7 +666,7 @@ def has_formula(order, kind):
         dist.moment(i, 'standardized')  # build up the cache
         check(i, 'central', 'normalize', ref)
 
-    ### Check Standard Moments ###
+    ### Check Standardized Moments ###
 
     var = dist.moment(2, 'central', method='quadrature')
     dist.reset_cache()
@@ -847,22 +951,6 @@ def test_tol(self):
         assert_allclose(res2, ref, rtol=X2.tol)
         assert abs(res2 - ref) > abs(res1 - ref)
 
-        # Test the tolerance logic in one dispatch method
-        # When tol is set, quadrature should be used -> correct entropy.
-        # When tol is not set, logexp should be used -> incorrect entropy.
-        wrong_entropy = 1.23456
-
-        class TestDist(ContinuousDistribution):
-            _variable = _RealParameter('x', domain=_RealDomain(endpoints=(0, 0.5)))
-            def _logpdf_formula(self, x, *args, **kwargs):
-                return np.full_like(x, np.log(2.))
-            def _entropy_formula(self, *args, **kwargs):
-                return wrong_entropy
-
-        X0 = _Uniform(a=0., b=0.5)
-        assert_allclose(TestDist(tol=1e-10).logentropy(), X0.logentropy())
-        assert_allclose(TestDist().logentropy(), np.log(wrong_entropy))
-
 
     def test_iv_policy(self):
         X = _Uniform(a=0, b=1)


Code review comments:
diff_hunk: @@ -240,6 +240,105 @@ def test_plot(self):
         ax = X.plot()
         assert ax == plt.gca()
 
+    @pytest.mark.parametrize('method_name', ['cdf', 'ccdf'])
+    def test_complement_safe(self, method_name):
+        X = stats.Normal()
+        X.tol = 1e-12
+        p = np.asarray([1e-4, 1e-3])
+        func = getattr(X, method_name)
+        ifunc = getattr(X, 'i'+method_name)
+        x = ifunc(p, method='formula')
+        p1 = func(x, method='complement_safe')
+        p2 = func(x, method='complement')
+        assert_equal(p1[1], p2[1])
+        assert p1[0] != p2[0]
+        assert_allclose(p1[0], p[0], rtol=X.tol)
+
+    @pytest.mark.parametrize('method_name', ['cdf', 'ccdf'])
+    def test_icomplement_safe(self, method_name):
+        X = stats.Normal()
+        X.tol = 1e-12
+        p = np.asarray([1e-4, 1e-3])
+        func = getattr(X, method_name)
+        ifunc = getattr(X, 'i'+method_name)
+        x1 = ifunc(p, method='complement_safe')
+        x2 = ifunc(p, method='complement')
+        assert_equal(x1[1], x2[1])
+        assert x1[0] != x2[0]
+        assert_allclose(func(x1[0]), p[0], rtol=X.tol)
+
+    def test_subtraction_safe(self):
+        X = stats.Normal()
+        X.tol = 1e-12
+
+        # Regular subtraction is fine in either tail (and of course, across tails)
+        x = [-11, -10, 10, 11]
+        y = [-10, -11, 11, 10]
+        p0 = X.cdf(x, y, method='quadrature')
+        p1 = X.cdf(x, y, method='subtraction_safe')
+        p2 = X.cdf(x, y, method='subtraction')
+        assert_equal(p2, p1)
+        assert_allclose(p1, p0, rtol=X.tol)
+
+        # Safe subtraction is needed in special cases
+        x = np.asarray([-1e-20, -1e-21, 1e-20, 1e-21, -1e-20])
+        y = np.asarray([-1e-21, -1e-20, 1e-21, 1e-20, 1e-20])
+        p0 = X.pdf(0)*(y-x)
+        p1 = X.cdf(x, y, method='subtraction_safe')
+        p2 = X.cdf(x, y, method='subtraction')
+        assert_equal(p2, 0)
+        assert_allclose(p1, p0, rtol=X.tol)
+
+    def test_logentropy_safe(self):
+        # simulate an `entropy` calculation over/underflowing with extreme parameters
+        class _Normal(stats.Normal):
+            def _entropy_formula(self, **params):
+                out = np.asarray(super()._entropy_formula(**params))
+                out[0] = 0
+                out[-1] = np.inf
+                return out
+
+        X = _Normal(sigma=[1, 2, 3])
+        with np.errstate(divide='ignore'):
+            res1 = X.logentropy(method='logexp_safe')
+            res2 = X.logentropy(method='logexp')
+        ref = X.logentropy(method='quadrature')
+        i_fl = [0, -1]  # first and last
+        assert np.isinf(res2[i_fl]).all()
+        assert res1[1] == res2[1]
+        assert res1[1] != ref[1]
path: scipy/stats/tests/test_continuous.py
text: @author:
```suggestion
        # quadrature happens to be perfectly accurate on some platforms
        # assert res1[1] != ref[1]
```