R21. Negative Binomial Distribution in R

In negative binomial distribution, we find probability of k successes in n trials, with the requirement that the last trial be a success. This requirement is due to the fact that the total probability is 1, and we not not want any double counting.

Thus, if we have 4 trials and we want 3 successes we can have these combinations (p = success, q = failure):

p p q p

p q p p

q p p p

The above can be found from choose(3,2) = 3.

Thus, if n is the number of trials, and k is the number of success, we have choose(n-1,k-1) combinations. This is due to the loss of freedom of the final trial. Each of these combinations has the probability p^k*q^(n-k), where q, failure, is 1-p.

We can rewrite the formulas in terms of failures, r, where r = n - k. This distribution goes from r = 0, to infinity. However, it goes to zero very fast. The reason we change the formulas to depend on k and r is so we can compare it to the R function for negative binomial distribution, dnbinom.

# ex21.R
N <- 12
p <- 0.55
k <- 3 
suc <- numeric(N)
for (r in 0:N-1) {
  suc[r+1] <- choose(r+k-1,k-1)*p^k*(1-p)^r
}
dis <- dnbinom(0:(N-1), size = k, prob = p)
# dis == suc
err <- sum(dis-suc)
cat("err =",err)
names(suc) <- 0:(N-1)
barplot(suc, xlab = 'r fails',
        ylab = paste(k,'suc with 0:',N-1,'fails',collapse=''),
        main = 'Negative Binomial Distribution',
        sub = paste('p =',p, collapse=''))
# err = -1.539567e-16

Output: