Repository: scipy/scipy
Pull request: BUG: stats: fix wrong shape output of burr and fisk distributions
#### Reference issue

closes gh-13234

#### What does this implement/fix?

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.

#### Additional information

Maybe wait for gh-12197 that adds support for array arguments in the `moment` method. I think this fix should work with array arguments too (but there may be edge cases I haven't considered). I hope my approach isn't contentious!

Problem statement:
BUG: stats: Wrong shape of burr.moment() and fisk.moment() when given scalars.
The `moment` method of the distribution `scipy.stats.burr` returns a 1-d array of length 1 instead of returning a scalar:
```
In [36]: from scipy.stats import burr

In [37]: burr.moment(1, 5, 4)  # Result should be a scalar, but we get a 1-d array.
Out[37]: array([1.50509474])
```
Compare to, for example, `beta` and `gamma`:
```
In [38]: from scipy.stats import beta, gamma

In [39]: beta.moment(1, 3, 2)
Out[39]: 0.6

In [40]: gamma.moment(1, 4.5)
Out[40]: 4.5
```
SciPy version:
```
In [41]: import scipy

In [42]: scipy.__version__                                                                                            
Out[42]: '1.7.0.dev0+ec0541e'
```

Hints from issue discussion:
The distribution `fisk` has the same problem:
```
In [134]: from scipy.stats import fisk

In [135]: fisk.moment(1, 8)                                                                                           
Out[135]: array([1.02617215])
```

Commit message: MAINT: stats: apply suggestions from code review

Remove the call to `np.broadcast_arrays` and use `np.ndim` instead to get the dimensions of the argument.

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

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index cc7e090a7867..96b791df9d04 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/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):
diff --git a/scipy/stats/tests/test_continuous_basic.py b/scipy/stats/tests/test_continuous_basic.py
index 45e086c657ab..62633656d8db 100644
--- a/scipy/stats/tests/test_continuous_basic.py
+++ b/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)
+    npt.assert_(isinstance(vals0, float))
+
+    vals1 = stats.fisk.moment(1, 8)
+    npt.assert_(isinstance(vals1, float))


Code review comments:
diff_hunk: @@ -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)
+    npt.assert_(isinstance(vals0, float))
+
+    vals1 = stats.fisk.moment(1, 8)
+    npt.assert_(isinstance(vals1, float))
path: scipy/stats/tests/test_continuous_basic.py
text: @user1:
One more suggestion, that I'll push myself. These days we can use a plain `assert` instead of `npt.assert_`.
```suggestion
    assert isinstance(vals0, float)

    vals1 = stats.fisk.moment(1, 8)
    assert isinstance(vals1, float)
```

@author:
Ah, thanks for letting me know!