Generates pseudorandom numbers from a negative binomial distribution.

Routine RNNBN generates pseudorandom numbers from a negative binomial distribution with parameters RK and P. RK and P must be positive and P must be less than 1. The probability function (with r = RK and p = P) is

for x = 0, 1, 2, ….

If r is an integer, the distribution is often called the Pascal distribution and can be thought of as modeling the length of a sequence of Bernoulli trials until r successes are obtained, where p is the probability of getting a success on any trial. In this form, the random variable takes values r, r + 1, r + 2, … and can be obtained from the negative binomial random variable defined above by adding r to the negative binomial variable. This latter form is also equivalent to the sum of r geometric random variables defined as taking values 1, 2, 3, … .

If rp/(1 ‑ p) is less than 100 and (1 ‑ p)r is greater than the machine epsilon, RNNBN uses the inverse CDF technique; otherwise, for each negative binomial deviate, RNNBN generates a gamma (r, p/(1 ‑ p)) deviate Y and then generates a Poisson deviate with parameter Y.

In this example, RNNBN is used to generate five pseudorandom deviates from a negative binomial (Pascal) distribution with parameter r equal to 4 and p equal to 0.3.