Repository: scipy/scipy
Pull request: MAINT: stats.truncnorm/stats.betaprime: fix _munp for higher moments, add tests
#### Reference issue
Closes gh-18634

#### What does this implement/fix?
gh-18634 reported 
> When using the stats.truncnorm.moment, I always get an error, if I try to compute a moment greater than [or equal to] 5. 

This fixes the bug and adds a distribution-specific tests. 

Since there was no existing test coverage, I wrote a generic test for all distributions that override `_munp`. This discovered that `betaprime` raised a `NotImplementedError` for higher moments. Now, rather than raising an error, it falls back to a generic implementation. 

#### Additional information
I'd like to backport at least the `truncnorm` fix; this does not require any specialized `stats` knowledge to review. The `_munp` function was relying on `_pdf` to accept a list; the fix was as simple as turning the list into an array.
I'd be happy to remove the `betaprime` part and generic test to make review easier, since nobody has complained about `betaprime`.

Problem statement:
BUG: stats.truncnorm.moments yields error for moment order greater than 5
### Describe your issue.

When using the stats.truncnorm.moment, I always get an error, if I try to compute a moment greater than 5

### Reproducing Code Example

```python
scipy.stats.truncnorm.moment(5, -2, 3)
```


### Error message

```shell
----> 1 truncnorm.moment(8, -2, 3)

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/stats/_distn_infrastructure.py:1344, in rv_generic.moment(self, order, *args, **kwds)
   1342     mu, mu2, g1, g2 = self._stats(*shapes, **mdict)
   1343 val = np.empty(loc.shape)  # val needs to be indexed by loc
-> 1344 val[...] = _moment_from_stats(n, mu, mu2, g1, g2, self._munp, shapes)
   1346 # Convert to transformed  X = L + S*Y
   1347 # E[X^n] = E[(L+S*Y)^n] = L^n sum(comb(n, k)*(S/L)^k E[Y^k], k=0...n)
   1348 result = zeros(i0.shape)

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/stats/_distn_infrastructure.py:394, in _moment_from_stats(n, mu, mu2, g1, g2, moment_func, args)
    392         val = mu4+4*mu*mu3+6*mu*mu*mu2+mu*mu*mu*mu
    393 else:
--> 394     val = moment_func(n, *args)
    396 return val

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8795, in truncnorm_gen._munp(self, n, a, b)
   8792         moments.append(mk)
   8793     return moments[-1]
-> 8795 return _lazywhere((n >= 0) & (a == a) & (b == b), (n, a, b),
   8796                   np.vectorize(n_th_moment, otypes=[np.float64]),
   8797                   np.nan)

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/_lib/_util.py:69, in _lazywhere(cond, arrays, f, fillvalue, f2)
     67 tcode = np.mintypecode([a.dtype.char for a in arrays])
     68 out = np.full(np.shape(arrays[0]), fill_value=fillvalue, dtype=tcode)
---> 69 np.place(out, cond, f(*temp))
     70 if f2 is not None:
     71     temp = tuple(np.extract(~cond, arr) for arr in arrays)

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/numpy/lib/function_base.py:2304, in vectorize.__call__(self, *args, **kwargs)
   2301     vargs = [args[_i] for _i in inds]
   2302     vargs.extend([kwargs[_n] for _n in names])
-> 2304 return self._vectorize_call(func=func, args=vargs)

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/numpy/lib/function_base.py:2387, in vectorize._vectorize_call(self, func, args)
   2384 # Convert args to object arrays first
   2385 inputs = [asanyarray(a, dtype=object) for a in args]
-> 2387 outputs = ufunc(*inputs)
   2389 if ufunc.nout == 1:
   2390     res = asanyarray(outputs, dtype=otypes[0])

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8781, in truncnorm_gen._munp.<locals>.n_th_moment(n, a, b)
   8776 def n_th_moment(n, a, b):
   8777     """
   8778     Returns n-th moment. Defined only if n >= 0.
   8779     Function cannot broadcast due to the loop over n
   8780     """
-> 8781     pA, pB = self._pdf([a, b], a, b)
   8782     probs = [pA, -pB]
   8783     moments = [0, 1]

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8691, in truncnorm_gen._pdf(self, x, a, b)
   8690 def _pdf(self, x, a, b):
-> 8691     return np.exp(self._logpdf(x, a, b))

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8694, in truncnorm_gen._logpdf(self, x, a, b)
   8693 def _logpdf(self, x, a, b):
-> 8694     return _norm_logpdf(x) - _log_gauss_mass(a, b)

File ~/anaconda3/envs/test_scipy/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:296, in _norm_logpdf(x)
    295 def _norm_logpdf(x):
--> 296     return -x**2 / 2.0 - _norm_pdf_logC

TypeError: unsupported operand type(s) for ** or pow(): 'list' and 'int'
```


### SciPy/NumPy/Python version and system information

```shell
1.10.1 1.22.3 sys.version_info(major=3, minor=9, micro=16, releaselevel='final', serial=0)
lapack_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/include']
lapack_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/include']
blas_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/include']
blas_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/huangtu577/anaconda3/envs/test_scipy/include']
```

Hints from issue discussion:
Thanks for reporting. I'll take a look. Could be as simple as making the list in [this line](https://github.com/scipy/scipy/blame/1305da5985cbb1524aaa58b1296647bb122e6249/scipy/stats/_continuous_distns.py#L9379) an array. That code wasn't modified in the recent `truncnorm` changes, and it looks like it relied on `_pdf` accepting a list.

I'm always curious about what these things are used for. Are you able to describe what you're working on?

Commit message: MAINT: stats.truncnorm/stats.betaprime: fix _munp for higher moments, add tests

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 279b97efa836..3cf386261ab0 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -964,7 +964,10 @@ def _munp(self, n, a, b):
                               ((b - 4.0)*(b - 3.0)*(b - 2.0)*(b - 1.0))),
                 fillvalue=np.inf)
         else:
-            raise NotImplementedError
+            return _lazywhere(b > n, (a, b),
+                              # sup=super() prevents late-binding issue
+                              lambda a, b, sup=super(): sup._munp(n, a, b),
+                              fillvalue=np.inf)
 
 
 betaprime = betaprime_gen(a=0.0, name='betaprime')
@@ -9376,7 +9379,7 @@ def n_th_moment(n, a, b):
             Returns n-th moment. Defined only if n >= 0.
             Function cannot broadcast due to the loop over n
             """
-            pA, pB = self._pdf([a, b], a, b)
+            pA, pB = self._pdf(np.asarray([a, b]), a, b)
             probs = [pA, -pB]
             moments = [0, 1]
             for k in range(1, n+1):
diff --git a/scipy/stats/tests/common_tests.py b/scipy/stats/tests/common_tests.py
index e8a2fbb8435a..53718db339e8 100644
--- a/scipy/stats/tests/common_tests.py
+++ b/scipy/stats/tests/common_tests.py
@@ -88,6 +88,19 @@ def check_kurt_expect(distfn, arg, m, v, k, msg):
         npt.assert_(np.isnan(k))
 
 
+def check_munp_expect(dist, args, msg):
+    # If _munp is overridden, test a higher moment. (Before gh-18634, some
+    # distributions had issues with moments 5 and higher.)
+    if dist._munp.__func__ != stats.rv_continuous._munp:
+        res = dist.moment(5, *args)  # shouldn't raise an error
+        ref = dist.expect(lambda x: x ** 5, args, lb=-np.inf, ub=np.inf)
+        if not np.isfinite(res):  # could be valid; automated test can't know
+            return
+        # loose tolerance, mostly to see whether _munp returns *something*
+        assert_allclose(res, ref, atol=1e-10, rtol=1e-4,
+                        err_msg=msg + ' - higher moment / _munp')
+
+
 def check_entropy(distfn, arg, msg):
     ent = distfn.entropy(*arg)
     npt.assert_(not np.isnan(ent), msg + 'test Entropy is nan')
diff --git a/scipy/stats/tests/test_continuous_basic.py b/scipy/stats/tests/test_continuous_basic.py
index c25a0bb30389..1beac9747b47 100644
--- a/scipy/stats/tests/test_continuous_basic.py
+++ b/scipy/stats/tests/test_continuous_basic.py
@@ -17,7 +17,7 @@
                            check_meth_dtype, check_ppf_dtype,
                            check_cmplx_deriv,
                            check_pickling, check_rvs_broadcast,
-                           check_freezing,)
+                           check_freezing, check_munp_expect,)
 from scipy.stats._distr_params import distcont
 from scipy.stats._distn_infrastructure import rv_continuous_frozen
 
@@ -318,6 +318,8 @@ def test_moments(distname, arg, normalization_ok, higher_ok, moment_ok,
                    "The integral is probably divergent, or slowly convergent.")
         sup.filter(IntegrationWarning,
                    "The maximum number of subdivisions.")
+        sup.filter(IntegrationWarning,
+                   "The algorithm does not converge.")
 
         if is_xfailing:
             sup.filter(IntegrationWarning)
@@ -333,6 +335,7 @@ def test_moments(distname, arg, normalization_ok, higher_ok, moment_ok,
                 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_munp_expect(distfn, arg, distname)
 
         check_loc_scale(distfn, arg, m, v, distname)
 
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 9dafb8e8b772..54dac10f26d2 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -1494,6 +1494,15 @@ def test_logcdf_gh17064(self):
         assert_allclose(stats.truncnorm(a, b).logcdf(x), expected)
         assert_allclose(stats.truncnorm(-b, -a).logsf(-x), expected)
 
+    def test_moments_gh18634(self):
+        # gh-18634 reported that moments 5 and higher didn't work; check that
+        # this is resolved
+        res = stats.truncnorm(-2, 3).moment(5)
+        # From Mathematica:
+        # Moment[TruncatedDistribution[{-2, 3}, NormalDistribution[]], 5]
+        ref = 1.645309620208361
+        assert_allclose(res, ref)
+
 
 class TestGenLogistic:
 


Code review comments:
diff_hunk: @@ -964,7 +964,10 @@ def _munp(self, n, a, b):
                               ((b - 4.0)*(b - 3.0)*(b - 2.0)*(b - 1.0))),
                 fillvalue=np.inf)
         else:
-            raise NotImplementedError
+            return _lazywhere(b > n, (a, b),
+                              # sup=super() prevents late-binding issue
+                              lambda a, b, sup=super(): sup._munp(n, a, b),
+                              fillvalue=np.inf)
path: scipy/stats/_continuous_distns.py
text: @author:
We should just use the [general formula](https://en.wikipedia.org/wiki/Beta_prime_distribution) here.

![image](https://github.com/scipy/scipy/assets/6570539/f7002338-2939-4082-ae6b-ceaade090942)