Repository: scipy/scipy
Pull request: ENH: stats.rv_continuous.expect: split interval to improve reliability
#### Reference issue
Closes gh-8928

#### What does this implement/fix?
gh-8928 reported many examples where `rv_continuous.expect` gave inaccurate results. Since it relies on `quad`, we cannot expect `expect` to be perfect for all intervals and integrands, but it can be improved based on knowledge of where the distribution's probability mass is non-negligble. This PR attempts to make `expect` more reliable (albeit slower) by splitting the interval at the tails and evaluating the integral for the left tail, body, and right tail separately.

#### Additional information
I will add tests from the original issue. 

Hopefully, this addresses the `stats`-specific part of the issue. I think that is all we can really aim for here: take advantage of the stats-specific part of the problem to improve the results. Beyond that, this is an `integrate` issue.

I skipped CI in the initial commit because I need to check locally which tests break due to removing use of `np.errstate` from `expect`. We can silience the integrator warnings in the tests, but users should see when the integrator knows it's having trouble.

Perhaps instead of splitting the integral explicitly, we can use the `points` argument of `quad`. I'll try that.

Problem statement:
BUG: scipy.stats.norm wrong expected value of function when loc is large
The expected value of a normal distribution should be the center of the distribution, but scipy gives wrong results when loc/scale is large.
### Reproducing code example:
```
>>> import scipy.stats as stats
>>> stats.norm.expect(loc=10, scale=1)
10.000000000000004
>>> stats.norm.expect(loc=20, scale=1)
20.0
>>> stats.norm.expect(loc=30, scale=1)
30.0
>>> stats.norm.expect(loc=36, scale=1)
3.3637320653305685e-10
>>> stats.norm.expect(loc=40, scale=1)
1.5746385460369797e-26
>>> stats.norm.expect(loc=10, scale=0.1)
2.3378433179863228e-37
```


### Scipy/Numpy/Python version information:
```
Scipy 1.1.0
Numpy 1.14.3
```

Hints from issue discussion:
The `expect` method is implemented with `scipy.integrate.quad`.  When `func` is not given, it is assumed that `func(x) = x`.  `quad` is used to numerically integrate `func(x)*pdf(x)` over the support of the distribution.  For the normal distribution, the integral is over [-inf, inf].  Apparently this use of `quad` isn't robust enough to handle large `loc`.

The problem also occurs with other distributions.  For example,

```
In [52]: from scipy.stats import gamma, logistic

In [53]: gamma.expect(args=(145,))
Out[53]: 145.00000000000045

In [54]: gamma.expect(args=(148,))  # Result should be 148.
Out[54]: 5.607236961785617e-09

In [55]: logistic.expect(loc=80)
Out[55]: 80.00000000000004

In [56]: logistic.expect(loc=85)  # Result should be 85.
Out[56]: 6.13219419324224e-09
```

Perhaps it could fixed by integrating over multiple intervals:  a finite interval [c, d] where c and d are based on the given values of `loc`, `scale` and the shape parameters (for distributions that have them), plus integrals over [a, c] and [d, b], where [a, b] is the support of the distribution (e.g. a = -inf and b = inf for the normal distribution).

Or, perhaps the variable of integration could be changed to `z = (x - loc)/scale`.  This wouldn't help for the `gamma` example, where `loc=0` and `scale=1`.



Note that the correct value can be computed by manually setting the integration limits:

```
In [73]: norm.expect(loc=50, scale=1)  # Bad result.
Out[73]: 9.38568238272931e-98

In [74]: norm.expect(loc=50, scale=1, lb=-50, ub=150)
Out[74]: 50.0

In [75]: gamma.expect(args=(148,))  # Bad result.
Out[75]: 5.607236961785617e-09

In [76]: gamma.expect(args=(148,), lb=0, ub=400)
Out[76]: 148.000000000001

In [77]: logistic.expect(loc=85)  # Bad result.
Out[77]: 6.13219419324224e-09

In [78]: logistic.expect(loc=85, lb=-50, ub=250)
Out[78]: 84.99999999999999
```
Could there be a way of throwing a warning or similar when this kind of inconsistencies occur? At least users would be aware that there might be problems with the obtained results.
The same issue also occurs when using a small loc and a large upper or lower bound: 

```
>>> norm.expect(loc=1, scale=1, ub=10, conditional=True)
0.9999999999982133

>>> norm.expect(loc=1, scale=1, ub=20, conditional=True)
0.9999999999999999

>>> norm.expect(loc=1, scale=1, ub=30, conditional=True)
1.0

>>> norm.expect(loc=1, scale=1, ub=40, conditional=True)
3.4007656189526604e-22

>>> norm.expect(loc=1, scale=1, ub=50, conditional=True)
8.381344573482554e-89
```

