[issue37439] Add random.binomialvariate()

[Mark Dickinson]

Mark Dickinson <dickinsm at gmail.com> added the comment:

If you want to be able to switch to something more efficient for large `n`, Knuth describes a simple O(log(n)) algorithm in TAOCP volume 4 (and attributes it to J. H. Ahrens). It's essentially a bisection search on the value, using the fact that we can use the beta distribution to generate order statistics. Here's a (too) simple Python implementation. It still needs thorough testing, and could be optimised in many ways - e.g., using sensible crossover point for n and not recursing all the way to n = 0.


    def binomialvariate(n, p):
        if n == 0:
            return 0
        a, b = (1 + n)//2, 1 + n//2
        x = betavariate(a, b)
        if x < p:
            return a + binomialvariate(b - 1, (p - x)/(1 - x))
        else:
            return binomialvariate(a - 1, p/x)


>>> binomialvariate(10**10, 0.5)
4999944649

----------

_______________________________________
Python tracker <report at bugs.python.org>
<https://bugs.python.org/issue37439>
_______________________________________