BUG: stats: fix wrong shape output of _stats method for burr and fisk…

… distributions (#13238)

As described in gh-13234, 'burr' and 'fisk' distributions returned a bad shaped array when scalar arguments were passed to the 'moment' method leading to an API inconsistency. This was because the call to '_stats' method from the 'moment' method when n<5 and n>1 returned a bad shaped array. The '_stats' method (in 'burr' distribution) has been fixed to return floating points when scalar arguments are passed. Regression tests for it have also been added.

Closes gh-13234

Co-authored-by: Warren Weckesser <warren.weckesser@gmail.com>

scipy/stats/_continuous_distns.py

@@ -957,6 +957,8 @@ def _stats(self, c, d):
             lambda c, e1, e2, e3, e4, mu2_if_c: (
                 ((e4 - 4*e3*e1 + 6*e2*e1**2 - 3*e1**4) / mu2_if_c**2) - 3),
             fillvalue=np.nan)
+        if np.ndim(c) == 0:
+            return mu.item(), mu2.item(), g1.item(), g2.item()
         return mu, mu2, g1, g2
 
     def _munp(self, n, c, d):

scipy/stats/tests/test_continuous_basic.py

@@ -654,3 +654,11 @@ def test_methods_with_lists(method, distname, args):
     npt.assert_allclose(result,
                         [f(*v) for v in zip(x, *shape2, loc, scale)],
                         rtol=1e-14, atol=5e-14)
+
+
+def test_burr_fisk_moment_gh13234_regression():
+    vals0 = stats.burr.moment(1, 5, 4)
+    assert isinstance(vals0, float)
+
+    vals1 = stats.fisk.moment(1, 8)
+    assert isinstance(vals1, float)