do we have any update on the bug fixes? 
I'm submitting a PR that breaks up the integration interval into lower tail, body, and upper tail; it also un-silences warnings raised by the integrator. It will make `expect` slower, but more reliable.
Wow, actually, this tends to speed things up, at least on average among the test cases:
```python3

import scipy.stats as stats
from scipy.stats import gamma, logistic

def f():
    a = stats.norm.expect(loc=10, scale=1)
    b = stats.norm.expect(loc=20, scale=1)
    c = stats.norm.expect(loc=30, scale=1)
    d = stats.norm.expect(loc=36, scale=1)
    e = stats.norm.expect(loc=40, scale=1)
    f = stats.norm.expect(loc=10, scale=0.1)
    g = gamma.expect(args=(145,))
    h = gamma.expect(args=(148,))  # Result should be 148.
    i = logistic.expect(loc=80)
    j = logistic.expect(loc=85)  # Result should be 85.
    return a, b, c, d, e, f, g, h, i, j

print(f())  # half incorrect in main, all correct in PR
%timeit f()  # 362ms in main, 207 ms in PR
```

Surprisingly, it even speeds things up for the problems that `main` currently gets correct.

```python3
def f():
    a = stats.norm.expect(loc=10, scale=1)
    b = stats.norm.expect(loc=20, scale=1)
    c = stats.norm.expect(loc=30, scale=1)
    # d = stats.norm.expect(loc=36, scale=1)
    # e = stats.norm.expect(loc=40, scale=1)
    # f = stats.norm.expect(loc=10, scale=0.1)
    g = gamma.expect(args=(145,))
    # h = gamma.expect(args=(148,))  # Result should be 148.
    i = logistic.expect(loc=80)
    # j = logistic.expect(loc=85)  # Result should be 85.
    # return a, b, c, d, e, f, g, h, i, j
    return a, b, c, g, i

%timeit f()  # 322 ms in main, 101 ms in PR
```

So I take back what I wrote; I'm in favor of a fix like this. I'll see if I can speed it up more with calling `quad` with `points` and a single call to `ppf`.

`points` won't work for functions with infinite support.
```ValueError: Infinity inputs cannot be used with break points.```

Commit message: TST: stats.rv_continuous_expect: add tests to demonstrate improvement

