Binomial Sampling By Discrete Inversion

Why invert a binomial distribution?

A binomial count records the number of successes in a fixed number of independent trials, each with the same success probability. One way to simulate it is to generate a Bernoulli outcome for every trial and add the results. Discrete inversion offers another route: generate one uniform variate and use the binomial cumulative probabilities to select the count. This method is especially useful when the probability masses can be computed sequentially, without evaluating a separate binomial coefficient for every possible count.

The essential task is to locate the first cumulative probability that reaches the uniform variate. The probability mass recurrence makes this search practical: once the mass at zero is known, each subsequent mass follows from the previous one by a simple ratio. The method must still respect the discrete inverse convention; selecting the wrong side of a cumulative threshold can shift the simulated outcome.

Definition: If \(X\sim\operatorname{Bin}(n,p)\), then \(X\) takes values \(0,1,\ldots,n\), with \(P(X=k)=\binom{n}{k}p^k(1-p)^{n-k}\). Write \(q=1-p\). For a uniform variate \(U\) on \((0,1)\), discrete inversion returns the smallest integer \(k\) for which \(P(X\leq k)\geq U\).

Recursive probability masses

For \(0<p<1\), the probability of zero successes is \(P(X=0)=q^n\). For \(k=0,\ldots,n-1\), adjacent binomial masses satisfy

P(X=k+1)=P(X=k)\frac{n-k}{k+1}\frac{p}{q}.

This relation lets us compute the masses in order. Define \(a_k=P(X=k)\) and \(C_k=P(X\leq k)\). Set \(a_0=q^n\), then calculate \(a_{k+1}\) from \(a_k\) and update the cumulative total by \(C_{k+1}=C_k+a_{k+1}\). Stop when \(C_k\geq U\); the selected count is \(k\).

Derivation: The binomial probability formula gives \(a_{k+1}=\binom{n}{k+1}p^{k+1}q^{n-k-1}\) and \(a_k=\binom{n}{k}p^kq^{n-k}\). Dividing the first expression by the second yields \[ \frac{a_{k+1}}{a_k} =\frac{\binom{n}{k+1}}{\binom{n}{k}}\frac{p}{q} =\frac{n-k}{k+1}\frac{p}{q}. \] Multiplying by \(a_k\) gives the recurrence. The ratio of binomial coefficients follows directly from their factorial definitions.

The cumulative sum is updated one mass at a time, so there is no need to calculate every mass if the threshold is reached early. A compact procedure is:

Check the parameters.
Use the specified integer \(n\) and probability \(p\). If \(p=0\), return 0; if \(p=1\), return \(n\).

Initialize.
For \(0<p<1\), set \(a_0=q^n\), \(C_0=a_0\), and \(k=0\).

Compare and advance.
If \(C_k\geq U\), return \(k\). Otherwise calculate \(a_{k+1}=a_k\frac{n-k}{k+1}\frac{p}{q}\), set \(C_{k+1}=C_k+a_{k+1}\), and increase \(k\) by one.

Finish by \(n\).
The masses sum to one, so the cumulative total at \(n\) is one. A valid \(U\in(0,1)\) must therefore select some value no greater than \(n\).

The stopping comparison is part of the method, not a coding detail. If \(U\) equals a cumulative probability exactly, the smallest \(k\) with \(C_k\geq U\) is selected. For an ideal continuous uniform variate, landing exactly on a particular threshold has probability zero, but a computer's finite set of possible values can make boundary handling relevant.

Worked applications

Worked Example: A small success probability

A hypothetical loss-prevention test is applied independently to \(n=5\) devices, with probability \(p=0.20\) of detecting a fault on each device. Use discrete inversion with \(U=0.70\) to simulate the number of detections. The uniform variate is supplied as part of the problem; assume it is on \((0,1)\).

Method. Use the binomial mass recurrence, starting with \(q=1-0.20=0.80\). The initial mass is

$$ a_0=q^n=0.80^5=0.32768. $$

The recurrence multiplier from mass \(a_k\) to \(a_{k+1}\) is \(\frac{5-k}{k+1}\frac{0.20}{0.80}\). Thus:

\begin{aligned} a_1&=0.32768\left(\frac{5}{1}\right)(0.25)=0.40960,\\ a_2&=0.40960\left(\frac{4}{2}\right)(0.25)=0.20480,\\ a_3&=0.20480\left(\frac{3}{3}\right)(0.25)=0.05120,\\ a_4&=0.05120\left(\frac{2}{4}\right)(0.25)=0.00640,\\ a_5&=0.00640\left(\frac{1}{5}\right)(0.25)=0.00032. \end{aligned}

The cumulative probabilities are \(C_0=0.32768\), \(C_1=0.32768+0.40960=0.73728\), \(C_2=0.94208\), \(C_3=0.99328\), \(C_4=0.99968\), and \(C_5=1.00000\). Since \(C_0<0.70\) and \(C_1\geq0.70\), the simulated count is \(1\). As a check, the six masses sum to \(0.32768+0.40960+0.20480+0.05120+0.00640+0.00032=1.00000\).

Common mistake avoided. Returning zero because \(U\) is larger than \(C_0\) ignores that inversion continues through the cumulative distribution. The first threshold that reaches \(0.70\) is \(C_1\), so the correct outcome is one detection.

Worked Example: A symmetric binomial count

A hypothetical inspection samples \(n=4\) parts. Each part independently fails inspection with probability \(p=0.50\). For \(U=0.60\), find the inverted binomial outcome.

Method. Here \(q=0.50\), and \(a_0=0.50^4=0.0625\). Since \(p/q=1\), each recurrence multiplier is just \((4-k)/(k+1)\). Therefore

