<p dir="auto"><a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom.html#scipy.stats.binom" rel="nofollow">scipy.stats.binom.sf</a> and <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom_test.html" rel="nofollow">scipy.stats.binom_test</a> returns incorrect values for large inputs.  The bounds on the inputs is not documented, and similar functions in other libraries (e.g. R's <code class="notranslate">binom.test</code>) do not have this problem.</p>
<h4 dir="auto">Reproducing code example:</h4>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; np.seterr(all='warn')
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}
&gt;&gt;&gt; scipy.stats.binom_test(1e7, 2e7, 0.5, 'greater')
0.5094024848633918
&gt;&gt;&gt; scipy.stats.binom_test(1e8, 2e8, 0.5, 'greater')
0.6860874231168742
&gt;&gt;&gt; scipy.stats.binom_test(1e9, 2e9, 0.5, 'greater')
0.8854076875623922
&gt;&gt;&gt; scipy.stats.binom_test(1e10, 2e10, 0.5, 'greater')
nan"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; np.seterr(all='warn')
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}
&gt;&gt;&gt; scipy.stats.binom_test(1e7, 2e7, 0.5, 'greater')
0.5094024848633918
&gt;&gt;&gt; scipy.stats.binom_test(1e8, 2e8, 0.5, 'greater')
0.6860874231168742
&gt;&gt;&gt; scipy.stats.binom_test(1e9, 2e9, 0.5, 'greater')
0.8854076875623922
&gt;&gt;&gt; scipy.stats.binom_test(1e10, 2e10, 0.5, 'greater')
nan
</code></pre></div>
<h4 dir="auto">Error message:</h4>
<p dir="auto">None.</p>
<h4 dir="auto">Expected behaviour:</h4>
<p dir="auto">All calls in the repro example should return a value near 0.5.  For example, in R:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt; binom.test(1e7, 2e7, 0.5, 'greater')$p.value
[1] 0.5000892
&gt; binom.test(1e8, 2e8, 0.5, 'greater')$p.value
[1] 0.5000282
&gt; binom.test(1e9, 2e9, 0.5, 'greater')$p.value
[1] 0.5000089
&gt; binom.test(1e10, 2e10, 0.5, 'greater')$p.value
[1] 0.5000028"><pre class="notranslate"><code class="notranslate">&gt; binom.test(1e7, 2e7, 0.5, 'greater')$p.value
[1] 0.5000892
&gt; binom.test(1e8, 2e8, 0.5, 'greater')$p.value
[1] 0.5000282
&gt; binom.test(1e9, 2e9, 0.5, 'greater')$p.value
[1] 0.5000089
&gt; binom.test(1e10, 2e10, 0.5, 'greater')$p.value
[1] 0.5000028
</code></pre></div>
<h4 dir="auto">Notes:</h4>
<p dir="auto">Internally, for the "greater" alternative, scipy calls <code class="notranslate">binom.sf(x - 1, n, p)</code>, which then calls <code class="notranslate">scipy.special.bdtrc(floor(x), n, p)</code>.  This forwards to <a href="http://www.netlib.org/cephes/doubldoc.html#bdtrc" rel="nofollow">Cephes' bdtrc</a>, according to scipy's docs.</p>
<p dir="auto">We can verify that the bug lies within <code class="notranslate">bdtrc</code>:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; scipy.special.bdtrc(1e10, 2e10 - 1e10 + 1, 0.5)
nan"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; scipy.special.bdtrc(1e10, 2e10 - 1e10 + 1, 0.5)
nan
</code></pre></div>
<p dir="auto">Non-broken implementations of the regularized incomplete beta function (such as TensorFlow's <code class="notranslate">tf.math.betainc</code>) will return the expected value (0.5).</p>
<p dir="auto">This also affects anything else that calls <code class="notranslate">bdtrc</code>, including <code class="notranslate">binom.sf</code>.  There may be a case that <code class="notranslate">bdtrc</code> is Working as Intended (because the documentation is explicit that it's just a wrapper for Cephes' buggy implementation), but I think for <code class="notranslate">binom.sf</code> and <code class="notranslate">binom_test</code> it's clear that there is at least a documentation bug if not an opportunity to make things better.</p>
<p dir="auto">I think there are 4 viable solutions, in what I think is most-preferred to least-preferred:</p>
<ol dir="auto">
<li>File a bug with Cephes and wait for them to fix it, documenting the breakage in the meantime.</li>
<li>Use a normal approximation in <code class="notranslate">binom.sf</code> when "k" and "n" are large (perhaps checking that "p" isn't too biased).  Fortunately the circumstances when <code class="notranslate">bdtrc</code> breaks usually coincide with when this approximation is good. :)</li>
<li>Implement a regularized incomplete beta function within scipy that isn't broken, instead of using Cephes.</li>
<li>Claim this is Working as Intended, but document the limitation.</li>
</ol>
<p dir="auto">I'm happy to send Pull Requests given some direction about which of the choices above to go with.</p>
<h4 dir="auto">Scipy/Numpy/Python version information:</h4>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
1.5.4 1.19.4 sys.version_info(major=3, minor=6, micro=8, releaselevel='final', serial=0)"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
1.5.4 1.19.4 sys.version_info(major=3, minor=6, micro=8, releaselevel='final', serial=0)
</code></pre></div>
<p dir="auto">The following command correctly returns around 0.5.<br>
scipy.stats.binom.cdf(3000000, 10000000, .3)</p>
<p dir="auto">The probability that a binomial RV with p=0.3 is less than or equal to 3,000,000 in 10,000,000 draws is around 0.5.</p>
<p dir="auto">The same should be true for P(X&lt;=30,000,000) in n=100,000,000 draws. (these numbers are 10 times higher than the preceding inputs)<br>
scipy.stats.binom.cdf(30000000, 100000000, .3)<br>
should return around 0.5, but it returns 0.57036.</p>
<p dir="auto">I also tested scipy.special.betainc, since that's what is used for calculating Binomial's CDF, and it exhibits the same problem.</p>