Patch to review:
diff --git a/scipy/stats/_distn_infrastructure.py b/scipy/stats/_distn_infrastructure.py
index f3d90e044f89..c2dfd45b9356 100644
--- a/scipy/stats/_distn_infrastructure.py
+++ b/scipy/stats/_distn_infrastructure.py
@@ -2885,6 +2885,13 @@ def expect(self, func=None, args=(), loc=0, scale=1, lb=None, ub=None,
         finite. For example ``cauchy(0).mean()`` returns ``np.nan`` and
         ``cauchy(0).expect()`` returns ``0.0``.
 
+        Likewise, the accuracy of results are not verified by the function.
+        `scipy.integrate.quad` is typically reliable for integrals that are
+        numerically favorable, but it is not guaranteed to converge
+        to a correct value for all possible intervals and integrands. This
+        function is provided for convenience; for critical applications,
+        check results against other integration methods.
+
         The function is not vectorized.
 
         Examples
@@ -2928,9 +2935,16 @@ def fun(x, *args):
         else:
             invfac = 1.0
         kwds['args'] = args
-        # Silence floating point warnings from integration.
-        with np.errstate(all='ignore'):
-            vals = integrate.quad(fun, lb, ub, **kwds)[0] / invfac
+        # split the interval to help integrator; see gh-8928
+        alpha = 0.05  # split body from tails at `alpha` probability mass
+        # considered separate calls to _ppf/_isf, but accuracy is not needed
+        c, d = loc + self._ppf(np.array([alpha, 1-alpha]), *args)*scale
+        # Do not silence warnings from integration.
+        lbc = integrate.quad(fun, lb, c, **kwds)[0]
+        cd = integrate.quad(fun, c, d, **kwds)[0]
+        dub = integrate.quad(fun, d, ub, **kwds)[0]
+        vals = (lbc + cd + dub) / invfac
+
         return vals
 
     def _param_info(self):
diff --git a/scipy/stats/tests/common_tests.py b/scipy/stats/tests/common_tests.py
index 5e9526e8ff0d..cb4ec4985810 100644
--- a/scipy/stats/tests/common_tests.py
+++ b/scipy/stats/tests/common_tests.py
@@ -24,9 +24,7 @@ def check_normalization(distfn, args, distname):
     norm_moment = distfn.moment(0, *args)
     npt.assert_allclose(norm_moment, 1.0)
 
-    # this is a temporary plug: either ncf or expect is problematic;
-    # best be marked as a knownfail, but I've no clue how to do it.
-    if distname == "ncf":
+    if distname == "rv_histogram_instance":
         atol, rtol = 1e-5, 0
     else:
         atol, rtol = 1e-7, 1e-7
@@ -66,10 +64,10 @@ def check_mean_expect(distfn, arg, m, msg):
 
 
 def check_var_expect(distfn, arg, m, v, msg):
+    kwargs = {'rtol': 5e-6} if msg=="rv_histogram_instance" else {}
     if np.isfinite(v):
         m2 = distfn.expect(lambda x: x*x, arg)
-        npt.assert_almost_equal(m2, v + m*m, decimal=5, err_msg=msg +
-                            ' - 2st moment (expect)')
+        npt.assert_allclose(m2, v + m*m, **kwargs)
 
 
 def check_skew_expect(distfn, arg, m, v, s, msg):
diff --git a/scipy/stats/tests/test_continuous_basic.py b/scipy/stats/tests/test_continuous_basic.py
index 40bb0e204f49..2d023d3b56de 100644
--- a/scipy/stats/tests/test_continuous_basic.py
+++ b/scipy/stats/tests/test_continuous_basic.py
@@ -248,8 +248,8 @@ def test_levy_stable_random_state_property():
 
 
 def cases_test_moments():
-    fail_normalization = set(['vonmises'])
-    fail_higher = set(['vonmises', 'ncf'])
+    fail_normalization = set()
+    fail_higher = set(['ncf'])
 
     for distname, arg in distcont[:] + [(histogram_test_instance, tuple())]:
         if distname == 'levy_stable':
@@ -286,19 +286,23 @@ def test_moments(distname, arg, normalization_ok, higher_ok, is_xfailing):
     with npt.suppress_warnings() as sup:
         sup.filter(IntegrationWarning,
                    "The integral is probably divergent, or slowly convergent.")
+        sup.filter(IntegrationWarning,
+                   "The maximum number of subdivisions.")
+
         if is_xfailing:
             sup.filter(IntegrationWarning)
 
         m, v, s, k = distfn.stats(*arg, moments='mvsk')
 
-        if normalization_ok:
-            check_normalization(distfn, arg, distname)
+        with np.errstate(all="ignore"):
+            if normalization_ok:
+                check_normalization(distfn, arg, distname)
 
-        if higher_ok:
-            check_mean_expect(distfn, arg, m, distname)
-            check_skew_expect(distfn, arg, m, v, s, distname)
-            check_var_expect(distfn, arg, m, v, distname)
-            check_kurt_expect(distfn, arg, m, v, k, distname)
+            if higher_ok:
+                check_mean_expect(distfn, arg, m, distname)
+                check_skew_expect(distfn, arg, m, v, s, distname)
+                check_var_expect(distfn, arg, m, v, distname)
+                check_kurt_expect(distfn, arg, m, v, k, distname)
 
         check_moment(distfn, arg, m, v, distname)
 
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index f57acd4ddf46..3dff3483dce0 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -4369,9 +4369,10 @@ def func(x):
             return x
 
         gm = stats.gamma(a=2, loc=3, scale=4)
-        gm_val = gm.expect(func, lb=1, ub=2, conditional=True)
-        gamma_val = stats.gamma.expect(func, args=(2,), loc=3, scale=4,
-                                       lb=1, ub=2, conditional=True)
+        with np.errstate(invalid="ignore", divide="ignore"):
+            gm_val = gm.expect(func, lb=1, ub=2, conditional=True)
+            gamma_val = stats.gamma.expect(func, args=(2,), loc=3, scale=4,
+                                           lb=1, ub=2, conditional=True)
         assert_allclose(gm_val, gamma_val)
 
         p = stats.poisson(3, loc=4)
@@ -4535,6 +4536,15 @@ def poiss_moment5(mu):
             m5 = stats.poisson.moment(5, mu)
             assert_allclose(m5, poiss_moment5(mu), rtol=1e-10)
 
+    def test_challenging_cases_gh8928(self):
+        # Several cases where `expect` failed to produce a correct result were
+        # reported in gh-8928. Check that these cases have been resolved.
+        assert_allclose(stats.norm.expect(loc=36, scale=1.0), 36)
+        assert_allclose(stats.norm.expect(loc=40, scale=1.0), 40)
+        assert_allclose(stats.norm.expect(loc=10, scale=0.1), 10)
+        assert_allclose(stats.gamma.expect(args=(148,)), 148)
+        assert_allclose(stats.logistic.expect(loc=85), 85)
+
 
 class TestNct:
     def test_nc_parameter(self):


Code review comments:
diff_hunk: @@ -66,10 +64,10 @@ def check_mean_expect(distfn, arg, m, msg):
 
 
 def check_var_expect(distfn, arg, m, v, msg):
+    kwargs = {'rtol': 5e-6} if msg=="rv_histogram_instance" else {}
path: scipy/stats/tests/common_tests.py
text: @author:
```suggestion
    kwargs = {'rtol': 5e-6} if msg == "rv_histogram_instance" else {}
```