Repository: scipy/scipy
Pull request: MAINT: stats._moment: warn when catastrophic cancellation occurs
#### Reference issue
Closes gh-15554
gh-13254
gh-13245
gh-11086
gh-10896

#### What does this implement/fix?
We get reports every once in a while of `skew` and/or `kurtosis` being unreliable with degenerate input. This PR would report a warning rather than trying to produce a correct result. Alternatively, we could replace the values that suffer from severe precision loss with `nan`. Re-reading gh-13254, I still don't think `0` is necessarily correct.


Problem statement:
BUG: stats.skew: unreliable when all values in axis slice are identical
### Describe your issue.

computing the skewness of a numpy ndarray will produce two result with following method for a column in ndarray (x):
1. skew(x[:, i])
2. skew(x)[i]

the issues happens when this column has equal float value


### Reproducing Code Example

```python
import numpy as np
from scipy.stats import skew

a = np.random.randn(100, 10)
a[:, 0] = 1.01

print(f'result in matrix value {skew(a)[0]}')   # -1.0
print(f'result in array value {skew(a[:, 0])}') # 0.0
```


### Error message

```shell
No error message
```


### SciPy/NumPy/Python version information

1.7.3 1.22.0 sys.version_info(major=3, minor=8, micro=5, releaselevel='final', serial=0)

Hints from issue discussion:
Hi @ronzhou, thanks for reporting. I can reproduce and indeed there is an issue here.

It's a resolution issue when we compute the moments with `_moment`.

```python
from scipy.stats._stats_py import _moment

_moment(a[:, 0], 2, axis, mean=a[:, 0].mean(axis=0, keepdims=True))
# 4.930380657631324e-32

_moment(a, 2, axis, mean=a.mean(axis=0, keepdims=True))[0]
# 1.7749370367472766e-30
```

And this slight difference	makes it fail here:

```python
zero = (m2 <= (np.finfo(m2.dtype).resolution * mean.squeeze(axis))**2)
vals = np.where(zero, 0, m3 / m2**1.5)
```

And the root cause is from NumPy actually. 

```python
import numpy.testing as npt
npt.assert_equal(a.mean(axis=0, keepdims=True)[0][0], a[:, 0].mean(axis=0, keepdims=True)[0])

# Items are not equal:
#  ACTUAL: 1.0100000000000013
#  DESIRED: 1.0100000000000002

```
From NumPy's doc

> For floating point numbers the numerical precision of sum (and np.add.reduce) is in general limited by directly adding each number individually to the result causing rounding errors in every step. However, often numpy will use a numerically better approach (partial pairwise summation) leading to improved precision in many use-cases. This improved precision is always provided when no axis is given. When axis is given, it will depend on which axis is summed. Technically, to provide the best speed possible, the improved precision is only used when the summation is along the fast axis in memory. Note that the exact precision may vary depending on other parameters.

So in the end we have to update our check. I am not sure where this (`np.finfo(m2.dtype).resolution * mean.squeeze(axis))**2`) is coming from. 

@peterbell10 I see you touch this area recently. Would you know? Can we bump the tolerance a bit here? I always thought it does not make much sense to go bellow the machine precision.
Hi @tupui , thanks for your response. 

In a 3 dim array, this is not an issue.

```
import numpy as np
from scipy.stats import skew

a = np.random.randn(10, 5, 3)
a[0, :, 0] = 1.01
print(skew(a, axis =1)[0, 0]) # 0.0
```

I use scipy.stats.skew to compute skewness along axis 0 in a 2d array. Currently, I use the method below to avoid this issue,

```
a = np.random.randn(100, 10)
a[:, 0] = 1.01

result = skew(a) 
result[np.nanmin(a, axis = 0) == np.nanmax(a, axis = 0)] = 0
```

Is there any better solution to avoid it temporarily? Thanks for your help! 
> Is there any better solution to avoid it temporarily? Thanks for your help!

Better I don't know but you can try playing with the resolution if your application allows that: `skew(a.astype('float32'))`

Commit message: MAINT: stats: ignore invalid comparison w/ nans

[skip actions][skip circleci]

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index db49bdcadee1..7ae69eac2482 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -1013,6 +1013,16 @@ def _moment(a, moment, axis, *, mean=None):
         # Starting point for exponentiation by squares
         mean = a.mean(axis, keepdims=True) if mean is None else mean
         a_zero_mean = a - mean
+
+        eps = np.finfo(a_zero_mean.dtype).resolution * 1e6
+        rel_diff = np.max(np.abs(a_zero_mean), axis=axis, keepdims=True)/mean
+        with np.errstate(invalid='ignore'):
+            precision_loss = np.any((0 < rel_diff) & (rel_diff < eps))
+        if precision_loss:
+            message = ("Precision loss occurred in moment calculation due to "
+                       "catastrophc cancellation. Results may be unreliable.")
+            warnings.warn(message, RuntimeWarning, stacklevel=4)
+
         if n_list[-1] == 1:
             s = a_zero_mean.copy()
         else:
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index c73637ca5541..9552383dcb61 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -3001,15 +3001,16 @@ def test_skew_propagate_nan(self):
     def test_skew_constant_value(self):
         # Skewness of a constant input should be zero even when the mean is not
         # exact (gh-13245)
-        a = np.repeat(-0.27829495, 10)
-        assert stats.skew(a) == 0.0
-        assert stats.skew(a * float(2**50)) == 0.0
-        assert stats.skew(a / float(2**50)) == 0.0
-        assert stats.skew(a, bias=False) == 0.0
+        with pytest.warns(RuntimeWarning, match="Precision loss occurred"):
+            a = np.repeat(-0.27829495, 10)
+            assert stats.skew(a) == 0.0
+            assert stats.skew(a * float(2**50)) == 0.0
+            assert stats.skew(a / float(2**50)) == 0.0
+            assert stats.skew(a, bias=False) == 0.0
 
-        # similarly, from gh-11086:
-        assert stats.skew([14.3]*7) == 0.0
-        assert stats.skew(1 + np.arange(-3, 4)*1e-16) == 0
+            # similarly, from gh-11086:
+            assert stats.skew([14.3]*7) == 0.0
+            assert stats.skew(1 + np.arange(-3, 4)*1e-16) == 0
 
     def test_kurtosis(self):
         # Scalar test case
@@ -3066,6 +3067,16 @@ def test_moment_accuracy(self):
         assert_allclose(np.power(tc_no_mean, 42).mean(),
                             stats.moment(self.testcase_moment_accuracy, 42))
 
+    def test_precision_loss_gh15554(self):
+        # gh-15554 was one of several issues that have reported problems with
+        # constant or near-constant input. We can't always fix these, but
+        # make sure there's a warning.
+        with pytest.warns(RuntimeWarning, match="Precision loss occurred"):
+            rng = np.random.default_rng(34095309370)
+            a = rng.random(size=(100, 10))
+            a[:, 0] = 1.01
+            stats.skew(a)[0]
+
 
 class TestStudentTest:
     X1 = np.array([-1, 0, 1])


Code review comments:
diff_hunk: @@ -1013,6 +1013,16 @@ def _moment(a, moment, axis, *, mean=None):
         # Starting point for exponentiation by squares
         mean = a.mean(axis, keepdims=True) if mean is None else mean
         a_zero_mean = a - mean
+
+        eps = np.finfo(a_zero_mean.dtype).resolution * 1e6
+        rel_diff = np.max(np.abs(a_zero_mean), axis=axis, keepdims=True)/mean
+        with np.errstate(invalid='ignore'):
+            precision_loss = np.any((0 < rel_diff) & (rel_diff < eps))
+        if precision_loss:
+            message = ("Precision loss occurred in moment calculation due to "
+                       "catastrophc cancellation. Results may be unreliable.")
path: scipy/stats/_stats_py.py
text: @user1:
I think the message is just missing an information for the user. Is this actionable or not for them. If you get the warning, then what? I would add a line in that sense.

@author:
I'm not sure what to recommend to them, though. Do you have a suggestion for what to do?

@user1:
Maybe we can recommend to play with dtype or say in which condition this usually happen?

@author:
> Precision loss occurred in moment calculation due to catastrophc cancellation. This occurs when the data are nearly identical. Results may be unreliable.

?

That is the condition in which it usually happens. I don't know that there is any relationship with `dtype` - seems like it could happen with `np.float32` just as much as `np.float64`.

@user1:
Better yes 👍 

