BUG: fix moments method to support arrays and list

.github/workflows/linux.yml

@@ -36,7 +36,7 @@ jobs:
         run: |      
           python3.7-dbg -c 'import sys; print("Python debug build:", hasattr(sys, "gettotalrefcount"))'
           python3.7-dbg -m pip install --upgrade pip setuptools wheel
-          python3.7-dbg -m pip install --upgrade numpy cython pytest pytest-xdist pybind11
+          python3.7-dbg -m pip install --upgrade numpy cython pytest hypothesis pytest-xdist pybind11
           python3.7-dbg -m pip install --upgrade mpmath gmpy2 pythran
           python3.7-dbg -m pip uninstall -y nose
           cd ..
@@ -81,7 +81,7 @@ jobs:
     - name: Install packages
       run: |
         pip install --user git+https://github.com/numpy/numpy.git
-        python3.10 -m pip install --user setuptools wheel cython pytest pybind11 pytest-xdist
+        python3.10 -m pip install --user setuptools wheel cython pytest hypothesis pybind11 pytest-xdist
         pip install --user git+https://github.com/serge-sans-paille/pythran.git
         python3.10 -m pip install -r mypy_requirements.txt
 

.github/workflows/macos.yml

@@ -77,7 +77,7 @@ jobs:
     - name: Install packages
       run: |
         pip install ${{ matrix.numpy-version }}
-        pip install setuptools wheel cython pytest pytest-xdist pybind11 pytest-xdist mpmath gmpy2 pythran
+        pip install setuptools wheel cython pytest hypothesis pytest-xdist pybind11 pytest-xdist mpmath gmpy2 pythran
 
     - name: Test SciPy
       run: |

.gitignore

@@ -93,6 +93,10 @@ pip-wheel-metadata
 .cache/*
 .pytest_cache/*
 
+# hypothesis cache #
+####################
+.hypothesis/*
+
 # mypy cache #
 ##############
 .mypy_cache/*

scipy/stats/_distn_infrastructure.py

@@ -1272,14 +1272,14 @@ def moment(self, n, *args, **kwds):
         """
         shapes, loc, scale = self._parse_args(*args, **kwds)
         args = np.broadcast_arrays(*(*shapes, loc, scale))
-        shapes, loc, scale = args[:-2], args[-2], args[-1]
+        *shapes, loc, scale = args
 
         i0 = np.logical_and(self._argcheck(*shapes), scale > 0)
         i1 = np.logical_and(i0, loc == 0)
         i2 = np.logical_and(i0, loc != 0)
 
         args = argsreduce(i0, *shapes, loc, scale)
-        shapes, loc, scale = args[:-2], args[-2], args[-1]
+        *shapes, loc, scale = args
 
         if (floor(n) != n):
             raise ValueError("Moment must be an integer.")
@@ -1316,7 +1316,7 @@ def moment(self, n, *args, **kwds):
                     j+=1
             mu, mu2, g1, g2 = mom
             args = argsreduce(loc != 0, *shapes, loc, scale, val)
-            shapes, loc, scale, val = args[:-3], args[-3], args[-2], args[-1]
+            *shapes, loc, scale, val = args
 
             res2 = zeros(loc.shape, dtype='d')
             fac = scale / loc

scipy/stats/tests/test_continuous_basic.py

@@ -2,6 +2,8 @@
 import numpy as np
 import numpy.testing as npt
 import pytest
+import hypothesis.extra.numpy as npst
+from hypothesis import given, strategies as st
 from pytest import raises as assert_raises
 from scipy.integrate import IntegrationWarning
 
@@ -741,7 +743,7 @@ def test_moments_with_array_gh12192_regression():
     vals3 = stats.norm.moment(n=2, loc=0, scale=-4)
     expected3 = np.nan
     npt.assert_equal(vals3, expected3)
-    npt.assert_(isinstance(vals3, expected3.__class__))
+    assert isinstance(vals3, expected3.__class__)
 
     # array loc with 0 entries and scale with invalid entries
     vals4 = stats.norm.moment(n=2, loc=[1, 0, 2], scale=[3, -4, -5])
@@ -762,7 +764,7 @@ def test_moments_with_array_gh12192_regression():
     vals7 = stats.chi.moment(n=2, df=1, loc=0, scale=0)
     expected7 = np.nan
     npt.assert_equal(vals7, expected7)
-    npt.assert_(isinstance(vals7, expected7.__class__))
+    assert isinstance(vals7, expected7.__class__)
 
     # array args, scalar loc, and scalar scale
     vals8 = stats.chi.moment(n=2, df=[1, 2, 3], loc=0, scale=0)
@@ -808,23 +810,54 @@ def test_broadcasting_in_moments_gh12192_regression():
     vals0 = stats.norm.moment(n=1, loc=np.array([1, 2, 3]), scale=[[1]])
     expected0 = np.array([[1., 2., 3.]])
     npt.assert_equal(vals0, expected0)
-    npt.assert_(vals0.shape == expected0.shape)
+    assert vals0.shape == expected0.shape
 
     vals1 = stats.norm.moment(n=1, loc=np.array([[1], [2], [3]]),
                               scale=[1, 2, 3])
     expected1 = np.array([[1., 1., 1.], [2., 2., 2.], [3., 3., 3.]])
     npt.assert_equal(vals1, expected1)
-    npt.assert_(vals1.shape == expected1.shape)
+    assert vals1.shape == expected1.shape
 
     vals2 = stats.chi.moment(n=1, df=[1., 2., 3.], loc=0., scale=1.)
     expected2 = np.array([0.79788456, 1.25331414, 1.59576912])
     npt.assert_allclose(vals2, expected2, rtol=1e-8)
-    npt.assert_(vals2.shape == expected2.shape)
+    assert vals2.shape == expected2.shape
 
     vals3 = stats.chi.moment(n=1, df=[[1.], [2.], [3.]], loc=[0., 1., 2.],
                              scale=[-1., 0., 3.])
     expected3 = np.array([[np.nan, np.nan, 4.39365368],
                           [np.nan, np.nan, 5.75994241],
                           [np.nan, np.nan, 6.78730736]])
     npt.assert_allclose(vals3, expected3, rtol=1e-8)
+    assert vals3.shape == expected3.shape
     npt.assert_(vals3.shape == expected3.shape)
+
+
+# A strategy to generate the shape of a one-dimensional array of length [0..7]
+# Subsequent draws will use the same shape to ensure compatibility.  For fancier
+# tricks see `npst.broadcastable_shapes()` or `mutually_broadcastable_shapes()`.
+shared_shapes = st.shared(st.tuples(st.integers(0, 7)))
+
+
+def scalar_or_array(**kwargs):
+    return st.floats(**kwargs) | npst.arrays(
+        dtype=float, shape=shared_shapes, elements=kwargs
+    )
+
+
+@given(
+    # Our moment `n`, between one and four (inclusive)
+    n=st.integers(1, 4),
+    # Our locations and scale are each either a scalar (Python-native) float,
+    # or an ndarray of floats.  If both are arrays, they will have the same shape.
+    loc=scalar_or_array(min_value=0, max_value=1000),  # what should the bounds be?
+    scale=scalar_or_array(min_value=-1, max_value=1000),
+)
+def test_moments_gh12192_regression(n, loc, scale):
+    got = stats.norm.moment(n, loc=loc, scale=scale)
+    ref = np.vectorize(stats.norm.moment, otypes=[float])(n, loc=loc, scale=scale)
+
+    if isinstance(got, np.ndarray):
+        assert got.shape == ref.shape
+        assert got.dtype == ref.dtype
+    np.testing.assert_allclose(got, ref)