Binomial
A single Bernoulli trial gives a \(0\) or a \(1\); counting how many successes occur across \(n\) such trials is the next-most-common discrete distribution.
A Binomial random variable counts the number of successes in \(n\) independent Bernoulli trials, each with success probability \(\theta\). Its PMF is \(\text{Prob}(X_{i} = x) = \binom{n}{x} \theta^{x} (1-\theta)^{n-x}\) for \(x~\text{in}~\{0, 1, \ldots, n\}\), mean \(n\theta\), and variance \(n\theta(1-\theta)\).
The Binomial distribution is useful for any situation that aggregates independent yes/no events: defective items in a batch, customers who enter a store, voters who pick a candidate, free throws that go in. The full specification gathers sample space, PMF, and CDF in one place: \[\begin{aligned}
X_{i} &~\text{in}~\{0,1,\ldots,n\} \\
\text{Prob}(X_{i}=x) &= \binom{n}{x}\theta^{x}(1-\theta)^{n-x} \\
F(x) = \text{Prob}(X_{i} \leq x) &= \sum_{k=0}^{\lfloor x \rfloor} \text{Prob}(X_{i}=x)
\end{aligned}\]
The simplest case: how many heads will I get when I flip a fair coin twice? We can first see that there are four potential outcomes: \(\{ (1,1), (1,0), (0,1), (0,0) \}\), each of which are equally likely. The corresponding sums of \(\{2,1,1,0\}\) are equally likely, which leads to probabilities \(\text{Prob}(X_{i}=0)=1/4\), \(\text{Prob}(X_{i}=1)=2/4\), and \(\text{Prob}(X_{i}=2)=1/4\).
If the coin is unfair, with probability of heads equal to \(0.25\), then the four outcomes \(\{(1,1), (1,0), (0,1), (0,0)\}\) are no longer equally likely: they have probabilities \(\{0.0625, 0.1875, 0.1875, 0.5625\}\). Collecting the matching sums, \(\text{Prob}(X_{i}=2)=0.0625\), \(\text{Prob}(X_{i}=1)=0.1875+0.1875=0.375\), and \(\text{Prob}(X_{i}=0)=0.5625\).
When I flip a coin three times, the potential outcomes are \[\begin{gathered}
\{
(1,1,1), (1,1,0), (1,0,1), (1,0,0),\\
(0,1,1), (0,1,0), (0,0,1), (0,0,0)
\}
\end{gathered}\] Assuming the coin is fair, compute the probabilities for each possible sum. Assuming the coin is unfair, with probability of heads equal to \(0.25\), compute the probabilities for each possible sum.
In general, we use the Binomial distribution. The \(\binom{n}{x}\) term is called the binomial coefficient and it counts all the different groups that lead to the same sum. The other terms intuitively reflect the probabilities of getting \(x\) successes and \(n-x\) failures.
For example, suppose \(n=15\) people walk past a store each hour and each has a \(\theta=30\%\) chance of entering. The number of customers who enter is \(X_{i} \sim \text{Binom}(15, 0.3)\).
Code
## Store: 15 passers-by, 30% chance each enters
n <- 15
theta <- 0.3
# PMF and CDF side by side
par(mfrow=c(1, 2))
# PMF plot
x <- seq(0, n)
f_x <- dbinom(x, n, theta)
plot(x, f_x, type = 'h', col = rgb(0, 0, 1, .8),
main=NA, xlab = 'Customers entering', ylab = 'Prob(X = x)')
title(bquote(paste('Binom(', .(n), ', ', .(theta), ')')), font.main=1)
points(x, f_x, pch = 16, col = rgb(0, 0, 1, .8))
## Theoretical CDF
F_x <- pbinom(x, n, theta)
plot(x, F_x, type='s', col=rgb(0, 0, 1, .8),
xlab = 'Customers entering', ylab = 'F(x)', main=NA)
points(x, F_x, pch=16, col=rgb(0, 0, 1, .8))
title(bquote(paste('Binom(', .(n), ', ', .(theta), ') CDF')), font.main=1)
A quality inspector checks \(n=20\) items from a production line with a \(5\%\) defect rate (\(\theta=0.05\)). What is the probability of finding \(2\) or fewer defectives?
We need \(\text{Prob}(X_{i} \leq 2) = F(2)\) where \(X_{i} \sim \text{Binom}(20, 0.05)\). By hand, the first few terms are \[\begin{aligned}
\text{Prob}(X_{i}=0) &= \binom{20}{0}(0.05)^{0}(0.95)^{20} \approx 0.36 \\
\text{Prob}(X_{i}=1) &= \binom{20}{1}(0.05)^{1}(0.95)^{19} \approx 0.38 \\
\text{Prob}(X_{i}=2) &= \binom{20}{2}(0.05)^{2}(0.95)^{18} \approx 0.19
\end{aligned}\] So \(\text{Prob}(X_{i} \leq 2) \approx 0.36 + 0.38 + 0.19 = 0.93\). There is roughly a \(93\%\) chance the inspector finds \(2\) or fewer defectives.
Code
# Individual terms
x <- c(0, 1, 2)
px <- dbinom(x, size=20, prob=0.05)
round(px, 2)
## [1] 0.36 0.38 0.19
# CDF: Prob(X <= 2)
pbinom(2, size=20, prob=0.05)
## [1] 0.9245163
Suppose that employees at a company are \(30\%\) female and \(70\%\) male. If we select a random sample of eight employees, what is the probability that \(2\) in the sample are female? That \(\geq 2\) are female?
Try conducting a simulation first,
Code
n <- 8
theta <- 0.3
x_sim <- rbinom(1000, n, theta)
For \(n\) coin flips, show that \(\mathbb{E}[X_{i}]=n\theta\) and \(\mathbb{V}[X_{i}]=n\theta(1-\theta)\). Build on the case \(n=2\) above, then consider case \(n=3\), and then the general case of any \(n\).
Code
# Compare empirical simulation vs theoretical
n <- 2
theta <- 3/4
x_sim <- rbinom(1000, size=n, prob=theta)
# Mean
c(empirical = mean(x_sim), theoretical = n*theta)
## empirical theoretical
## 1.511 1.500
# Variance
c(empirical = var(x_sim), theoretical = n*theta*(1-theta))
## empirical theoretical
## 0.3782573 0.3750000
The shape of the Binomial changes with both \(n\) and \(\theta\). For the store, compare a quiet hour with \(n=15\) passers-by to a busy hour with \(n=100\), and a \(10\%\) chance of entering to a \(50\%\) chance. With \(\theta=0.5\) the PMF is symmetric around \(n/2\), while a small \(\theta\) piles the probability near \(0\) and gives a long right tail. As \(n\) grows, the PMF spreads out and becomes more bell-shaped, even when \(\theta\) is small.
Code
# Shape varies with both n and theta
# Columns vary the success probability theta, rows vary the number of trials n
par(mfrow=c(2, 2))
for(n in c(15, 100)) {
for(theta in c(0.1, 0.5)) {
x <- seq(0, n)
f_x <- dbinom(x, n, theta)
plot(x, f_x, type='h', col=rgb(0, 0, 1, .8),
xlab='Customers entering', ylab='Prob(X = x)', main=NA)
points(x, f_x, pch=16, col=rgb(0, 0, 1, .8))
title(bquote(paste('Binom(', .(n), ', ', .(theta), ')')), font.main=1)
}
}
Poisson
When the trials are too many and individually too rare to count, the Binomial limits to a simpler model parameterized by a single rate.
A Poisson random variable counts events in a fixed interval when events are rare, independent, and arrive at a constant rate \(\lambda\). Its PMF is \(\text{Prob}(X_{i} = x) = e^{-\lambda} \lambda^{x} / x!\) for \(x~\text{in}~\{0, 1, 2, \ldots\}\), and the mean and variance are both \(\lambda\).
The Poisson distribution is useful for count outcomes where events are sparse but the opportunities are many: arrivals at a store, decay events from a radioactive source, typos per page. The single parameter \(\lambda\) is the average number of events per interval, and the full specification adds the CDF: \[\begin{aligned}
X_{i} &~\text{in}~\{0,1,\ldots\} \\
\text{Prob}(X_{i}=x) &= e^{-\lambda}\lambda^x/x! \\
F(x) = \text{Prob}(X_{i} \leq x) &= \sum_{k=0}^{\lfloor x \rfloor} \text{Prob}(X_{i}=x)
\end{aligned}\] Note the infinite discrete support \(\{0,1,2,\ldots\}\): there is no upper bound on the count, though large values become vanishingly unlikely.
For example, a small store averages \(\lambda = 3.5\) customers per hour. In any given hour the actual count might be 2, or 5, or 7. For our store, \(\lambda = 3.5\) means we expect 3.5 customers per hour on average and the following assumptions hold
- Arrivals are independent of each other (one customer entering does not make the next more or less likely)
- The average rate \(\lambda\) is constant over the interval
- Two customers cannot arrive at exactly the same instant
We can then compute a few probabilities by hand. For example, the chance that zero customers arrive in a given hour is \[\begin{aligned}
\text{Prob}(X_{i}=0) &= e^{-3.5} \cdot 3.5^{0} / 0! = e^{-3.5} \approx 0.030
\end{aligned}\] The chance of exactly 2 customers is \[\begin{aligned}
\text{Prob}(X_{i}=2) &= e^{-3.5} \cdot 3.5^{2} / 2! = e^{-3.5} \cdot 6.125 \approx 0.185
\end{aligned}\]
Code
## Store: lambda = 3.5 customers per hour
lambda <- 3.5
# PMF and CDF side by side
par(mfrow=c(1, 2))
# PMF plot
x <- seq(0, 15)
f_x <- dpois(x, lambda)
plot(x, f_x, type='h', col=rgb(0, 0, 1, .8),
xlab = 'Customers per hour', ylab = 'Prob(X = x)', main=NA)
points(x, f_x, pch = 16, col = rgb(0, 0, 1, .8))
title(bquote(paste('Pois(', .(lambda), ')')), font.main=1)
# CDF plot
F_x <- ppois(x, lambda)
plot(x, F_x, type='s', col=rgb(0, 0, 1, .8),
xlab = 'Customers per hour', ylab = 'F(x)', main=NA)
points(x, F_x, pch=16, col=rgb(0, 0, 1, .8))
title(bquote(paste('Pois(', .(lambda), ') CDF')), font.main=1)
Show that \(\mathbb{E}[X_{i}] = \mathbb{V}[X_{i}]= \lambda\). That is, the Poisson mean and variance are both equal to the rate parameter.
Code
# Compare empirical simulation vs theoretical
lambda <- 3.5
x_sim <- rpois(1e4, lambda)
c(emp_mean = mean(x_sim), th_mean = lambda)
## emp_mean th_mean
## 3.5075 3.5000
c(emp_var = var(x_sim), th_var = lambda)
## emp_var th_var
## 3.513495 3.500000
During the lunch rush, the store averages \(\lambda=4.5\) customers per hour. What is the probability that more than \(7\) customers arrive in a given hour?
We need \(\text{Prob}(X_{i} > 7) = 1 - \text{Prob}(X_{i} \leq 7) = 1 - F(7)\). \[\begin{aligned}
\text{Prob}(X_{i} > 7) &= 1 - \sum_{k=0}^{7} \frac{e^{-4.5} \cdot 4.5^{k}}{k!}
\end{aligned}\]
Code
# Prob(X > 7)
1 - ppois(7, lambda=4.5)
## [1] 0.08658647
# Verify: sum the PMF from 0 to 7
sum(dpois(0:7, lambda=4.5))
## [1] 0.9134135
There is roughly a \(9\%\) chance of more than \(7\) arrivals in a given hour.
The Poisson assumes the mean equals the variance. In real count data the variance is often larger than the mean (called overdispersion), because real-world arrivals are bunchier than the Poisson assumes. Before using a Poisson, compute the ratio \(\text{var}(x)/\text{mean}(x)\) on the data; a ratio well above \(1\) is a warning sign.
The shape of the Poisson also changes with \(\lambda\). A quiet specialty shop (\(\lambda=1\)) has a strongly right-skewed arrival pattern, while a busy supermarket (\(\lambda=15\)) looks nearly Normal.
Code
# Shape varies with lambda
par(mfrow=c(1, 3))
for(lam in c(1, 5, 15)) {
x <- seq(0, 30)
f_x <- dpois(x, lam)
plot(x, f_x, type='h', col=rgb(0, 0, 1, .8),
xlab = 'Customers per hour', ylab = 'Prob(X = x)', main=NA)
points(x, f_x, pch=16, col=rgb(0, 0, 1, .8))
title(bquote(paste('Pois(', .(lam), ')')), font.main=1)
}
