Repository: statsmodels/statsmodels (Python)

Issue:
proportion_confint inaccurate for method="binom_test"
There is a symmetry for binomial proportion confidence interval therefore following property should holds (approximately):
proportion_confint(count, nobs, method="binom_test")[0] + proportion_confint(nobs-count, nobs, method="binom_test")[1] ==1.0
but for some cases this property does not hold, for example, the following code:
```python
from statsmodels.stats.proportion import proportion_confint
eps = 1e-8
for nobs in [47, 50]:
    for count in range(nobs+1):
        p1 = proportion_confint(count, nobs, method="binom_test")
        p2 = proportion_confint(nobs-count, nobs, method="binom_test")
        if abs(p1[0]+p2[1]-1.0)>eps:
            print(f"{nobs=} {count=} {p1[0]+p2[1]}")
```
gives ouput:
```
nobs=47 count=3 0.9899913901078568
nobs=47 count=6 0.9857792597568958
nobs=47 count=12 0.9864319871901812
nobs=47 count=16 0.9851553703015263
nobs=47 count=23 0.9802691065584829
nobs=47 count=24 0.980269106558583
nobs=47 count=31 0.978547082683319
nobs=47 count=35 0.9718877308702878
nobs=47 count=41 0.978042999703059
nobs=47 count=44 0.9680701802192785
nobs=50 count=21 0.9796386543796607
nobs=50 count=29 0.9798450168011505
```
Which indicate that there is some bug.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``
<details>

INSTALLED VERSIONS
------------------
Python: 3.8.5.final.0

statsmodels
===========

Installed: 0.13.1 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 1.17.3 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\numpy)
scipy: 1.3.3 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\scipy)
pandas: 0.25.3 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\pandas)
    dateutil: 2.8.1 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\dateutil)
patsy: 0.5.2 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\patsy)

Optional Dependencies
=====================

matplotlib: 3.2.0rc1 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\matplotlib)
    backend: TkAgg
cvxopt: Not installed
joblib: Not installed

Developer Tools
================

IPython: 7.17.0 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\IPython)
    jinja2: 2.11.2 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\jinja2)
sphinx: Not installed
    pygments: 2.6.1 (C:\Users\Admin\AppData\Local\Programs\Python\Python38-32\lib\site-packages\pygments)
pytest: Not installed
virtualenv: Not installed

</details>


Hints:
Thanks for reporting

I can confirm this, but don't see any obvious problem.
It's either because of some discreteness in scipy.stats.binom_test or in scipy.optimize.brentq.
The differences in the confint are nontrivial, but inverting binom_test is not one of the recommended confints (both in the literature and in terms of the quickly done implementation).

But I don't have time to dig into this right now (It's too far away from what I'm currently working on.)

label `bug-wrong` because we should have consistent results for symmetry.

However my guess is that they choose different solution because of non-monotonicities or discreteness, which are often a problem with hypothesis tests for proportions.
Requirement for size is only `<= alpha` but not `== alpha` because of the discreteness of the sample space.

Patch:
diff --git a/statsmodels/distributions/copula/transforms.py b/statsmodels/distributions/copula/transforms.py
--- a/statsmodels/distributions/copula/transforms.py
+++ b/statsmodels/distributions/copula/transforms.py
@@ -8,6 +8,7 @@
 License: BSD-3
 
 """
+import warnings
 
 import numpy as np
 from scipy.special import expm1
@@ -24,7 +25,10 @@ class TransfFrank(Transforms):
 
     def evaluate(self, t, theta):
         t = np.asarray(t)
-        return - (np.log(-expm1(-theta*t)) - np.log(-expm1(-theta)))
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", RuntimeWarning)
+            val = -(np.log(-expm1(-theta*t)) - np.log(-expm1(-theta)))
+        return val
         # return - np.log(expm1(-theta*t) / expm1(-theta))
 
     def inverse(self, phi, theta):
diff --git a/statsmodels/distributions/tools.py b/statsmodels/distributions/tools.py
--- a/statsmodels/distributions/tools.py
+++ b/statsmodels/distributions/tools.py
@@ -6,6 +6,7 @@
 License: BSD-3
 
 """
+import warnings
 
 import numpy as np
 from scipy import interpolate, stats
@@ -324,13 +325,19 @@ def _eval_bernstein_1d(x, fvals, method="binom"):
     n = k_terms - 1.
 
     if method.lower() == "binom":
-        poly_base = stats.binom.pmf(k, n, xx[..., None])
+        # Divide by 0 RuntimeWarning here
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", RuntimeWarning)
+            poly_base = stats.binom.pmf(k, n, xx[..., None])
         bp_values = (fvals * poly_base).sum(-1)
     elif method.lower() == "bpoly":
         bpb = interpolate.BPoly(fvals[:, None], [0., 1])
         bp_values = bpb(x)
     elif method.lower() == "beta":
-        poly_base = stats.beta.pdf(xx[..., None], k + 1, n - k + 1) / (n + 1)
+        # Divide by 0 RuntimeWarning here
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", RuntimeWarning)
+            poly_base = stats.beta.pdf(xx[..., None], k + 1, n - k + 1) / (n + 1)
         bp_values = (fvals * poly_base).sum(-1)
     else:
         raise ValueError("method not recogized")
diff --git a/statsmodels/stats/proportion.py b/statsmodels/stats/proportion.py
--- a/statsmodels/stats/proportion.py
+++ b/statsmodels/stats/proportion.py
@@ -8,19 +8,21 @@
 """
 
 from statsmodels.compat.python import lzip
+
 import numpy as np
-from scipy import stats, optimize
-from sys import float_info
+from scipy import optimize, stats
 
-from statsmodels.stats.base import AllPairsResults
-from statsmodels.tools.sm_exceptions import HypothesisTestWarning
+from statsmodels.stats.base import AllPairsResults, HolderTuple
 from statsmodels.stats.weightstats import _zstat_generic2
-from statsmodels.stats.base import HolderTuple
+from statsmodels.tools.sm_exceptions import HypothesisTestWarning
 from statsmodels.tools.testing import Holder
 
+FLOAT_INFO = np.finfo(float)
+
 
 def proportion_confint(count, nobs, alpha=0.05, method='normal'):
-    '''confidence interval for a binomial proportion
+    """
+    Confidence interval for a binomial proportion
 
     Parameters
     ----------
@@ -76,17 +78,28 @@ def proportion_confint(count, nobs, alpha=0.05, method='normal'):
     Brown, Lawrence D.; Cai, T. Tony; DasGupta, Anirban (2001). "Interval
         Estimation for a Binomial Proportion",
         Statistical Science 16 (2): 101–133. doi:10.1214/ss/1009213286.