@author:
Committed. Thanks! That will be helpful.

diff_hunk: @@ -1013,6 +1013,16 @@ def _moment(a, moment, axis, *, mean=None):
         # Starting point for exponentiation by squares
         mean = a.mean(axis, keepdims=True) if mean is None else mean
         a_zero_mean = a - mean
+
+        eps = np.finfo(a_zero_mean.dtype).resolution * 1e6
+        rel_diff = np.max(np.abs(a_zero_mean), axis=axis, keepdims=True)/mean
+        with np.errstate(invalid='ignore'):
+            precision_loss = np.any((0 < rel_diff) & (rel_diff < eps))
+        if precision_loss:
+            message = ("Precision loss occurred in moment calculation due to "
+                       "catastrophc cancellation. Results may be unreliable.")
path: scipy/stats/_stats_py.py
text: @author:
```suggestion
                       "catastrophic cancellation. Results may be unreliable.")
```

diff_hunk: @@ -1013,6 +1013,16 @@ def _moment(a, moment, axis, *, mean=None):
         # Starting point for exponentiation by squares
         mean = a.mean(axis, keepdims=True) if mean is None else mean
         a_zero_mean = a - mean
+
+        eps = np.finfo(a_zero_mean.dtype).resolution * 1e6
+        rel_diff = np.max(np.abs(a_zero_mean), axis=axis, keepdims=True)/mean
path: scipy/stats/_stats_py.py
text: @author:
@user1 thanks for your help today. This issue is related to that one (gh-10821) in that it also involves catastrophic cancellation when computing a moment from data.

In gh-10821, you looked for the variance being zero (`< np.finfo(dtype).tiny`) and did something to compensate. When this happened, you didn't really need to calculate the variance very accurately because it signaled that the optimization need not consider that value of `lmbda`.

This is general n-th order moment computation code, so there's not much we can do to compensate when precision loss occurs. So here, I'm looking for the relative difference between data and the mean being less than some tolerance, and I just want to warn the user that some sort of precision loss has occurred. I could probably be more specific about how much precision loss and exactly where it has occurred, but that can be left to a follow-up.

Do you have any thoughts on this PR?

@user1:
It's a trade off to me. Adding a precision loss check adds more computational, but it gives warnings when something looks wrong. My concern is the heuristic for checking the precision loss and making sure that does not warn for valid input. For example, how was `1e6` decided when computing `eps`?

As for the PR, if you want the bare mimumum change, one can directly detect the near constant case in `skew` and raise a specific warning for that.

@author:
Thanks @user1!  We already have an exact equality check in `skew` and `kurtosis`, but afterward we got a another report for which that check was insufficient. Also, we've gotten the same sort of issue in a few contexts, so I thought I'd treat more broadly. 
`1e6` was arbitrary, but it corresponds with a loss of approximately 10 decimal digits of precision for `float64`. This may be too noisy, so I'll reduce it so that this only warns when there is almost complete loss of precision.
Thanks!

diff_hunk: @@ -1013,6 +1013,16 @@ def _moment(a, moment, axis, *, mean=None):
         # Starting point for exponentiation by squares
         mean = a.mean(axis, keepdims=True) if mean is None else mean
         a_zero_mean = a - mean
+
+        eps = np.finfo(a_zero_mean.dtype).resolution * 1e6
+        rel_diff = np.max(np.abs(a_zero_mean), axis=axis, keepdims=True)/mean
+        with np.errstate(invalid='ignore'):
+            precision_loss = np.any((0 < rel_diff) & (rel_diff < eps))
path: scipy/stats/_stats_py.py
text: @author:
```suggestion
            precision_loss = np.any(rel_diff < eps)
```
We probably _do_ want this warning to show when all the data is exactly numerically equal. Related: I don't think it was correct in gh-13254 to set the `skew`/`kurtosis` to zero when all data were equal; it should be `nan`, and there should be a warning.

@user1:
> We probably do want this warning to show when all the data is exactly numerically equal.

Yes, agreed.

> I don't think it was correct in [gh-13254](https://github.com/scipy/scipy/pull/13254) to set the `skew`/`kurtosis` to zero when all data were equal; it should be `nan`, and there should be a warning.

I agree: 0 is not the correct skew/kurtosis for near-constant inputs.