\begin{aligned} a_1&=0.0625\left(\frac{4}{1}\right)=0.25,\\ a_2&=0.25\left(\frac{3}{2}\right)=0.375,\\ a_3&=0.375\left(\frac{2}{3}\right)=0.25,\\ a_4&=0.25\left(\frac{1}{4}\right)=0.0625. \end{aligned}

Accumulate in order: \(C_0=0.0625\), \(C_1=0.3125\), and \(C_2=0.6875\). Because \(C_1<0.60\leq C_2\), the returned outcome is \(2\). The mass check is \(0.0625+0.25+0.375+0.25+0.0625=1\), as required.

Common mistake avoided. The count is not the index of the first mass that exceeds \(U\). The masses themselves are \(a_k\), while inversion compares \(U\) with cumulative probabilities \(C_k\). Here \(a_2=0.375\) is not the selection threshold; \(C_2=0.6875\) is.

Worked Example: Stopping before the upper tail

A hypothetical monitoring system records whether each of six independent events triggers an alert. The probability of an alert per event is \(p=0.40\). Use \(U=0.80\) to generate the alert count.

Method. Set \(q=0.60\), so \(a_0=0.60^6=0.046656\). The recurrence ratio includes \(p/q=2/3\). Calculate successive masses and cumulative values until the threshold is reached:

\begin{aligned} a_1&=0.046656\left(6\right)\left(\frac{2}{3}\right)=0.186624, & C_1&=0.046656+0.186624=0.233280,\\ a_2&=0.186624\left(\frac{5}{2}\right)\left(\frac{2}{3}\right)=0.311040, & C_2&=0.233280+0.311040=0.544320,\\ a_3&=0.311040\left(\frac{4}{3}\right)\left(\frac{2}{3}\right)=0.276480, & C_3&=0.544320+0.276480=0.820800. \end{aligned}

The earlier cumulative values are \(C_0=0.046656\), \(C_1=0.233280\), and \(C_2=0.544320\), all below \(0.80\). Since \(C_3=0.820800\geq0.80\), return \(3\). There is no need to calculate \(a_4,a_5,\) or \(a_6\) for this draw. As a check on the stopping decision, \(C_2<0.80\leq C_3\), so three is the smallest eligible count.

Common mistake avoided. Continuing to calculate all masses is unnecessary for this inversion. The method stops as soon as the cumulative total reaches the specified uniform variate; later masses cannot change the selected count.

Conditions, checks, and implementation traps

The recurrence assumes \(0<p<1\), so that division by \(q\) is valid. At the endpoints the binomial distribution is degenerate: \(P(X=0)=1\) when \(p=0\), and \(P(X=n)=1\) when \(p=1\). These cases should be handled directly rather than passed through the recurrence. The trial count \(n\) must be a nonnegative integer and \(p\) must be a probability in \([0,1]\).

In exact arithmetic, the masses generated by the recurrence are nonnegative and sum to one. A practical check is to compute all masses for a small example and confirm that their total is one. For a single draw, however, cumulative inversion only requires calculating through the first \(k\) for which \(C_k\geq U\). The outcome can also be checked by verifying both inequalities: the previous cumulative probability is below \(U\), and the selected cumulative probability is at least \(U\). For \(k=0\), only the second condition is needed.

A frequent implementation error is to update the mass with \(\frac{n-k}{k}\) instead of \(\frac{n-k}{k+1}\). For example, with \(n=4\) and \(p=q=0.5\), the correct step from \(a_0=0.0625\) to \(a_1\) multiplies by \(4/1\), giving \(a_1=0.25\). Using \(4/0\) is undefined; using the wrong index elsewhere also breaks the probability sequence. The denominator is \(k+1\) because the recurrence moves from the mass at \(k\) to the mass at \(k+1\).

Another trap is to compare the uniform variate with \(a_k\) rather than \(C_k\), or to return the first \(k\) for which \(C_k>U\) without considering the stated inverse convention. The standard rule here is the smallest \(k\) satisfying \(C_k\geq U\). A computer implementation should also avoid rounding intermediate masses or cumulative totals too aggressively: premature rounding can move a threshold comparison across \(U\). Keep adequate numerical precision during the recursion and round only for display.

Key takeaway: For binomial inversion, initialize with \(a_0=(1-p)^n\), recursively calculate each next mass, and add it to the cumulative probability. Return the first \(k\) with \(C_k\geq U\); handle \(p=0\) and \(p=1\) separately.

This method connects the inverse transform principle to a familiar discrete count model. In simulation work, correct distributional logic matters as much as speed: a single uniform variate can produce a binomial count, but only if the masses, cumulative totals, and boundary comparison are handled consistently. These steps are aligned with the SOA Exam ASTAM (Advanced Short-Term Actuarial Mathematics) syllabus.

Check Your Understanding

For each problem, use the recursive masses and the rule that the selected value is the smallest \(k\) with \(C_k\geq U\). Probabilities are shown exactly or to enough precision to identify the threshold.

A hypothetical diagnostic device is used on \(n=3\) independent samples, with probability \(p=0.20\) of a positive result per sample. Given \(U=0.50\), what count is returned? (Answer: 0)
A hypothetical quality check has \(n=4\) independent opportunities for a defect, each with probability \(p=0.30\). Given \(U=0.65\), what count is returned? (Answer: 1)
A hypothetical alert process has \(n=5\) independent events, each triggering an alert with probability \(p=0.40\). Given \(U=0.75\), what count is returned? (Answer: 3)
A hypothetical fleet has \(n=6\) independent units, each with failure probability \(p=0.50\). Given \(U=0.90\), use the cumulative binomial masses to find the simulated failure count. (Answer: 5)
A hypothetical sensor registers an event independently on each of \(n=3\) attempts, with probability \(p=0.75\) per attempt. Given \(U=0.58\), what count is returned? (Answer: 3)

What You'll Learn