Numerical algorithm to generate numbers from Binomial distribution

Question

I need to generate random numbers from Binomial(n, p) distribution.

A Binomial(n, p) random variable is sum of n uniform variables which take 1 with probability p. In pseudo code, x=0; for(i=0; I<n; ++I) x+=(rand()<p?1:0); will generate a Binomial(n, p).

I need to generate this for small as well as really large n, for example n = 10^6 and p=0.02. Is there any fast numerical algorithm to generate it?

EDIT -

Right now this is what I have as approximation (along with functions for exact Poisson and Normal distribution)-

    public long Binomial(long n, double p) {
        // As of now it is an approximation
        if (n < 1000) {
            long result = 0;
            for (int i=0; i<n; ++i)
                if (random.NextDouble() < p) result++;
            return result;
        }
        if (n * p < 10) return Poisson(n * p);
        else if (n * (1 - p) < 10) return n - Poisson(n * p);
        else {
            long v = (long)(0.5 + nextNormal(n * p, Math.Sqrt(n * p * (1 - p))));
            if (v < 0) v = 0;
            else if (v > n) v = n;
            return v;
        }
    }

Does the Poisson distribution work well for `n * p < 10` or for `n * (1 - p) < 10`? How comes that you choose that distribution? — HelloGoodbye, Mar 06 '14 at 13:44
Yes, for large n. Binomial(n,lambda/n) converges to Poisson(lambda), as n goes to infinity. — KalEl, Mar 27 '14 at 03:31

score 5 · Answer 1 · answered Nov 19 '09 at 15:03

Another option would be to sample from Normal or Poisson as you do and then add a Metropolis-Hastings step to accept or reject your sample. If you accept you are done, if you reject, you have to completely resample again. My guess is that because the approximation is so close, you will almost always get an accept step, once in a while you might reject.

Also Luc Devroye's book has some great algorithms for Binomial sampling.

PS If you end up with a good algorithm; would you mind sharing it at Math.Net Numerics?

Richie Cotton · Accepted Answer · 2009-11-13T12:07:04.830

4

If you are willing to pay, then take a look at NMath by Centerspace.

Otherwise, the C code used by the Stats program R is here, and should be straightforward to port to C#.

EDIT: There are details (inc. code) on creating a method for this on p178 of Practical Numerical Methods with C# by Jack Xu.

ANOTHER EDIT: A free C# library that does what you want.

edited Nov 13 '09 at 12:07

answered Nov 13 '09 at 11:59

Richie Cotton

118,240
47
247
360

score 3 · Answer 3 · answered Nov 13 '09 at 11:53

3

There's no obvious way to do this efficiently. For small n, you might as well just us the formula to calculate the inverse PDF. For larger n, you're probably best off using one of the approximations to other distributions that are easier to calculate.

answered Nov 13 '09 at 11:53

Noldorin

144,213
56
264
302

I'm actually using some adhoc approximation routine I wrote that uses the Poisson and Normal distributions. However I'm not sure how far off the numbers it generates will statistically be from a true Binomial. I've added my code as an edit. – KalEl Nov 13 '09 at 12:19

Numerical algorithm to generate numbers from Binomial distribution

3 Answers3

Linked