Rather than using mathematical libraries, how would you sample from a binomial random variable efficiently?

Given the binomial random variable X, where $k$ are the number of successes in $n$ trials with a success probability of $p$ $$ P(X=k|n, p) = \binom{n}{k} p^k(1-p)^{n-k}, $$ how could I obtain $N$ number of samples of $X$?

The naïve approach is to decompose $X$ into $X = Y_1 + Y_2 + ... + Y_n$, where $Y$ are the bernoulli experiments $$ Y_i = P(Y_i = 1) = p. $$

Or, in other words, I can test $n$ times if some random value is above $p$ and count how many times this was a success. This is, however, terribly inefficient when $n$ is large.

It is possible to sample a continuous random variable by finding the inverse CDF ($F^{-1}(x)$), sampling from the uniform distribution $ u = U(0,1)$ and calculating the value of the sample in the inverse CDF $F^{-1}(u)$. Given that this is a discrete distribution, how could I apply this? Which other methods are available?

  • $\begingroup$ It seems you mixed up $N$ and $n$ a couple of times? $\endgroup$ – joriki Sep 8 '15 at 23:19
  • $\begingroup$ Yes, I mixed them up. Question edited to correct these. $\endgroup$ – UndeadKernel Sep 9 '15 at 9:25

It makes no difference if the distribution from which you wish to sample is continuous or discrete. For example, suppose you wish to sample from $$X \sim \operatorname{Binomial}(n = 5, p = 0.7).$$ Then $$\begin{align*} \Pr[X < 0] &= 0 \\ \Pr[X \le 0] &= 0.00243 \\ \Pr[X \le 1] &= 0.03078 \\ \Pr[X \le 2] &= 0.16308 \\ \Pr[X \le 3] &= 0.47178 \\ \Pr[X \le 4] &= 0.83193 \\ \Pr[X \le 5] &= 1. \end{align*}$$ Then you draw realizations $u_i$ from a continuous uniform distribution on $[0,1]$, and you will select $X_i = x_i$ if $$\Pr[X < x_i] < u_i \le \Pr[X \le x_i].$$ This is the inverse CDF method. You'd create a lookup table for a given $n$ and $p$, and locate the index of the interval in which your random number lies between.

I can't speak to which method is faster, but I can say for certain that the inverse CDF method for large $n$ requires a lot more memory to create the lookup table, and is probably not as efficient. By the time your $n$ is so large that creating individual realizations is not computationally feasible, you are likely to be better off using a normal approximation to the binomial (unless $p$ is extremely close to $0$ or $1$).

  • $\begingroup$ I suppose you meant: $Pr[X<x_i]<ui≤Pr[X≤x_{i+1}]$. Notice $x_{i+1}$ in the second inequality. $\endgroup$ – UndeadKernel Sep 9 '15 at 9:27
  • $\begingroup$ @UndeadKernel No, I did mean what I typed. Note that your definition does not permit the choice of $X_i = 5$, for instance, if $u_i$ is to be drawn from $[0,1]$, since $\Pr[U = 1] = 0$. My definition works nicely because it is more naturally extensible to the continuous case; as then $\Pr[X < x_i] = u_i = F(x_i)$ hence $x_i = F^{-1}(u_i)$. $\endgroup$ – heropup Sep 9 '15 at 11:04
  • 1
    $\begingroup$ @UndeadKernel Cf. math.stackexchange.com/questions/1420248/… $\endgroup$ – Ian Sep 9 '15 at 11:21

A method for the general discrete case is called the alias method. This performs a preprocessing procedure, which can be done in $O(n)$ time*. This converts the $n$ values and probabilities into $n$ "bins", each of which contains either $1$ or $2$ values and some of their probabilities. The construction ensures that each bin has equal total probability. Thus, to sample, you uniformly choose a bin. Then you sample from the two values using their relative probabilities within the bin.

Using the alias method, the sampling time is almost entirely accounted for by the time to generate a sufficiently precise uniform random variable in order to resolve the probabilities within the bins. This is pretty fast; for most practical purposes it can be considered constant time. In particular, it is faster than the $\log(n)$ time that would be required for binary search with a CDF lookup table.

* Note that a naive implementation of the preprocessor chooses the value with the smallest probability (say $m$) to be one half of the bin and the value with the largest probability to be the maximum probability (say $M$) to be the largest one. This works well for theory because it makes it more obvious that you can always continue the procedure, but it slows down the implementation because you have to figure out where to put $M$ now that you've subtracted probability from it. Doing that bumps the complexity of the preprocessor up to $O(n\log(n))$. But this is not necessary: all that matters is that the values with probability less than $1/n$ are in one place while the values with probability more than $1/n$ are in another. Thus you do need to find a place to put $M$, but it is easy to do so, just check whether the new value of $P(M)$ is less than $1/n$ or greater than $1/n$.

Reference: http://www.jstor.org/stable/2683739 (many other references can be found).

  • $\begingroup$ Very nice method; I'd never heard of it. Here's a concise explanation and implementation in Python. $\endgroup$ – joriki Sep 9 '15 at 0:48

You could use the acceptance rejection method when $n$ is large

A pseudo code is as below

Define $c_{max} = \underset{i}{max} (n+1) {n \choose i} p^i (1-p)^{n-i} $

  1. Generate $U_1 \sim U[0,1]$. Let $I = Int((n+1)*U_1)+1$

  2. Generate $U_2 \sim U[0,1]$.

  3. If $U_2 < {n \choose I-1} p^{I-1} (1-p)^{n-I+1}/c_{max}$ output $I-1$ else go to step 1

  • 2
    $\begingroup$ I didn't understand what you were doing until I realized this was the Metropolis algorithm, with the initial transition matrix having entries $\frac{1}{n+1}$ everywhere. $\endgroup$ – Ian Sep 9 '15 at 10:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.