BUG: stats.truncnorm.moments yields error for moment order greater than 5

[huangtu577]

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

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

Error message

----> 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

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']

[mdhaber]

Thanks for reporting. I'll take a look. Could be as simple as making the list in this line 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?

[tupui]

Do we really want to go into that direction and consider this to be on scope? To me higher moment is quite niche. Usually it's already hard to converge a 4th order, so going higher?

If we really do that, maybe we should add some warnings both in the doc and raise some warnings if we think moments are not converged.

[mdhaber]

It's not a new direction. truncnorm was a simple bug fix. betaprime was the only distribution that raised NotImplementedError for higher moments. The integrator will warn if it has trouble.
As I mentioned, I can remove the betaprime implementation and generic test if the reviewer would prefer not to review them (but it's not a matter of scope). I was conflicted about adding those - on one hand, it concerned me that a basic test like that was missing, and on the other (of course) I don't want to spend time on fixing something I'll replace soon. What tipped the scales is that I thought I would take very little time, so I went ahead.

[tupui]

Mmm ok but if we don't test and document a behavior, we can argue that it's a fortunate accident that something extra is working.

[mdhaber]

What is not tested/documented?

[tupui]

I am confused now, I thought that's what you were saying 😅