20  Distributions

A handful of distributions cover most of what you will meet. They are not arbitrary — each one is the answer to a question about a mechanism, and knowing which question tells you which distribution to reach for.

This chapter is a field guide. The connections between the families matter as much as the families themselves, because most of them are limits or special cases of one another.

20.1 A map of the common distributions

Distribution Answers Support
Bernoulli did it happen? \(\{0,1\}\)
Binomial how many out of \(n\)? \(\{0,\dots,n\}\)
Poisson how many in an interval? \(\{0,1,2,\dots\}\)
Uniform all values equally likely \([a,b]\)
Normal sums of many small effects \(\mathbb{R}\)
Exponential how long until the next event? \([0,\infty)\)
Gamma how long until the \(k\)-th? \([0,\infty)\)
Beta a proportion, or a belief about one \([0,1]\)
Student’s \(t\) a mean with unknown variance \(\mathbb{R}\)
Chi-squared a sum of squared normals \([0,\infty)\)

They are heavily interrelated: the binomial is a sum of Bernoullis, the Poisson is a limit of binomials, the gamma is a sum of exponentials, chi-squared is a special gamma, and the normal is the limit of almost everything (Section 22.4).

20.2 Bernoulli and binomial

A Bernoulli trial is a single yes/no event with success probability \(p\):

\[ \mathbb{E}[X] = p, \qquad \operatorname{Var}(X) = p(1-p) \]

The variance is largest at \(p = 0.5\) and vanishes at \(p = 0\) or \(1\) — a certain outcome has no spread.

Count the successes in \(n\) independent trials and you get the binomial:

\[ P(X = k) = \binom{n}{k}p^k(1-p)^{n-k}, \qquad \mathbb{E}[X] = np, \quad \operatorname{Var}(X) = np(1-p) \tag{20.1}\]

The \(\binom{n}{k}\) counts the orderings (Section 18.4); the rest is the probability of any one of them.

c(
  mean = 10 * 0.3,
  variance = 10 * 0.3 * 0.7,
  P_exactly_3 = dbinom(3, size = 10, prob = 0.3),
  P_at_most_3 = pbinom(3, size = 10, prob = 0.3)
)
       mean    variance P_exactly_3 P_at_most_3 
  3.0000000   2.1000000   0.2668279   0.6496107
NoteIn machine learning

Logistic regression models a Bernoulli outcome, with \(p\) depending on the features through the logistic function (Section 3.8). Its loss — binary cross-entropy — is exactly the negative log-likelihood of that Bernoulli (Section 21.5), which is why the two are always described together.

20.3 Poisson

Counts of events in a fixed interval, when events occur independently at a constant average rate \(\lambda\):

\[ P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \qquad \mathbb{E}[X] = \operatorname{Var}(X) = \lambda \tag{20.2}\]

Mean equals variance — a strong and checkable claim. Real count data is often overdispersed, with variance exceeding the mean, which is the signal to reach for a negative binomial instead.

The Poisson is the limit of the binomial when \(n\) is large, \(p\) small, and \(np = \lambda\) stays fixed — many chances, each unlikely.

