Repository: scipy/scipy
Pull request: ENH: stats.moment: enable non-central moment calculation
#### Reference issue
Closes gh-17749

## What does this implement/fix?

Described in gh-17749.

#### Additional information

No unit tests and documentation added, just opened for reviewers to comment on the API - hence it's a draft.

Problem statement:
ENH: Compute non centraled moments with `stats.moment`?
### Is your feature request related to a problem? Please describe.

To compute a non-central moment of order `n>2`, one can't use the current implementation, but has to reside to implement their own.

### Describe the solution you'd like.

I'd like the `stats.moment` function to perhaps support a `central: Bool` argument, that will tell it whether the input vectors should be centralized around their average or not.

### Describe alternatives you've considered.

Implement this myself:

```python
def nmoment(arr, n):
	return np.sum(arr**n)/len(arr)
```

### Additional context (e.g. screenshots, GIFs)

_No response_

Hints from issue discussion:
Haven't tested this yet, but it seems a fix wouldn't be that hard:

```patch
diff --git i/scipy/stats/_mstats_basic.py w/scipy/stats/_mstats_basic.py
index 145323c08..d25a6da5d 100644
--- i/scipy/stats/_mstats_basic.py
+++ w/scipy/stats/_mstats_basic.py
@@ -2549,7 +2549,7 @@ def winsorize(a, limits=None, inclusive=(True, True), inplace=False,
                                    upinc, contains_nan, nan_policy)
 
 
-def moment(a, moment=1, axis=0):
+def moment(a, moment=1, centralized=True, axis=0):
     """
     Calculates the nth moment about the mean for a sample.
 
@@ -2559,13 +2559,15 @@ def moment(a, moment=1, axis=0):
        data
     moment : int, optional
        order of central moment that is returned
+    centralized : bool, optional
+       If True, then the moments are computed centralized. Default is True.
     axis : int or None, optional
        Axis along which the central moment is computed. Default is 0.
        If None, compute over the whole array `a`.
 
     Returns
     -------
-    n-th central moment : ndarray or float
+    n-th (central) moment : ndarray or float
        The appropriate moment along the given axis or over all values if axis
        is None. The denominator for the moment calculation is the number of
        observations, no degrees of freedom correction is done.
@@ -2590,7 +2592,10 @@ def moment(a, moment=1, axis=0):
 
     # for array_like moment input, return a value for each.
     if not np.isscalar(moment):
-        mean = a.mean(axis, keepdims=True)
+        if centralized:
+            mean = a.mean(axis, keepdims=True)
+        else:
+            mean = 0
         mmnt = [_moment(a, i, axis, mean=mean) for i in moment]
         return ma.array(mmnt)
     else:
diff --git i/scipy/stats/_stats_py.py w/scipy/stats/_stats_py.py
index e3e5f2bff..00d8c9f6f 100644
--- i/scipy/stats/_stats_py.py
+++ w/scipy/stats/_stats_py.py
@@ -1046,14 +1046,14 @@ def moment(a, moment=1, axis=0, nan_policy='propagate'):
 
     Notes
     -----
-    The k-th central moment of a data sample is:
+    The k-th moment of a data sample is:
 
     .. math::
 
-        m_k = \frac{1}{n} \sum_{i = 1}^n (x_i - \bar{x})^k
+        m_k = \frac{1}{n} \sum_{i = 1}^n (x_i)^k
 
-    Where n is the number of samples and x-bar is the mean. This function uses
-    exponentiation by squares [1]_ for efficiency.
+    Where n is the number of samples. This function uses exponentiation by
+    squares [1]_ for efficiency.
 
     Note that, if `a` is an empty array (``a.size == 0``), array `moment` with
     one element (`moment.size == 1`) is treated the same as scalar `moment`
```
Hi, can I take up this issue if no one's working on it? I might need some help since this would be my first contribution to scipy.
May I suggest something simpler for a first contribution? For example, https://github.com/scipy/scipy/issues/17683 would involve writing unit tests for existing functions, which would be just one part of what you'd need to do for this issue. Also, can you find other example like the one reported in https://github.com/scipy/scipy/issues/17999 in which a stats function documentation refers to numbers that aren't produced by the code? Then, you could make an adjustment to the documentation, which is another component of what you'd need to do for this issue. Once you have some experience with those things, tackling a more complex issue like this would be easier.
Okay, thank you, I'll start with these then!
@mdhaber if I'd go and file a PR with the patch I pasted above, as a solution to this issue, do you think a separate unit test is needed to solve the functionality added here?
Unit tests and documentation are required for all new features. This is why I suggested starting with an issue that requires tackling just one of these components before something that requires implementation, documentation, and unit tests.  
I took a closer look. The suggested patch would remove the mention of "central" from the documentation. We wouldn't want to do that, since the default behavior must remain unchanged.

Instead of a boolean "centralized", why not a float `center`? Then the user can pass the mean if they already have it or 0 if they want the non-central moment. Then the only change to the documentation would be to redefine the x-bar symbol as the desired center, which is the mean by default.

Either way, an API change like this requires an email to the [mailing list](https://svn.scipy.org/scipylib/mailing-lists.html). Please suggest the proposed change there.
Link for reference: https://mail.python.org/archives/list/scipy-dev@python.org/thread/FGVXUAA5Q2326LSZHXHLDE6I2MTTFFVO/
Sounds good to me to add a float to define the mean. As for naming `center` seems ok. Alternative, `mode`? 
`mode` is too ambiguous IMO, `None` / float `mean` sounds best to me in the meantime, but we can also use the same `None` / float logic but for an argument named `center`.

Commit message: stats.moment: update docstring phrasing

Co-authored-by: Matt Haberland <mhaberla@calpoly.edu>

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index 290118c41b6e..ffbeb1e9dc35 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -1009,7 +1009,7 @@ def _moment_result_object(*args):
     _moment_result_object, n_samples=1, result_to_tuple=lambda x: (x,),
     n_outputs=_moment_outputs
 )