-
-    '''
+    """
 
     pd_index = getattr(count, 'index', None)
     if pd_index is not None and callable(pd_index):
         # this rules out lists, lists have an index method
         pd_index = None
-    count = np.asarray(count)
-    nobs = np.asarray(nobs)
 
-    q_ = count * 1. / nobs
+    def _check(x: np.ndarray, name: str) -> np.ndarray:
+        if np.issubdtype(x.dtype, np.integer):
+            return x
+        y = x.astype(np.int64, casting="unsafe")
+        if np.any(y != x):
+            raise ValueError(
+                f"{name} must have an integral dtype. Found data with "
+                f"dtype {x.dtype}"
+            )
+        return y
+
+    count = _check(np.asarray(count), "count")
+    nobs = _check(np.asarray(nobs), "nobs")
+
+    q_ = count / nobs
     alpha_2 = 0.5 * alpha
 
     if method == 'normal':
@@ -95,18 +108,53 @@ def proportion_confint(count, nobs, alpha=0.05, method='normal'):
         ci_low = q_ - dist
         ci_upp = q_ + dist
 
+
+
     elif method == 'binom_test':
         # inverting the binomial test
         def func(qi):
-            return stats.binom_test(q_ * nobs, nobs, p=qi) - alpha
+            return stats.binom_test(count, nobs, p=qi) - alpha
+
+        def _bound(qi, lower=True):
+            """
+            Try hard to find a bound different from eps or 1 - eps
+
+            Parameters
+            ----------
+            qi : float
+                The empirical success rate
+            lower : bool
+                Whether to fund a lower bound for the left side of the CI
+
+            Returns
+            -------
+            float
+                The coarse bound
+            """
+            default = FLOAT_INFO.eps if lower else 1.0 - FLOAT_INFO.eps
+
+            def step(v):
+                return v / 8 if lower else v + (1.0 - v) / 8
+
+            x = step(qi)
+            w = func(x)
+            cnt = 1
+            while w > 0 and cnt < 10:
+                x = step(x)
+                w = func(x)
+                cnt += 1
+            return x if cnt < 10 else default
+
         if count == 0:
             ci_low = 0
         else:
-            ci_low = optimize.brentq(func, float_info.min, q_)
+            lower_bnd = _bound(q_, lower=True)
+            ci_low = optimize.brentq(func, lower_bnd, q_)
         if count == nobs:
             ci_upp = 1
         else:
-            ci_upp = optimize.brentq(func, q_, 1. - float_info.epsilon)
+            upper_bnd = _bound(q_, lower=False)
+            ci_upp = optimize.brentq(func, q_, upper_bnd)
 
     elif method == 'beta':
         ci_low = stats.beta.ppf(alpha_2, count, nobs - count + 1)


Test patch:
diff --git a/statsmodels/distributions/tests/test_bernstein.py b/statsmodels/distributions/tests/test_bernstein.py
--- a/statsmodels/distributions/tests/test_bernstein.py
+++ b/statsmodels/distributions/tests/test_bernstein.py
@@ -21,6 +21,9 @@
 
 def test_bernstein_distribution_1d():
     grid = dt._Grid([501])
+    loc = grid.x_flat == 0
+    grid.x_flat[loc] = grid.x_flat[~loc].min() / 2
+    grid.x_flat[grid.x_flat == 1] = 1 - grid.x_flat.min()
 
     distr = stats.beta(3, 5)
     cdf_g = distr.cdf(np.squeeze(grid.x_flat))
diff --git a/statsmodels/stats/tests/test_proportion.py b/statsmodels/stats/tests/test_proportion.py
--- a/statsmodels/stats/tests/test_proportion.py
+++ b/statsmodels/stats/tests/test_proportion.py
@@ -7,24 +7,29 @@
 """
 import warnings
 
-import pytest
 import numpy as np
+from numpy.testing import (
+    assert_allclose,
+    assert_almost_equal,
+    assert_array_less,
+    assert_equal,
+    assert_raises,
+)
 import pandas as pd
-from numpy.testing import (assert_almost_equal, assert_equal,
-                           assert_array_less, assert_raises, assert_allclose)
-
-from statsmodels.stats.proportion import (proportion_confint,
-                                          confint_proportions_2indep,
-                                          multinomial_proportions_confint,
-                                          score_test_proportions_2indep,
-                                          power_proportions_2indep,
-                                          samplesize_proportions_2indep_onetail,
-                                          )
+import pytest
+
 import statsmodels.stats.proportion as smprop
+from statsmodels.stats.proportion import (
+    confint_proportions_2indep,
+    multinomial_proportions_confint,
+    power_proportions_2indep,
+    proportion_confint,
+    samplesize_proportions_2indep_onetail,
+    score_test_proportions_2indep,
+)
 from statsmodels.tools.sm_exceptions import HypothesisTestWarning
 from statsmodels.tools.testing import Holder
 
-
 probci_methods = {'agresti_coull': 'agresti-coull',
                   'normal': 'asymptotic',
                   'beta': 'exact',
@@ -909,3 +914,34 @@ def test_power_2indep():
                                                ratio=1, alpha=0.05, value=0,
                                                alternative='two-sided')
     assert_allclose(n2, n, rtol=1e-13)
+
+
+@pytest.mark.parametrize("count", np.arange(10, 90, 5))
+@pytest.mark.parametrize("method", list(probci_methods.keys()) + ["binom_test"])
+def test_ci_symmetry(count, method):
+    n = 100
+    a = proportion_confint(count, n, method=method)
+    b = proportion_confint(n - count, n, method=method)
+    assert_allclose(np.array(a), 1.0 - np.array(b[::-1]))
+
+
+@pytest.mark.parametrize("nobs", [47, 50])
+@pytest.mark.parametrize("count", np.arange(48))
+def test_ci_symmetry_binom_test(nobs, count):
+    a = proportion_confint(count, nobs, method="binom_test")
+    b = proportion_confint(nobs - count, nobs, method="binom_test")
+    assert_allclose(np.array(a), 1.0 - np.array(b[::-1]))
+
+
+def test_int_check():
+    with pytest.raises(ValueError):
+        proportion_confint(10.5, 20)
+    with pytest.raises(ValueError):
+        proportion_confint(10, 20.5)
+    with pytest.raises(ValueError):
+        proportion_confint(np.array([10.3]), 20)
+    a = proportion_confint(21.0, 47, method="binom_test")
+    b = proportion_confint(21, 47, method="binom_test")
+    c = proportion_confint(21, 47.0, method="binom_test")
+    assert_allclose(a, b)
+    assert_allclose(a, c)