Repository: scipy/scipy
Pull request: BUG: loglaplace moment should be non-negative.
The current implementation of scipy.stats.loglaplace.moments incorrectly returns a negative number when the requested moment does not exist. This commit fixes this bug and correctly returns np.inf instead.

A test case is added for parameter combinations that lead to non-finite moments. Non-central moment is allowed to return inf or nan as it is not well-defined when a lower-order non-central moment is non-finite.

The generic stats() implementation is updated to suppress floating point warnings caused by (inf-inf) operation.

<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Depending on your changes, you can skip CI operations and save time and energy: 
http://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
Closes gh-20001.

#### What does this implement/fix?
Return `np.inf` when the moment doesn't exist for certain parameter combinations.

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
BUG: scipy.stats.loglaplace may return negative moments
### Describe your issue.

When the $r$-th moment of `scipy.stats.LogLaplace` diverges, the moment function returns a negative number instead.

The error is due to the `_munp` method not checking arguments. The fix is to return `np.inf` when |r|>=c.

### Reproducing Code Example

```python
import scipy.stats
scipy.stats.loglaplace.mean(0.5)
# returns -0.3333333333333
```


### Error message

```shell
No error, but wrong result.
```


### SciPy/NumPy/Python version and system information

```shell
Current main branch.
```

Hints from issue discussion:

Commit message: BUG: loglaplace moment should be non-negative.

The current implementation of scipy.stats.loglaplace.moments incorrectly
returns a negative number when the requested moment does not exist. This
commit fixes this bug and correctly returns np.inf instead.

A test case is added for parameter combinations that lead to non-finite
moments. Non-central moment is allowed to return inf or nan as it is not
well-defined when a lower-order non-central moment is non-finite.

The generic stats() implementation is updated to suppress floating point
warnings caused by (inf-inf) operation.

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 779077fb8c1a..86c96f72f734 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -6373,7 +6373,9 @@ def _isf(self, q, c):
         return np.where(q > 0.5, (2.0*(1.0 - q))**(1.0/c), (2*q)**(-1.0/c))
 
     def _munp(self, n, c):
-        return c**2 / (c**2 - n**2)
+        # Return inf if abs(n) >= c
+        with np.errstate(divide='ignore'):
+            return c**2 / np.maximum(c**2 - n**2, 0)
 
     def _entropy(self, c):
         return np.log(2.0/c) + 1.0
diff --git a/scipy/stats/_distn_infrastructure.py b/scipy/stats/_distn_infrastructure.py
index 1d2f1304838b..6ca41ccc4629 100644
--- a/scipy/stats/_distn_infrastructure.py
+++ b/scipy/stats/_distn_infrastructure.py
@@ -1154,7 +1154,8 @@ def stats(self, *args, **kwds):
                         mu = self._munp(1, *goodargs)
                     if mu2 is None:
                         mu2p = self._munp(2, *goodargs)
-                        mu2 = mu2p - mu * mu
+                        with np.errstate(invalid='ignore'):
+                            mu2 = mu2p - mu * mu
                     with np.errstate(invalid='ignore'):
                         mu3 = (-mu*mu - 3*mu2)*mu + mu3p
                         g1 = mu3 / np.power(mu2, 1.5)
@@ -1169,7 +1170,8 @@ def stats(self, *args, **kwds):
                         mu = self._munp(1, *goodargs)
                     if mu2 is None:
                         mu2p = self._munp(2, *goodargs)
-                        mu2 = mu2p - mu * mu
+                        with np.errstate(invalid='ignore'):
+                            mu2 = mu2p - mu * mu
                     if g1 is None:
                         mu3 = None
                     else:
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index f05780198778..f0f68a890e2c 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -3078,6 +3078,23 @@ def test_isf(self):
                1151387.578354072, 1640845512466.0906]
         assert_allclose(stats.loglaplace.isf(q, c), ref, rtol=1e-14)
 
+    @pytest.mark.parametrize('r, c',
+                             [(1, [0.5, 1.0]),
+                              (2, [0.5, 1.0, 1.5, 2.0]),
+                              (3, [0.5, 1.0, 1.5, 2.0, 2.5, 3.0]),
+                              (4, [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0])])
+    def test_moment(self, r, c):
+        # r-th non-central moment is infinite if |r| >= c.
+        assert_almost_equal(stats.loglaplace.moment(r, c), np.inf)
+
+    @pytest.mark.parametrize('mom, c',
+                             [('m', [0.5, 1.0]),
+                              ('v', [0.5, 1.0, 1.5, 2.0]),
+                              ('s', [0.5, 1.0, 1.5, 2.0, 2.5, 3.0]),
+                              ('k', [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0])])
+    def test_stats(self, mom, c):
+        # r-th non-central moment is non-finite (inf or nan) if r >= c.
+        assert not np.any(np.isfinite(stats.loglaplace.stats(c, moments=mom)))
 
 class TestPowerlaw:
 


Code review comments:
diff_hunk: @@ -6373,7 +6373,9 @@ def _isf(self, q, c):
         return np.where(q > 0.5, (2.0*(1.0 - q))**(1.0/c), (2*q)**(-1.0/c))
 
     def _munp(self, n, c):
-        return c**2 / (c**2 - n**2)
+        # Return inf if abs(n) >= c
+        with np.errstate(divide='ignore'):
+            return c**2 / np.maximum(c**2 - n**2, 0)
path: scipy/stats/_continuous_distns.py
text: @user1:
Test for equality `c == abs(n)` if they are integers and return explicit `np.inf` instead of delegating to numpy for inf handling. 

@author:
Thanks for reviewing. I made a new commit with that change.