Negative Binomial Distribution
Probability Distributions
Waiting for the r-th Success — Overdispersed Counts
The negative binomial distribution generalizes the geometric: instead of waiting for the first success, we wait for the -th success. It is the go-to model for overdispersed count data.
- Insurance — claims per policyholder (variance > mean)
- Epidemiology — disease cases per region (heterogeneous rates)
- Ecology — species counts per quadrat (aggregated populations)
- Transportation — passengers per bus (bursty arrivals)
When the Poisson's mean-equals-variance assumption fails, the negative binomial saves the day.
Core Concepts
The negative binomial distribution generalizes the geometric distribution: instead of waiting for the first success, we wait for the -th success. It arises naturally as a sum of independent geometric random variables and serves as a flexible model for overdispersed count data.
PMF Derivation
Mean and Variance
Cumulant Generating Function
The cumulant generating function is:
yielding cumulants:
The Negative Binomial as a Poisson-Gamma Mixture
Overdispersion in Practice
Worked Example: Call Center Modeling
Python Implementation
import numpy as np
from scipy import stats
np.random.seed(42)
# Negative binomial parameters
r, p = 6, 0.4
n = 10000
# Simulate
samples = np.random.negative_binomial(r, p, size=n)
# Verify mean and variance
mean_theory = r * (1 - p) / p
var_theory = r * (1 - p) / p**2
print(f"NB(r={r}, p={p}):")
print(f" Theoretical mean: {mean_theory:.4f}, variance: {var_theory:.4f}")
print(f" Empirical mean: {np.mean(samples):.4f}, variance: {np.var(samples, ddof=0):.4f}")
# Show relationship to geometric sum
geom_samples = np.random.geometric(p, size=(n, r))
geom_sum = geom_samples.sum(axis=1) - r # convert from "trials" to "failures"
print(f"\n Sum of {r} Geometric(p={p}): mean={np.mean(geom_sum):.4f}, var={np.var(geom_sum, ddof=0):.4f}")
print(f" (Should match NB values above)")
# Show Poisson-Gamma mixture
lam = np.random.gamma(shape=r, scale=(1-p)/p, size=n)
poisson_samples = np.random.poisson(lam)
print(f"\n Poisson-Gamma mixture: mean={np.mean(poisson_samples):.4f}, var={np.var(poisson_samples, ddof=0):.4f}")
Python Implementation: Overdispersion Detection
import numpy as np
from scipy import stats
np.random.seed(42)
# Generate overdispersed count data (NB instead of Poisson)
true_mu = 5
true_alpha = 2.0 # dispersion parameter: r = 1/alpha
r = 1 / true_alpha
p = r / (r + true_mu)
n = 500
data = np.random.negative_binomial(r, p, size=n)
sample_mean = np.mean(data)
sample_var = np.var(data, ddof=1)
dispersion_ratio = sample_var / sample_mean
print(f"Overdispersion Test")
print(f" Sample mean: {sample_mean:.4f}")
print(f" Sample var: {sample_var:.4f}")
print(f" Var/Mean: {dispersion_ratio:.4f}")
print(f" (Var/Mean > 1 suggests overdispersion; Poisson requires Var/Mean ≈ 1)")
# Method of moments estimates for NB
p_hat = sample_mean / sample_var
r_hat = sample_mean * p_hat / (1 - p_hat)
print(f"\n Method of moments estimates:")
print(f" p̂ = {p_hat:.4f}, r̂ = {r_hat:.4f}")
print(f" Estimated mean: {r_hat*(1-p_hat)/p_hat:.4f}")
print(f" Estimated var: {r_hat*(1-p_hat)/p_hat**2:.4f}")