R Notebook

R Functions dnbinom, pnbinom, and rnbinom

Random variable X is distributed X∼NB(r,p) with mean \[μ=\frac{r}{p}\] and variance \[σ2=r(1−p)/p2\] if X is the count of independent Bernoulli trials required to achieve the rth successful trial when the probability of success is constant p.

The probability of X=n trials is

\[ f(X=n)= {n−1 \choose r−1} p^r (1−p)^{n−r} \]

1. R function dnbinom(x, size, prob) is the probability of x failures prior to the rth success (note the difference) when the probability of success is prob.

2. R function pgeom(q, prob, lower.tail) is the cumulative probability (lower.tail = TRUE for left tail, lower.tail = FALSE for right tail) of less than or equal to q failures prior to success.

Example

An oil company has a p = 0.20 chance of striking oil when drilling a well. What is the probability the company drills x = 7 wells to strike oil r = 3 times?

r = 3
p = 0.20
n = 7 - r
# exact
dnbinom(x = n, size = r, prob = p)

[1] 0.049152

# simulated
set.seed(1)
mean(rnbinom(n = 10000000, size = r, prob = p) == n)

[1] 0.0491621

library(dplyr)
library(ggplot2)

data.frame(x = 0:25, prob = dnbinom(x = 0:25, size = r, prob = p)) %>%
  mutate(Failures = ifelse(x == n, n, "other")) %>%
ggplot(aes(x = factor(x), y = prob, fill = Failures)) +
  geom_col() +
  geom_text(
    aes(label = round(prob,2), y = prob + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0
  ) +
  labs(title = "Probability of r = 3 Successes in X = 7 Trials",
       subtitle = "NB(3,.2)",
       x = "Failed Trials (X - r)",
       y = "Probability") 

#Example

What is the expected number of trials to achieve r = 3 successes when the probability of success is p = 0.2 ?

r = 3
p = 0.20
# mean
# exact
r / p

[1] 15

simulated

mean(rnbinom(n = 1000000, size = 3, prob = p)) + r

[1] 15.00988

# Variance
# exact
r * (1 - p) / p^2

[1] 60

simulated

var(rnbinom(n = 10000000, size = r, prob = p))

[1] 60.03182

library(dplyr)
library(ggplot2)

data.frame(x = 1:20, 
           pmf = dnbinom(x = 1:20, size = r, prob = p),
           cdf = pnbinom(q = 1:20, size = r, prob = p, lower.tail = TRUE)) %>%
ggplot(aes(x = factor(x), y = cdf , fill = cdf)) +
  geom_col() +
  geom_text(
    aes(label = round(cdf,2), y = cdf + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0
  ) +
  labs(title = "Cumulative Probability of X = x failed trials to achieve 3rd success",
       subtitle = "NB(3,.2)",
       x = "Failed Trials (x)",
       y = "probability") 

#Question

You are surveying people exiting from a polling booth and asking them if they voted independent. The probability (p) that a person voted independent is 20%. What is the probability that 70 people must be asked before you can find 5 people who voted independent?

Code: We use dnbinom(x, size, prob) - x= number of failures that need to happen before reaching the required number of successes.

• size= number of successes.

• prob= the probability of success in each trial. Infinite and missing values are not allowed. To answer the above question:

p =.2
r = 5 
n = 70 -5

dnbinom(x= fail,size= r, prob = p)

[1] 6.058753e-05

To plot the probability mass function for a negative binomial function in R, we can use the following functions:

x <- 65
size <- 5
prob = .2
dnbinom(x, size, prob) # is used to create the probability mass function

[1] 0.00013892

plot(0:x, dnbinom(0:x, size, prob) , type = 'h') # to plot the probability mass function, specifying the plot to be a histogram (type=’h’)

NA
NA

Referance

https://rpubs.com/mpfoley73/458738