-def moment(a, moment=1, axis=0, nan_policy='propagate'):
+def moment(a, moment=1, axis=0, nan_policy='propagate', *, center=None):
     r"""Calculate the nth moment about the mean for a sample.
 
     A moment is a specific quantitative measure of the shape of a set of
@@ -1033,9 +1033,14 @@ def moment(a, moment=1, axis=0, nan_policy='propagate'):
           * 'raise': throws an error
           * 'omit': performs the calculations ignoring nan values
 
+    center : float or None, optional
+       The point about which moments are taken. This can be the sample mean,
+       the origin, or any other be point. If `None` (default) compute the
+       center as the sample mean.
+
     Returns
     -------
-    n-th central moment : ndarray or float
+    n-th moment about the `center` : ndarray or float
        The appropriate moment along the given axis or over all values if axis
        is None. The denominator for the moment calculation is the number of
        observations, no degrees of freedom correction is done.
@@ -1046,14 +1051,15 @@ def moment(a, moment=1, axis=0, nan_policy='propagate'):
 
     Notes
     -----
-    The k-th central moment of a data sample is:
+    The k-th moment of a data sample is:
 
     .. math::
 
-        m_k = \frac{1}{n} \sum_{i = 1}^n (x_i - \bar{x})^k
+        m_k = \frac{1}{n} \sum_{i = 1}^n (x_i - c)^k
 
-    Where n is the number of samples and x-bar is the mean. This function uses
-    exponentiation by squares [1]_ for efficiency.
+    Where `n` is the number of samples, and `c` is the center around which the
+    moment is calculated. This function uses exponentiation by squares [1]_ for
+    efficiency.
 
     Note that, if `a` is an empty array (``a.size == 0``), array `moment` with
     one element (`moment.size == 1`) is treated the same as scalar `moment`
@@ -1074,19 +1080,14 @@ def moment(a, moment=1, axis=0, nan_policy='propagate'):
     """
     a, axis = _chk_asarray(a, axis)
 
-    contains_nan, nan_policy = _contains_nan(a, nan_policy)
-
-    if contains_nan and nan_policy == 'omit':
-        a = ma.masked_invalid(a)
-        return mstats_basic.moment(a, moment, axis)
-
+    if center is None:
+        center = a.mean(axis, keepdims=True)
     # for array_like moment input, return a value for each.
     if not np.isscalar(moment):
-        mean = a.mean(axis, keepdims=True)
-        mmnt = [_moment(a, i, axis, mean=mean) for i in moment]
+        mmnt = [_moment(a, i, axis, mean=center) for i in moment]
         return np.array(mmnt)
     else:
-        return _moment(a, moment, axis)
+        return _moment(a, moment, axis, mean=center)
 
 
 # Moment with optional pre-computed mean, equal to a.mean(axis, keepdims=True)
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 79f0dfb2fe48..fb28e7e45823 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -3001,6 +3001,27 @@ def _assert_equal(self, actual, expect, *, shape=None, dtype=None):
             dtype = expect.dtype
         assert actual.dtype == dtype
 
+    def test_moment_center(self):
+        # mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))
+        y = stats.moment(self.testcase, 1, center=0)
+        assert_approx_equal(y, 2.5)
+        y = stats.moment(self.testcase, 1, center=2)
+        assert_approx_equal(y, 0.5)
+        y = stats.moment(self.testcase, 2, center=0)
+        assert_approx_equal(y, 7.5)
+        y = stats.moment(self.testcase, 2, center=1)
+        assert_approx_equal(y, 3.5)
+        y = stats.moment(self.testcase, 3, center=1)
+        assert_approx_equal(y, 9)
+        y = stats.moment(self.testcase, 3, center=2)
+        assert_approx_equal(y, 2)
+        # Compare against the naive np.sum implementation:
+        x = self.testmathworks
+        y = stats.moment(x, 3, center=0)
+        assert_approx_equal(y, np.sum(x**3)/len(x))
+        y = stats.moment(x, 2, center=0)
+        assert_approx_equal(y, np.sum(x**2)/len(x))
+
     def test_moment(self):
         # mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))
         y = stats.moment(self.scalar_testcase)


Code review comments:
diff_hunk: @@ -3001,6 +3001,27 @@ def _assert_equal(self, actual, expect, *, shape=None, dtype=None):
             dtype = expect.dtype
         assert actual.dtype == dtype
 
+    def test_moment_center(self):
+        # mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))
+        y = stats.moment(self.testcase, 1, center=0)
+        assert_approx_equal(y, 2.5)
+        y = stats.moment(self.testcase, 1, center=2)
path: scipy/stats/tests/test_stats.py
text: @user1:
To avoid repetition, please use [`pytest.mark.parametrize`](https://docs.pytest.org/en/6.2.x/parametrize.html).

diff_hunk: @@ -3001,6 +3001,27 @@ def _assert_equal(self, actual, expect, *, shape=None, dtype=None):
             dtype = expect.dtype
         assert actual.dtype == dtype
 
+    def test_moment_center(self):
+        # mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))
+        y = stats.moment(self.testcase, 1, center=0)
+        assert_approx_equal(y, 2.5)
path: scipy/stats/tests/test_stats.py
text: @user1:
Please use `np.assert_allclose`.

![image](https://user-images.githubusercontent.com/6570539/230737728-be003b48-3a6e-486a-a4ab-d08342a05926.png)