Repository: scipy/scipy
Pull request: ENH: stats.geometric.entropy: implement analytical formula
<!-- 
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

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
Fixes #18226

#### What does this implement/fix?
<!--Please explain your changes.-->

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
ENH: stats.geometric.entropy: implement analytical formula 
The geometric distribution entropy is calculated numerically, and is sometimes incorrect.
```
ss.geom(0.0146).entropy()  # Should be: 5.219397961962304
/home/neil/.cache/pypoetry/virtualenvs/efax-ZssyUsLU-py3.11/lib/python3.11/site-packages/scipy/stats/_distn_infrastructure.py:3814: RuntimeWarning: expect(): sum did not converge
  warnings.warn('expect(): sum did not converge', RuntimeWarning)
Out[2]: 2.9609120114214393
```
It should be calculated analytically as:
```python
-np.log(p) - np.log1p(-p) * (1.0-p) / p
```

Hints from issue discussion:
This is a good suggestion and the analytical formula looks like a straight forward improvement. Would you be interested in creating a PR?

We had many PRs related to entropy calculations lately, see e.g. https://github.com/scipy/scipy/pull/17857 or https://github.com/scipy/scipy/pull/17742, they could help you especially for the creation of test cases.
Thanks for reporting. I changed the title and label because the method is performing as expected, but an enhancement PR would be most welcome.

```python3
# no problem when shape parameter is less extreme
import numpy as np
from scipy import stats
res = stats.geom(0.146).entropy()
ref = -np.log(p) - np.log1p(-p) * (1.0-p) / p
np.testing.assert_allclose(res, ref)

# warns about failure to converge by design
res = stats.geom(0.0146).entropy()
```
It might be preferable for it to raise an error rather than emit a warning, but that is a more general issue.

I'm curious - can you share the context in which you encountered this shortcoming? I've been interested to see how `entropy` is used given all the work on it lately.
> I changed the title and label because the method is performing as expected

Thanks, that makes sense.

> I'm curious - can you share the context in which you encountered this shortcoming? I've been interested to see how entropy is used given all the work on it lately.

Sure!  I have my own statistical library for exponential family distribution, [efax](https://github.com/NeilGirdhar/efax/).  I [test it against scipy](https://github.com/NeilGirdhar/efax/blob/7dd2fb05c720d35a6bae6ebc641bfe41194b0d50/tests/test_match_scipy.py#L21).

> This is a good suggestion and the analytical formula looks like a straight forward improvement. Would you be interested in creating a PR?

Okay, I'll see if I find the time.

Commit message: ENH: stats.geometric.entropy: implement analytical formula

Patch to review:
diff --git a/scipy/stats/_discrete_distns.py b/scipy/stats/_discrete_distns.py
index 9cfa6d76a04a..26e90e0b6ac3 100644
--- a/scipy/stats/_discrete_distns.py
+++ b/scipy/stats/_discrete_distns.py
@@ -441,6 +441,9 @@ def _stats(self, p):
         g2 = np.polyval([1, -6, 6], p)/(1.0-p)
         return mu, var, g1, g2
 
+    def _entropy(self, p):
+        return -np.log(p) - np.log1p(-p) * (1.0-p) / p
+
 
 geom = geom_gen(a=1, name='geom', longname="A geometric")
 
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 76c0d6892eda..765402ef17fd 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -999,6 +999,12 @@ def test_ppf_underflow(self):
         # this should not underflow
         assert_allclose(stats.geom.ppf(1e-20, 1e-20), 1.0, atol=1e-14)
 
+    def test_entropy(self):
+        # Simple tests of entropy.
+        h = stats.geom(np.array([0.0146, 0.2, 0.5, 0.8])).entropy()
+        expected_h = np.array([5.21939796, 2.50201212, 1.38629436, 0.62550303])
+        assert_allclose(h, expected_h)
+
 
 class TestPlanck:
     def setup_method(self):


Code review comments:
diff_hunk: @@ -999,6 +999,12 @@ def test_ppf_underflow(self):
         # this should not underflow
         assert_allclose(stats.geom.ppf(1e-20, 1e-20), 1.0, atol=1e-14)
 
+    def test_entropy(self):
+        # Simple tests of entropy.
+        h = stats.geom(np.array([0.0146, 0.2, 0.5, 0.8])).entropy()
+        expected_h = np.array([5.21939796, 2.50201212, 1.38629436, 0.62550303])
+        assert_allclose(h, expected_h)
path: scipy/stats/tests/test_distributions.py
text: @user1:
Let's target the test at the issue. The other cases would already pass in main, and `_check_private_entropy` in `test_discrete_basic` already checks for consistency between `_entropy` and `_pmf`.

```suggestion
    def test_entropy_gh18226(self):
        # gh-18226 reported that `geom.entropy` produced a warning and
        # inaccurate output for small p. Check that this is resolved.
        h = stats.geom(0.0146).entropy()
        assert_allclose(h, 5.219397961962308, rtol=1e-15)
```

Normally, I would include code that generates the reference value, but I will omit that here because I couldn't decide what the reference should be. Here is one possibility.

<details>

```python3
from mpmath import mp
mp.dps = 100
p = mp.mpf(0.0146)
def f(k):
    return p * (1 - p) ** (k - 1)
float(mp.nsum(lambda k: -f(k) * mp.log(f(k)), (1, mp.inf)))
```

</details>