Repository: scipy/scipy
Pull request: TST/DOC: stats: explain/check 100% interval for discrete distributions
#### Reference issue
Closes gh-11026

#### What does this implement/fix?
gh-11026 reported that discrete distributions returned incorrect values for `interval` with `confidence=1` (i.e. 100% confidence interval). They are not incorrect, per se, but it is somewhat surprising that the left end of the interval is beyond the support (as explained [here](https://github.com/scipy/scipy/issues/11026#issuecomment-1094200367)). This adds a note about this to the docs and adds a test to check for consistent behavior across discrete distributions. 


Problem statement:
rv_discrete.interval returning wrong values for alpha = 1
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

Scipy's rv_discrete.interval method, it appears that the lower bound is decreased by 1 when the alpha is set to exactly 1 (as far as I can tell, this is the most reliable way to get a generic rv's support's minimum)

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```
from scipy.stats import bernoulli  # also tested with randint and binom

rv = bernoulli(0.5)

print(rv.interval(1))  # expected (0,1), got (-1,1)
```

#### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
-->
1.3.1 1.17.3 sys.version_info(major=3, minor=8, micro=0, releaselevel='final', serial=0)

Hints from issue discussion:
According to the documentation, `interval`, produces a

> Confidence interval with equal areas around the median.

Here, presumably, "equal areas" means "equal mass".

If you use the definition that a *median* of a random variable $X$ is any number $m$ satisfying $P(X \geq m) \geq \frac{1}{2}$, $P(X \leq m) \geq \frac{1}{2}$, then the medians of `rv` are all numbers $0 < m < 1$. Then the above definition of course depends on the value of $m$, but in this particular case, any interval $[x, y]$ with $x \leq 0$, $y \geq 1$ will work for any median, and in particular $[-1, 1]$ is one such interval (as is $[0, 1]$); why would you expect $[0, 1]$ instead? It does have one nice property of being the intersection of all intervals with the given property, which could be useful if one were to create a definition of `interval`'s behavior that's well-defined.

Now for some reason `rv.median()` produces 0, and with that as a definition, *no* interval satisfies the condition in the documentation of `interval`, so following that, nothing should really be returned at all. But maybe then (without having consulted the implementation), $[-1, 1]$ is given since it happens to be symmetric around 0, and some weird notion of "area" is at play. That doesn't explain why `bernoulli.interval(0.999, 0.5)` produces `(0.0, 1.0)` though.

This is really an issue of the documentation being unclear, both in the case of `interval` and `median`.
`interval` is defined in terms of tail probabilities with at least the given coverage.

In the discrete case, if we need to include the lower bound, then the value is set to one lower (-1).

I always had problems with this, it's not intuitive and I don't remember it unless I go through the argument again.

The main reason is that `cdf` is defined as closed on the right, and `sf` as open on the left
i.e. cdf: Prob(x <= 0) > 0, or sf: Prob(x > 0) < 1
but cdf: Prob(x <= -1) == 0, or sf: Prob(x > -1) == 1

I don't think we can change this, because current computation with discrete distributions depend on this definition. AFAIR.


@josef-pkt Are you saying that the intervals produced by `interval` are open? Then its outputs are even more wrong in this case (so if other parts of `scipy.stats` depend on them, that should also be problematic). Or perhaps include the upper bound but not the lower? In any case, that really should be spelled out in the docs as well.
AFAIR, cdf and similar for discrete distributions are used in hypothesis tests for discrete data.
Those are verified in their results, so they use the current behavior correctly.

The main point in the definition is `Prob(x in A) >= coverage` rate

I don't remember or I'm not sure whether the interval is open on the left and closed on the right
i.e. Prob(x in A) = Prob(a < x <= b) = cdf(b) - cdf(a)

My guess would be that the usual definition for confidence intervals uses closed on both sides, based on the computation in related usage, i.e Prob(x in A) = cdf(b) - cdf(a - 1)

I haven't looked at this in the long time, so I don't remember details
`interval` was added later to the distributions, and I worked more with cdf and sf, e.g.
https://github.com/scipy/scipy/issues/2384

a quick check in statsmodels, for binom_test version I also only use cdf and sf.
I use `interval` only for continuous distributions.

So `interval` might not actually be used in followup code. 

In distributions like poisson, the interval can be in the interior for the sample space, That can be used to check the definition for other cases than coverage=1

I'm not sure interval for discrete distributions has been systematically reviewed, I don't think I ever did.

searching issues for "interval discrete"

#9706 docstring for interval
#10317 nbinom, poisson example seems to return 0 as lower endpoint, and not -1
#9762 (closed) poisson seems to include endpoints, argument Prob(x in A) = cdf(b) - cdf(a - 1)

----

There is a BUG for prob=1, I guess in some special case handling
```
stats.bernoulli.interval(1-1e-14, 0.5)
(0.0, 1.0)
```

same for poisson
```
stats.poisson.interval(1, 0.5)
(-1.0, inf)
stats.poisson.interval(1, 0)
(-1.0, inf)
stats.poisson.interval(0.999, 0.001)
(0.0, 1.0)
```
So, `interval` follows standard definition with closed interval on both sides

for discrete distributions 
`A = {x: a <= x <= b)`
`Prob(x in A) = cdf(b) - cdf(a - 1) >= coverage rate`


I didn't check the docs, but they should specify that in the discrete case both endpoints are in the interval

IIRC, then `interval` also requires  (equal tail)

alpha = 1 - coverage prob

lower tail prob
`cdf(a - 1) <= alpha / 2`
upper tail prob
`sf(b) <= alpha / 2`
Again, though, if you require that both end points of the interval be in the support, you will most likely not even have an interval with equal support on both sides on a median.
@WarrenWeckesser @mckib2 might this be another one for Boostinator?
@mdhaber Boost agrees with SciPy in all the cases discussed above, namely:
```
median: 0.0
interval(1): (-1.0, 1.0)
interval(.999): (0.0, 1.0)
interval(1-1e-14): (0.0, 1.0)
```

Commit message: TST/DOC: stats: explain/check 100% interval for discrete distributions

Patch to review:
diff --git a/scipy/stats/_distn_infrastructure.py b/scipy/stats/_distn_infrastructure.py
index 57dc3fa43ac8..a105523101f0 100644
--- a/scipy/stats/_distn_infrastructure.py
+++ b/scipy/stats/_distn_infrastructure.py
@@ -1532,6 +1532,15 @@ def interval(self, confidence=None, *args, **kwds):
             end-points of range that contain ``100 * alpha %`` of the rv's
             possible values.
 
+        Notes
+        -----
+        This is typically implemented as ``ppf([p_tail, 1-p_tail])``, where
+        ``ppf`` is the inverse cumulative distribution function and
+        ``p_tail = (1-confidence)/2``. Suppose ``[c, d]`` is the support of a
+        discrete distribution; then ``ppf([0, 1]) == (c-1, d)``. Therefore,
+        when ``confidence==1`` and the distribution is discrete, the left end
+        of the interval will be beyond the support of the distribution.
+
         """
         # This function was originally written with parameter `alpha`, but
         # `alpha` is also the name of a shape parameter of two distributions.
diff --git a/scipy/stats/tests/test_discrete_basic.py b/scipy/stats/tests/test_discrete_basic.py
index 70d3d1017dea..d02451d29e27 100644
--- a/scipy/stats/tests/test_discrete_basic.py
+++ b/scipy/stats/tests/test_discrete_basic.py
@@ -384,3 +384,19 @@ def test_frozen_attributes():
     frozen_binom = stats.binom(10, 0.5)
     assert isinstance(frozen_binom, rv_discrete_frozen)
     delattr(stats.binom, 'pdf')
+
+
+@pytest.mark.parametrize('distname, shapes', distdiscrete)
+def test_interval(distname, shapes):
+    # gh-11026 reported that `interval` returns incorrect values when
+    # `confidence=1`. The values were not incorrect, but it was not intuitive
+    # that the left end of the interval should extend beyond the support of the
+    # distribution. Confirm that this is the behavior for all distributions.
+    if isinstance(distname, str):
+        dist = getattr(stats, distname)
+    else:
+        dist = distname
+    a, b = dist.support(*shapes)
+    npt.assert_equal(dist.ppf([0, 1], *shapes), (a-1, b))
+    npt.assert_equal(dist.isf([1, 0], *shapes), (a-1, b))
+    npt.assert_equal(dist.interval(1, *shapes), (a-1, b))


Code review comments:
diff_hunk: @@ -1532,6 +1532,15 @@ def interval(self, confidence=None, *args, **kwds):
             end-points of range that contain ``100 * alpha %`` of the rv's
             possible values.
 
+        Notes
+        -----
+        This is typically implemented as ``ppf([p_tail, 1-p_tail])``, where
path: scipy/stats/_distn_infrastructure.py
text: @user1:
drop "typically"
it is implemented


@author:
Good idea. I meant that a distribution could override the generic implementation, but I don't think any do.