k <- 0:8
rbind(
  binomial = dbinom(k, size = 1000, prob = 0.003),
  poisson = dpois(k, lambda = 3)
)
               [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
binomial 0.04956308 0.1491367 0.2241537 0.2243786 0.1682839 0.1008691
poisson  0.04978707 0.1493612 0.2240418 0.2240418 0.1680314 0.1008188
               [,7]       [,8]        [,9]
binomial 0.05033337 0.02150653 0.008032594
poisson  0.05040941 0.02160403 0.008101512
draw_line(
  x = k,
  y = list(
    `Binomial(1000, 0.003)` = dbinom(k, 1000, 0.003),
    `Poisson(3)` = dpois(k, 3)
  ),
  xlab = "k",
  ylab = "probability"
)
Figure 20.1: Binomial(1000, 0.003) against Poisson(3). The two are nearly indistinguishable — which is why the Poisson is used for rare events in large populations, where \(n\) and \(p\) may not even be separately known.

The agreement is to three decimal places, and note what the Poisson needs: only \(\lambda\). You do not have to know \(n\) or \(p\) separately, which is usually the situation with arrivals, defects, or rare disease counts.

20.4 Uniform

Every value in \([a,b]\) equally likely:

\[ f(x) = \frac{1}{b-a}, \qquad \mathbb{E}[X] = \frac{a+b}{2}, \quad \operatorname{Var}(X) = \frac{(b-a)^2}{12} \]

c(mean = (2 + 8) / 2, variance = (8 - 2)^2 / 12)
    mean variance 
       5        3

Its main role is as raw material. Every other random number in R starts as a uniform and is transformed — feeding uniforms through an inverse CDF produces any distribution you like, which is the inverse transform method.

set.seed(1)
u <- runif(5)
rbind(
  uniform = round(u, 4),
  as_exponential = round(-log(1 - u) / 0.5, 4)
)
                 [,1]   [,2]   [,3]   [,4]   [,5]
uniform        0.2655 0.3721 0.5729 0.9082 0.2017
as_exponential 0.6172 0.9308 1.7013 4.7765 0.4505

20.5 Normal

The central distribution of statistics:

\[ f(x) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \tag{20.3}\]

with mean \(\mu\) and variance \(\sigma^2\). The \(\sqrt{2\pi}\) is the Gaussian integral from Equation 14.4, present purely to make the total probability 1.

Its dominance is not aesthetic — it is the central limit theorem: sums of many independent contributions tend to normality whatever their individual shapes (Section 22.4). Measurement error, biological variation and aggregate effects are all sums of many small causes.

c(
  within_1_sd = pnorm(1) - pnorm(-1),
  within_2_sd = pnorm(2) - pnorm(-2),
  within_3_sd = pnorm(3) - pnorm(-3)
)
within_1_sd within_2_sd within_3_sd 
  0.6826895   0.9544997   0.9973002

The 68–95–99.7 rule. Note it is 95.45% within two standard deviations; the familiar 95% interval uses \(\pm1.96\).

WarningWatch out

The normal has very light tails — probability falls off like \(e^{-x^2}\). A five-sigma event has probability about \(6\times10^{-7}\), so under a normal model it is essentially impossible.

Financial returns, network traffic, and city sizes all have far heavier tails than that, and modeling them as normal will badly understate the chance of extremes. Assuming normality is a substantive assumption, not a formality.

20.6 Exponential and gamma

The exponential models waiting time until the next event, when events arrive at constant rate \(\lambda\):

\[ f(x) = \lambda e^{-\lambda x}, \qquad \mathbb{E}[X] = \frac{1}{\lambda}, \quad \operatorname{Var}(X) = \frac{1}{\lambda^2} \]

It is the continuous partner of the Poisson: if counts in an interval are Poisson, the gaps between events are exponential.

Its defining property is memorylessness:

\[ P(X > s + t \mid X > s) = P(X > t) \]

Having waited already tells you nothing about how much longer you must wait.

rate <- 0.5
c(
  P_gt_3 = 1 - pexp(3, rate),
  P_gt_5_given_gt_2 =
    (1 - pexp(5, rate)) / (1 - pexp(2, rate))
)
           P_gt_3 P_gt_5_given_gt_2 
        0.2231302         0.2231302

Identical. This is realistic for radioactive decay and reasonable for some arrival processes; it is clearly wrong for anything that ages, like machine parts or patients, which is why survival analysis needs richer families.

Sum \(k\) independent exponentials and you get the gamma, the waiting time until the \(k\)-th event:

\[ \mathbb{E}[X] = \frac{k}{\lambda}, \qquad \operatorname{Var}(X) = \frac{k}{\lambda^2} \]

c(
  gamma_mean = 3 / 0.5,
  gamma_var = 3 / 0.5^2,
  check_mean = mean(rgamma(1e5, shape = 3, rate = 0.5))
)
gamma_mean  gamma_var check_mean 
    6.0000    12.0000     5.9944

20.7 Beta

Defined on \([0,1]\), so it models proportions — or beliefs about a probability:

\[ \mathbb{E}[X] = \frac{\alpha}{\alpha+\beta} \]

The two shape parameters behave like counts of prior successes and failures, which makes Beta remarkably flexible: \(\text{Beta}(1,1)\) is uniform, large equal parameters give a tight bump at \(0.5\), and values below 1 push mass to the ends.

c(
  mean = 2 / (2 + 5),
  from_r = 2 / 7,
  sd = sqrt(2 * 5 / ((7^2) * 8))
)
     mean    from_r        sd 
0.2857143 0.2857143 0.1597191
NoteIn machine learning

Beta is the conjugate prior for the binomial: start with \(\text{Beta}(\alpha,\beta)\), observe \(s\) successes and \(f\) failures, and the posterior is \(\text{Beta}(\alpha+s, \beta+f)\) — updating is just addition. That closed form is why it appears throughout Bayesian A/B testing and Thompson sampling (Section 21.8).

20.8 Student’s t

Shaped like the normal but with heavier tails, controlled by the degrees of freedom:

\[ \operatorname{Var}(T) = \frac{\nu}{\nu-2} \quad (\nu > 2) \]

As \(\nu \to \infty\) it converges to the normal. At small \(\nu\) the tails are much heavier — and at \(\nu \leq 2\) the variance is infinite.

dfs <- c(1, 3, 10, 30)
rbind(
  df = dfs,
  P_abs_gt_2 = round(2 * (1 - pt(2, dfs)), 6)
)
               [,1]     [,2]      [,3]      [,4]
df         1.000000 3.000000 10.000000 30.000000
P_abs_gt_2 0.295167 0.139326  0.073388  0.054625
round(2 * (1 - pnorm(2)), 6)
[1] 0.0455
draw_fun(
  list(
    `t(3)` = function(x) dt(x, 3),
    `normal` = dnorm
  ),
  from = -5, to = 5,
  ylab = "density"
)
Figure 20.2: Student’s \(t\) with 3 degrees of freedom against the standard normal. The peak is lower and the tails are visibly fatter — the same total probability, redistributed toward the extremes.

The \(t\) arises when estimating a mean with the variance also estimated from the data. The extra uncertainty about \(\sigma\) is what thickens the tails, and it is why small samples use \(t\) rather than normal critical values (Section 22.5).

20.9 Chi-squared and F

The chi-squared distribution with \(k\) degrees of freedom is the sum of \(k\) squared standard normals:

\[ \chi^2_k = \sum_{i=1}^{k} Z_i^2, \qquad \mathbb{E} = k, \quad \operatorname{Var} = 2k \]

set.seed(3)
sim <- rowSums(matrix(rnorm(20000 * 3), ncol = 3)^2)
c(
  simulated_mean = mean(sim), theory = 3,
  simulated_var = var(sim), theory_var = 6
)
simulated_mean         theory  simulated_var     theory_var 
      3.032056       3.000000       6.137380       6.000000

Since variances are sums of squares, chi-squared governs the sampling distribution of a sample variance — which is why it appears in variance tests and goodness-of-fit tests.

The F distribution is a ratio of two independent chi-squareds, each divided by its degrees of freedom. Comparing two variances is a ratio, which is why F is the distribution behind ANOVA and the overall significance test in regression.

20.10 The multivariate normal

The normal generalized to vectors, parameterized by a mean vector and a covariance matrix:

\[ \mathbf{X} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \]

\(\boldsymbol{\Sigma}\) must be symmetric positive semi-definite (Section 8.9) — a matrix cannot be a covariance matrix otherwise, since \(\mathbf{v}^\top\boldsymbol{\Sigma}\mathbf{v}\) is the variance of \(\mathbf{v}^\top\mathbf{X}\) and variances cannot be negative.

Its level sets are ellipses, and this is where the linear algebra pays off: the density depends on \(\mathbf{x}\) only through the quadratic form \((\mathbf{x}-\boldsymbol{\mu})^\top\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\), so the contours are exactly the ellipses of Section 8.8 — axes along the eigenvectors of \(\boldsymbol{\Sigma}\), lengths set by its eigenvalues.

Sigma <- matrix(c(1, 0.8, 0.8, 1), nrow = 2)
ev <- eigen(Sigma, symmetric = TRUE)
c(
  eigenvalues = ev$values,
  axis_ratio = sqrt(ev$values[1] / ev$values[2])
)
eigenvalues1 eigenvalues2   axis_ratio 
         1.8          0.2          3.0 
Sinv <- solve(Sigma)
mvn_q <- function(x, y) {
  Sinv[1, 1] * x^2 + 2 * Sinv[1, 2] * x * y + Sinv[2, 2] * y^2
}
draw_contour(
  mvn_q, c(-3, 3), c(-3, 3),
  xlab = "x1", ylab = "x2"
)
Figure 20.3: Level sets of the quadratic form in a bivariate normal’s exponent, with correlation 0.8 — the same ellipses as the density’s contours, since the density is a decreasing function of this form. They tilt along the direction of shared variation, elongated by the square root of the eigenvalue ratio, here 3 to 1.

Correlation 0.8 tilts the ellipses toward the diagonal. At correlation 0 they would be circles; approaching 1 they collapse toward a line — which is the collinearity of Section 6.11, seen as a distribution.

20.11 Mixtures

A mixture draws from one of several component distributions, chosen at random:

\[ f(x) = \sum_k \pi_k f_k(x), \qquad \sum_k \pi_k = 1 \]

Mixtures produce shapes no single standard family can, most obviously multiple modes.

mix <- function(x) {
  0.6 * dnorm(x, -2, 1) + 0.4 * dnorm(x, 2, 0.7)
}
draw_fun(mix, from = -6, to = 5, ylab = "density")
Figure 20.4: A mixture of two normals. No single normal has two peaks, so a bimodal histogram is strong evidence that the population contains distinct subgroups. 
c(
  mixture_mean = 0.6 * (-2) + 0.4 * 2,
  density_at_the_mean = mix(-0.4),
  density_at_left_peak = mix(-2)
)
        mixture_mean  density_at_the_mean density_at_left_peak 
         -0.40000000           0.06719118           0.23936539

Note the mean, \(-0.4\), sits in the valley between the two peaks — a region of low density. For a multimodal distribution the mean can be a value the data rarely takes, which is a good reason to plot before summarizing.

NoteIn machine learning

Gaussian mixture models are clustering: each component is a cluster, and fitting means estimating the weights, means and covariances. Because you never observe which component produced each point, the fit uses the EM algorithm rather than a closed form.

20.12 Choosing a distribution

Work from the mechanism, not the histogram.

Ask Suggests
binary outcome? Bernoulli
number of successes in \(n\) trials? binomial
count with no fixed upper limit? Poisson
count, variance far above the mean? negative binomial
waiting time, constant rate? exponential
waiting time until the \(k\)-th event? gamma
a proportion? beta
sum of many small effects? normal
heavy tails, occasional extremes? \(t\), or a heavy-tailed family
distinct subgroups? a mixture

Two habits worth forming. Check the support: a distribution on \([0,\infty)\) cannot model something that goes negative, and this rules out candidates instantly. And check the mean–variance relationship: the Poisson insists they are equal, the normal lets them vary freely, and comparing the two in your data is a fast diagnostic.

20.13 Sampling in R

R uses a consistent four-letter scheme for every distribution:

Prefix Gives Example
d density or mass dnorm(0)
p CDF, \(P(X \leq x)\) pnorm(1.96)
q quantile, the inverse CDF qnorm(0.975)
r random draws rnorm(10) 
set.seed(7)
c(
  density_at_0 = dnorm(0),
  cdf_at_1.96 = pnorm(1.96),
  quantile_97.5 = qnorm(0.975),
  one_draw = rnorm(1)
)
 density_at_0   cdf_at_1.96 quantile_97.5      one_draw 
    0.3989423     0.9750021     1.9599640     2.2872472

The q and p functions are inverses of each other, which is how critical values are found:

c(
  round_trip = pnorm(qnorm(0.975)),
  t_critical_df10 = qt(0.975, df = 10),
  normal_critical = qnorm(0.975)
)
     round_trip t_critical_df10 normal_critical 
       0.975000        2.228139        1.959964

The \(t\) critical value is larger, which is the heavier tails demanding a wider interval for the same confidence.

20.14 Summary

Distribution Mean Variance
Bernoulli(\(p\)) \(p\) \(p(1-p)\)
Binomial(\(n,p\)) \(np\) \(np(1-p)\)
Poisson(\(\lambda\)) \(\lambda\) \(\lambda\)
Uniform(\(a,b\)) \((a+b)/2\) \((b-a)^2/12\)
Normal(\(\mu,\sigma^2\)) \(\mu\) \(\sigma^2\)
Exponential(\(\lambda\)) \(1/\lambda\) \(1/\lambda^2\)
Gamma(\(k,\lambda\)) \(k/\lambda\) \(k/\lambda^2\)
Beta(\(\alpha,\beta\)) \(\alpha/(\alpha+\beta)\) —
\(\chi^2_k\) \(k\) \(2k\)
\(t_\nu\) \(0\) \(\nu/(\nu-2)\)

20.15 Exercises

1. A fair coin is flipped 20 times. Find \(P(\text{exactly }10)\) and \(P(\text{at least }15)\).

c(
  exactly_10 = dbinom(10, 20, 0.5),
  at_least_15 = 1 - pbinom(14, 20, 0.5)
)
 exactly_10 at_least_15 
 0.17619705  0.02069473

Only 17.6% for exactly ten — the single most likely count is still unlikely, because the probability is spread across 21 possible values. Note pbinom(14, ...) for “at least 15”: the CDF is \(P(X \leq x)\), so the complement needs \(14\), not \(15\).

2. A call center receives 4 calls per hour on average. What is the probability of exactly 6 in an hour, and of none in 30 minutes?

c(
  six_in_an_hour = dpois(6, lambda = 4),
  none_in_half_an_hour = dpois(0, lambda = 2)
)
      six_in_an_hour none_in_half_an_hour 
           0.1041956            0.1353353

The rate scales with the interval: half an hour has \(\lambda = 2\). The second answer is \(e^{-2}\), and it also equals \(P(\text{waiting time} > 0.5)\) under the matching exponential — the two views of the same process agree.

3. Verify that the exponential is memoryless for a different pair of times.

c(
  P_gt_4 = 1 - pexp(4, rate),
  P_gt_10_given_gt_6 =
    (1 - pexp(10, rate)) / (1 - pexp(6, rate)),
  P_gt_1 = 1 - pexp(1, rate),
  P_gt_7_given_gt_6 =
    (1 - pexp(7, rate)) / (1 - pexp(6, rate))
)
            P_gt_4 P_gt_10_given_gt_6             P_gt_1  P_gt_7_given_gt_6 
         0.1353353          0.1353353          0.6065307          0.6065307

Both pairs match. The exponential is the only continuous distribution with this property, which is a strong constraint — and the reason it is often too simple for real waiting times.

4. Compare tail probabilities for \(t\) with 5 df against the normal at 2, 3 and 4 standard deviations.

q <- c(2, 3, 4)
rbind(
  t5 = 2 * (1 - pt(q, 5)),
  normal = 2 * (1 - pnorm(q)),
  ratio = (2 * (1 - pt(q, 5))) / (2 * (1 - pnorm(q)))
)
             [,1]         [,2]         [,3]
t5     0.10193948  0.030099248 1.032342e-02
normal 0.04550026  0.002699796 6.334248e-05
ratio  2.24041511 11.148711678 1.629778e+02

The ratio grows rapidly: further out, the \(t\) is many times more likely to produce an extreme value. At 4 standard deviations it is well over a hundred times more likely.

Tail behavior diverges where it matters most — in the region you use for rare-event probabilities and significance thresholds.

5. Show that a Beta(1,1) is uniform, and find the posterior after 7 successes in 10 trials starting from that prior.

c(
  beta11_at_0.2 = dbeta(0.2, 1, 1),
  beta11_at_0.9 = dbeta(0.9, 1, 1),
  posterior_mean = (1 + 7) / (1 + 7 + 1 + 3),
  sample_proportion = 7 / 10
)
    beta11_at_0.2     beta11_at_0.9    posterior_mean sample_proportion 
        1.0000000         1.0000000         0.6666667         0.7000000

The Beta(1,1) density is flat at 1 everywhere — the uniform. After 7 successes and 3 failures the posterior is Beta(8,4), with mean \(8/12 = 0.667\).

It sits slightly below the sample proportion of 0.7, pulled toward the prior’s 0.5. That shrinkage is the same phenomenon as ridge regression (Section 6.11), arrived at from the Bayesian side.

6. For the covariance matrix used in Figure 20.3, what correlation would make the contours circular? What happens as correlation approaches 1?

ratio_at <- function(r) {
  S <- matrix(c(1, r, r, 1), nrow = 2)
  e <- eigen(S, symmetric = TRUE)$values
  sqrt(e[1] / e[2])
}
rhos <- c(0, 0.5, 0.8, 0.95, 0.99)
rbind(
  correlation = rhos,
  axis_ratio = round(sapply(rhos, ratio_at), 3)
)
            [,1]  [,2] [,3]  [,4]   [,5]
correlation    0 0.500  0.8 0.950  0.990
axis_ratio     1 1.732  3.0 6.245 14.107

At \(\rho = 0\) the ratio is 1 — circular contours, since the eigenvalues are equal. As \(\rho \to 1\) the ratio grows without bound and the ellipse collapses toward a line: the two variables become redundant, \(\boldsymbol{\Sigma}\) approaches singularity, and the distribution loses a dimension.

That is collinearity described distributionally, and it is why a near-singular covariance matrix is as much a modeling problem as a numerical one.