scipy.stats.binom_test / binom.sf return incorrect values for large x and n

[p00ya]

scipy.stats.binom.sf and scipy.stats.binom_test returns incorrect values for large inputs. The bounds on the inputs is not documented, and similar functions in other libraries (e.g. R's binom.test) do not have this problem.

Reproducing code example:

>>> np.seterr(all='warn')
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}
>>> scipy.stats.binom_test(1e7, 2e7, 0.5, 'greater')
0.5094024848633918
>>> scipy.stats.binom_test(1e8, 2e8, 0.5, 'greater')
0.6860874231168742
>>> scipy.stats.binom_test(1e9, 2e9, 0.5, 'greater')
0.8854076875623922
>>> scipy.stats.binom_test(1e10, 2e10, 0.5, 'greater')
nan

Error message:

None.

Expected behaviour:

All calls in the repro example should return a value near 0.5. For example, in R:

> binom.test(1e7, 2e7, 0.5, 'greater')$p.value
[1] 0.5000892
> binom.test(1e8, 2e8, 0.5, 'greater')$p.value
[1] 0.5000282
> binom.test(1e9, 2e9, 0.5, 'greater')$p.value
[1] 0.5000089
> binom.test(1e10, 2e10, 0.5, 'greater')$p.value
[1] 0.5000028

Notes:

Internally, for the "greater" alternative, scipy calls binom.sf(x - 1, n, p), which then calls scipy.special.bdtrc(floor(x), n, p). This forwards to Cephes' bdtrc, according to scipy's docs.

We can verify that the bug lies within bdtrc:

>>> scipy.special.bdtrc(1e10, 2e10 - 1e10 + 1, 0.5)
nan

Non-broken implementations of the regularized incomplete beta function (such as TensorFlow's tf.math.betainc) will return the expected value (0.5).

This also affects anything else that calls bdtrc, including binom.sf. There may be a case that bdtrc is Working as Intended (because the documentation is explicit that it's just a wrapper for Cephes' buggy implementation), but I think for binom.sf and binom_test it's clear that there is at least a documentation bug if not an opportunity to make things better.

I think there are 4 viable solutions, in what I think is most-preferred to least-preferred:

1. File a bug with Cephes and wait for them to fix it, documenting the breakage in the meantime.
2. Use a normal approximation in binom.sf when "k" and "n" are large (perhaps checking that "p" isn't too biased). Fortunately the circumstances when bdtrc breaks usually coincide with when this approximation is good. :)
3. Implement a regularized incomplete beta function within scipy that isn't broken, instead of using Cephes.
4. Claim this is Working as Intended, but document the limitation.

I'm happy to send Pull Requests given some direction about which of the choices above to go with.

Scipy/Numpy/Python version information:

>>> 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)

[chrisb83]

thanks for reporting it, this is bug which also impacts the distribution stats.binom

stats.binom.cdf(1e9, 2e9, 0.5) # 0.885...

[mdhaber]

Sounds like a potential job for boostinator!
@mckib2 Could you add this to the demo?

[mckib2]

@mdhaber Done. The boost implementation of binomial survival function appears to do the right thing here.

[rlucas7]

I think this is the same issue as this one: #5503

…

-Lucas Roberts

On Nov 15, 2020, at 2:32 PM, Christoph Baumgarten ***@***.***> wrote:  thanks for reporting it, this is bug which also impacts the distribution stats.binom stats.binom.cdf(1e9, 2e9, 0.5) # 0.885... — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[mdhaber]

I agree with @rlucas7 that this is essentially a duplicate of gh-5503, but I'll add a note there to look for additional test cases here.