8
$\begingroup$

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?

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

3 Answers 3

5
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ I suppose you meant: $Pr[X<x_i]<ui≤Pr[X≤x_{i+1}]$. Notice $x_{i+1}$ in the second inequality. $\endgroup$ Commented Sep 9, 2015 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
    Commented Sep 9, 2015 at 11:04
  • 1
    $\begingroup$ @UndeadKernel Cf. math.stackexchange.com/questions/1420248/… $\endgroup$
    – Ian
    Commented Sep 9, 2015 at 11:21
  • $\begingroup$ I don't think this method is actually better than generating Bernoullis and summing them up for large n. $\endgroup$
    – mhsnk
    Commented Dec 8, 2021 at 1:31
4
$\begingroup$

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 are either done (if the bin has 1 value) or you independently choose one of the 2 values in the bin 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).

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

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

$\endgroup$
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
    Commented Sep 9, 2015 at 10:35

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .