5. Common Probability Distributions

A small number of distributions show up again and again because each one comes from a simple mechanism: independent trials, waiting for an event, multiplying many factors, summing many small effects. This chapter covers the main discrete and continuous families, the normal distribution and the chi-square, $t$ and $F$ distributions derived from it, and a first look at stochastic processes, where these distributions describe how randomness unfolds over time. If you know the mechanism behind each distribution, you can pick a model for new data from how the data were generated instead of trying distributions until one fits.

<a id="51-discrete-distributions"></a>

5.1 Discrete Distributions

Discrete distributions model counts: how many successes, how many attempts until something happens, how many events in an interval. Each one comes with a short story, and the main skill is matching a situation to the right story.

Discrete Uniform

Why it matters. When there is no reason to prefer one outcome over another, as with a fair die, a shuffled deck or a random index into an array, every value is equally likely. The discrete uniform distribution is the model for "pick one at random". It is the starting point of classical probability (Section 3.2), the target of every random-integer function in a programming library, and a common source of subtle bugs when it is implemented incorrectly.

Definition. $X$ has a discrete uniform distribution on the integers $a, a + 1, \dots, b$, written $X \sim$ DiscreteUniform${a, \dots, b}$, if

$$p(k) = \frac{1}{n}, \qquad k = a, a + 1, \dots, b, \qquad n = b - a + 1.$$

Its CDF is a staircase with $n$ equal steps: $F(x) = (\lfloor x \rfloor - a + 1)/n$ for $a \le x \le b$.

Mean and variance. For ${1, \dots, n}$, the formulas for the sums of the first $n$ integers and squares give

$$E[X] = \frac{n + 1}{2}, \qquad E[X^2] = \frac{(n + 1)(2n + 1)}{6}, \qquad \operatorname{Var}(X) = \frac{(n + 1)(2n + 1)}{6} - \frac{(n + 1)^2}{4} = \frac{n^2 - 1}{12}.$$

Shifting by a constant does not change the variance (Section 4.3), so for ${a, \dots, b}$:

Quantity	Value
Mean	$\frac{a + b}{2}$
Variance	$\frac{n^2 - 1}{12}$, with $n = b - a + 1$
Median	$\frac{a + b}{2}$ (any value between the two middle values if $n$ is even)
Mode	every value in the support

Analogy. A lottery drum with $n$ identical balls, mixed well: each ball has the same chance of coming out. All the ways a "random" pick fails in practice come from balls that are not identical (some are heavier) or a drum that is not well mixed.

A classic implementation bug: modulo bias. A generator produces uniform integers in ${0, \dots, M - 1}$ and you need a digit in ${0, \dots, 9}$. Taking r % 10 is uniform only if 10 divides $M$. With $M = 256$ (one random byte), the values 250 to 255 form an incomplete block that maps to the digits 0 to 5, so those digits are slightly more likely. The fix is rejection: discard the values of the incomplete block and draw again.

  byte        0 .. 9 | 10 .. 19 | ... | 240 .. 249 | 250 .. 255
  byte % 10   0 .. 9 |  0 .. 9  | ... |   0 .. 9   |   0 .. 5
              '----------------------------------'  '---------'
              25 complete blocks, each digit 25 times  reject these

import numpy as np
from collections import Counter
from fractions import Fraction

# Exact mean and variance of a fair die, X ~ DiscreteUniform{1, ..., 6}
n = 6
mean = sum(Fraction(k, n) for k in range(1, n + 1))
var = sum(Fraction(k * k, n) for k in range(1, n + 1)) - mean ** 2
print(f"die: mean = {mean}, variance = {var}, (n^2 - 1)/12 = {Fraction(n * n - 1, 12)}")

# Modulo bias, exactly: each byte value 0..255 occurs once
counts = Counter(b % 10 for b in range(256))
print("byte % 10 :", " ".join(f"{d}:{counts[d]}/256" for d in range(10)))

# The same by simulation, and the rejection fix
rng = np.random.default_rng(5)
raw = rng.integers(0, 256, size=1_200_000)          # uniform random bytes
biased = raw[:1_000_000] % 10
fixed = raw[raw < 250][:1_000_000] % 10             # reject the incomplete block 250..255
print(f"P(digit <= 5): biased = {np.mean(biased <= 5):.4f} (exact 156/256 = {156 / 256:.4f}), "
      f"rejection = {np.mean(fixed <= 5):.4f} (target 0.6000)")

die: mean = 7/2, variance = 35/12, (n^2 - 1)/12 = 35/12
byte % 10 : 0:26/256 1:26/256 2:26/256 3:26/256 4:26/256 5:26/256 6:25/256 7:25/256 8:25/256 9:25/256
P(digit <= 5): biased = 0.6094 (exact 156/256 = 0.6094), rejection = 0.6000 (target 0.6000)

Library functions such as rng.integers(low, high) and Python's random.randrange already use rejection internally; the bug appears when you build your own mapping from raw random bits.

Key takeaways

• The discrete uniform gives each of $n$ consecutive integers probability $1/n$.
• Mean $(a + b)/2$, variance $(n^2 - 1)/12$ with $n = b - a + 1$.
• It models "no preference" and underlies classical probability by counting.
• Mapping random integers to a range with % introduces bias unless the range size divides the generator's range; use rejection or a library function.

🧪 Practice

1. $X$ is uniform on ${0, 1, \dots, 9}$. Compute $E[X]$, $\operatorname{Var}(X)$ and $P(X^2 > 20)$.
2. Let $X$ and $Y$ be independent fair dice. Show that $X + Y$ is not uniform, while $(X + Y) \bmod 6$ is uniform on ${0, \dots, 5}$. Explain the difference.
3. A generator returns uniform integers in ${0, \dots, 2^{31} - 1}$, and you map them to ${0, \dots, 999,999,999}$ with %. Which outputs are favored, and by what relative factor?
4. Interview: (German tank problem) Captured tanks carry serial numbers 1 to $N$, and you observe four of them: 19, 40, 42, 60. Estimate $N$. Hint: the sample maximum is always at most $N$, so it underestimates on average; think about the expected gap between consecutive observed serial numbers and add it to the maximum.

Bernoulli

Why it matters. Many random questions have a yes/no answer: did the user click, is the part defective, did the patient recover, did the request fail. The Bernoulli distribution is the model for one such trial. It is the building block of the binomial, geometric and negative binomial distributions, and every indicator random variable from Section 4.1 is Bernoulli.

Definition. $X \sim$ Bernoulli$(p)$, with $0 \le p \le 1$, if

$$P(X = 1) = p, \qquad P(X = 0) = 1 - p, \qquad \text{or compactly} \qquad p(x) = p^x (1 - p)^{1 - x}, \quad x \in {0, 1}.$$

The outcome 1 is called a success and 0 a failure, whatever the actual meaning ("success" can be "the server crashed").

Moments. Because $X \in {0, 1}$, $X^k = X$ for every $k \ge 1$, so every raw moment equals $p$:

$$E[X] = p, \qquad E[X^2] = p, \qquad \operatorname{Var}(X) = p - p^2 = p(1 - p).$$

The variance is largest at $p = 1/2$ (maximum uncertainty) and vanishes at $p = 0$ or $p = 1$ (no uncertainty at all).

$p$	0.01	0.1	0.3	0.5	0.7	0.9	0.99
$p(1 - p)$	0.0099	0.09	0.21	0.25	0.21	0.09	0.0099

Generating a Bernoulli variable. If $U \sim$ Uniform$(0, 1)$, then $X = \mathbf{1}{U < p}$ is Bernoulli$(p)$, because $P(U < p) = p$. This is inverse transform sampling (Section 4.4) in its simplest form.

The indicator method. Any count is a sum of indicators, and expectation is linear even without independence (Section 4.3). If $X = \mathbf{1}{A_1} + \dots + \mathbf{1}{A_n}$, then $E[X] = P(A_1) + \dots + P(A_n)$. This often avoids computing the distribution of $X$ at all.

import numpy as np

rng = np.random.default_rng(1)

# Bernoulli(p) from uniform numbers: X = 1{U < p}
p = 0.3
x = (rng.random(200_000) < p).astype(int)
print(f"mean = {x.mean():.4f} (p = {p}), variance = {x.var():.4f} (p(1-p) = {p * (1 - p):.4f}), "
      f"E[X^3] = {np.mean(x ** 3):.4f} (also p)")

# Indicator method: four ads shown to one user, with different click probabilities
click_p = np.array([0.20, 0.05, 0.10, 0.40])
clicks = (rng.random((200_000, 4)) < click_p).sum(axis=1)   # each row is one user
print(f"expected clicks: sum of p_i = {click_p.sum():.2f}, simulated = {clicks.mean():.4f}")
print(f"P(no click): product of (1 - p_i) = {np.prod(1 - click_p):.4f}, simulated = {np.mean(clicks == 0):.4f}")

mean = 0.3002 (p = 0.3), variance = 0.2101 (p(1-p) = 0.2100), E[X^3] = 0.3002 (also p)
expected clicks: sum of p_i = 0.75, simulated = 0.7496
P(no click): product of (1 - p_i) = 0.4104, simulated = 0.4116

The expected number of clicks needs only the sum of the probabilities. The probability of no click needs a product, and that step uses the independence of the four trials.

Key takeaways

• Bernoulli$(p)$ is a single trial with outcomes 1 (probability $p$) and 0.
• $E[X] = p$, $\operatorname{Var}(X) = p(1 - p)$, which is largest at $p = 1/2$.
• $\mathbf{1}{U < p}$ turns a uniform random number into a Bernoulli variable.
• Writing counts as sums of indicators gives expected values by adding probabilities, with or without independence.

🧪 Practice

1. Use calculus to show that $p(1 - p)$ is maximized at $p = 1/2$. What is the largest possible standard deviation of a Bernoulli variable?
2. $X \sim$ Bernoulli$(p)$. What are the distributions of $1 - X$, $X^2$ and $X(1 - X)$?
3. In a group of 30 people, let $X$ be the number of pairs who share a birthday. Use indicators to compute $E[X]$ without finding the distribution of $X$.
4. Interview: You have a coin with an unknown bias $p \in (0, 1)$. How can you use it to produce a perfectly fair random bit? Hint: toss the coin twice and compare the probabilities of HT and TH; decide what to do with HH and TT (von Neumann's trick).

Binomial

Why it matters. Count the successes in a fixed number of independent, identical yes/no trials: conversions among 1,000 visitors, defective chips in a batch of 50, correct answers when guessing on a quiz. The binomial distribution is the model for all of them. It is the basis of inference for proportions (Sections 8.3 and 9.2) and of A/B testing.

The story. Run $n$ independent Bernoulli$(p)$ trials and let $X$ be the number of successes. Four conditions must hold, remembered as BINS:

Letter	Condition	Violated when
Binary	each trial is success or failure	outcomes have more than two categories (use the multinomial, Section 6.4)
Independent	trials do not affect each other	sampling without replacement from a small population (use the hypergeometric)
Number fixed	$n$ is set in advance	you stop at the first success (use the geometric)
Same $p$	every trial has the same success probability	users differ in their click probabilities (Poisson binomial)

The PMF. Any particular sequence with $k$ successes and $n - k$ failures has probability $p^k (1 - p)^{n - k}$, by independence. There are $\binom{n}{k}$ such sequences (choose the positions of the successes, Section 1.3). Therefore

$$P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k}, \qquad k = 0, 1, \dots, n.$$

  n = 3 trials, k = 2 successes:  C(3, 2) = 3 sequences, each with probability p^2 (1 - p)

     S S F   ->  p * p * (1-p)
     S F S   ->  p * (1-p) * p        P(X = 2) = 3 p^2 (1 - p)
     F S S   ->  (1-p) * p * p

The binomial theorem shows the probabilities sum to 1: $\sum_k \binom{n}{k} p^k (1 - p)^{n - k} = (p + (1 - p))^n = 1$.

Mean and variance. Write $X = X_1 + \dots + X_n$ with independent $X_i \sim$ Bernoulli$(p)$. Then

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

Shape. The distribution is symmetric for $p = 1/2$, right-skewed for $p < 1/2$ and left-skewed for $p > 1/2$. Its mode is $\lfloor (n + 1)p \rfloor$, and as $n$ grows it looks more and more like a normal curve (Section 5.3). Sums of independent binomials with the same $p$ are binomial: Binomial$(n_1, p)$ + Binomial$(n_2, p)$ = Binomial$(n_1 + n_2, p)$.

Example. A student guesses all 10 questions of a multiple-choice quiz with 4 options each. The number of correct answers is $X \sim$ Binomial$(10, 0.25)$.

import numpy as np
from math import comb

n, p = 10, 0.25
pmf = [comb(n, k) * p ** k * (1 - p) ** (n - k) for k in range(n + 1)]

print(" k   P(X = k)")
for k, pk in enumerate(pmf):
    print(f"{k:2d}   {pk:.4f}   {'#' * round(pk * 100)}")      # 1 '#' per percentage point

print(f"P(X >= 6) = {sum(pmf[6:]):.4f}   (passing by pure guessing)")
print(f"mean = {sum(k * pk for k, pk in enumerate(pmf)):.2f} (np = {n * p}), "
      f"variance = {sum(k * k * pk for k, pk in enumerate(pmf)) - (n * p) ** 2:.4f} (np(1-p) = {n * p * (1 - p)})")

# Simulation: a binomial count is a sum of n Bernoulli trials
rng = np.random.default_rng(2)
sums = (rng.random((100_000, n)) < p).sum(axis=1)
direct = rng.binomial(n, p, size=100_000)
print(f"P(X >= 6): sum of Bernoullis = {np.mean(sums >= 6):.4f}, rng.binomial = {np.mean(direct >= 6):.4f}")

 k   P(X = k)
 0   0.0563   ######
 1   0.1877   ###################
 2   0.2816   ############################
 3   0.2503   #########################
 4   0.1460   ###############
 5   0.0584   ######
 6   0.0162   ##
 7   0.0031
 8   0.0004
 9   0.0000
10   0.0000
P(X >= 6) = 0.0197   (passing by pure guessing)
mean = 2.50 (np = 2.5), variance = 1.8750 (np(1-p) = 1.875)
P(X >= 6): sum of Bernoullis = 0.0199, rng.binomial = 0.0197

Key takeaways

• Binomial$(n, p)$ counts successes in $n$ independent trials with the same success probability $p$ (BINS).
• $P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k}$: the number of arrangements times the probability of one arrangement.
• $E[X] = np$ and $\operatorname{Var}(X) = np(1 - p)$, from writing $X$ as a sum of Bernoulli indicators.
• Sampling without replacement, a random number of trials, or varying $p$ all break the binomial model.

🧪 Practice

1. A fair coin is tossed 8 times. Compute $P(\text{exactly 4 heads})$ and $P(\text{at least 6 heads})$.
2. Show that $\frac{P(X = k + 1)}{P(X = k)} = \frac{(n - k),p}{(k + 1)(1 - p)}$, and use it to prove that the mode is $\lfloor (n + 1)p \rfloor$.
3. Each of 20 independent servers is down on a given day with probability 0.02. Find the probability that at least one is down, and the expected number that are down.
4. Interview: An A/B test variant has a historical conversion rate of 5%. Out of 1,000 visitors, 65 converted. Is that surprising? Hint: compute the mean and standard deviation of Binomial$(1000, 0.05)$ and express 65 in standard deviations from the mean; then say what else could explain the result besides a real improvement.

Geometric

Why it matters. Instead of fixing the number of trials and counting successes, you often wait for the first success: retries until a request goes through, sales calls until the first sale, die rolls until the first 6. The geometric distribution models this waiting time. It is the discrete counterpart of the exponential distribution (Section 5.2) and the only discrete distribution without memory.

Definition. Run independent Bernoulli$(p)$ trials and let $K$ be the number of trials up to and including the first success. The first $k - 1$ trials must fail and the $k$-th must succeed:

$$P(K = k) = (1 - p)^{k - 1} p, \qquad k = 1, 2, 3, \dots$$

The tail is especially simple: $K > k$ means the first $k$ trials all failed, so $P(K > k) = (1 - p)^k$.

Two conventions. Some authors and libraries count the failures before the first success, $Y = K - 1 \in {0, 1, 2, \dots}$. Always check which one is meant.

Convention	Support	PMF	Mean	Variance	Software
Trials until first success $K$	$1, 2, \dots$	$(1 - p)^{k-1} p$	$\frac{1}{p}$	$\frac{1 - p}{p^2}$	numpy rng.geometric, scipy geom
Failures before first success $Y$	$0, 1, \dots$	$(1 - p)^{y} p$	$\frac{1 - p}{p}$	$\frac{1 - p}{p^2}$	R dgeom

Mean by first-step analysis (Section 3.4). The first trial succeeds with probability $p$ and the wait is 1. Otherwise one trial is used up and the process starts afresh:

$$E[K] = p \cdot 1 + (1 - p)(1 + E[K]) \quad \Longrightarrow \quad E[K] = \frac{1}{p}.$$

A 6 on a die needs 6 rolls on average, and an event with probability 1% needs 100 attempts on average.

Memorylessness. For all $m, n \ge 0$,

$$P(K > m + n \mid K > m) = \frac{(1 - p)^{m + n}}{(1 - p)^m} = (1 - p)^n = P(K > n).$$

Having already failed $m$ times tells you nothing about how much longer you will wait. The die does not "owe" you a 6. Believing otherwise is the gambler's fallacy. The geometric is the only distribution on the positive integers with this property.

import numpy as np

rng = np.random.default_rng(3)
p = 1 / 6                                          # waiting for a 6 on a fair die
k = rng.geometric(p, size=1_000_000)               # numpy counts trials: support 1, 2, 3, ...

print(f"mean = {k.mean():.3f} (1/p = 6), variance = {k.var():.2f} ((1-p)/p^2 = {(1 - p) / p ** 2:.0f})")
print(f"P(K > 10) = {np.mean(k > 10):.4f}   exact (5/6)^10 = {(5 / 6) ** 10:.4f}")

# Memorylessness: after 5 failed rolls, the remaining wait has the original distribution
extra = k[k > 5] - 5
print(" j   P(K = j)   P(K - 5 = j | K > 5)")
for j in range(1, 6):
    print(f" {j}   {np.mean(k == j):.4f}     {np.mean(extra == j):.4f}")
print(f"mean of remaining wait after 5 failures = {extra.mean():.3f}")

mean = 6.004 (1/p = 6), variance = 29.97 ((1-p)/p^2 = 30)
P(K > 10) = 0.1616   exact (5/6)^10 = 0.1615
 j   P(K = j)   P(K - 5 = j | K > 5)
 1   0.1662     0.1671
 2   0.1388     0.1395
 3   0.1156     0.1156
 4   0.0967     0.0964
 5   0.0799     0.0801
mean of remaining wait after 5 failures = 5.990

Key takeaways

• Geometric$(p)$ is the number of independent trials until the first success: $P(K = k) = (1 - p)^{k - 1}p$, $P(K > k) = (1 - p)^k$.
• $E[K] = 1/p$, $\operatorname{Var}(K) = (1 - p)/p^2$.
• Two conventions exist (trials vs. failures); they differ by 1 in the support and the mean.
• It is memoryless: past failures do not shorten the remaining wait.

🧪 Practice

1. A die is rolled until the first 6. Find the probability that more than 10 rolls are needed, and the probability that the number of rolls is even.
2. A request succeeds with probability 0.2 per attempt, independently. What is the smallest number of attempts that gives at least a 99% chance of success?
3. Use the tail formula $E[K] = \sum_{k \ge 0} P(K > k)$ (Section 4.3) to rederive $E[K] = 1/p$.
4. Interview: (Coupon collector) How many rolls of a fair die do you need, on average, to see all six faces at least once? Hint: split the wait into stages: once you have seen $i$ distinct faces, the wait for a new one is geometric with success probability $(6 - i)/6$; add the expected stage lengths.

Negative Binomial

Why it matters. The negative binomial generalizes the geometric to waiting for the $r$-th success: how many games a team plays to win 4, how many calls a salesperson makes to close 3 deals. It also has a second life as the standard model for counts that vary more than the Poisson allows, such as the number of claims per customer or reads per gene in sequencing data, and it is the basis of negative binomial regression (Section 10.5).

The story. Run independent Bernoulli$(p)$ trials until the $r$-th success, and let $K$ be the number of trials. The event ${K = k}$ means that the $k$-th trial is a success and the first $k - 1$ trials contain exactly $r - 1$ successes:

$$P(K = k) = \binom{k - 1}{r - 1} p^r (1 - p)^{k - r}, \qquad k = r, r + 1, \dots$$

  r = 3: wait for the 3rd success

   trial:   1   2   3   4   5   6
            F   S   F   F   S   S      K = 6
            '---------------'   ^
            2 successes among   the 6th trial must be
            the first 5 trials  the 3rd success

The waiting time for the $r$-th success is the sum of $r$ independent geometric waits (one per success), so

$$E[K] = \frac{r}{p}, \qquad \operatorname{Var}(K) = \frac{r(1 - p)}{p^2}.$$

As with the geometric, many libraries count the failures $Y = K - r$ instead, with $P(Y = y) = \binom{y + r - 1}{y} p^r (1 - p)^y$ and mean $r(1 - p)/p$. numpy's rng.negative_binomial(r, p) returns failures.

Second life: overdispersed counts. A Poisson model (later in this section) forces the variance of a count to equal its mean. Real counts are often overdispersed: some customers file many more claims than others because their underlying rates differ. If each unit has its own Poisson rate $\Lambda$, and the rates vary across units according to a gamma distribution (Section 5.2), the resulting count is negative binomial. In this parameterization by the mean $\mu$ and a dispersion parameter $r > 0$ (not necessarily an integer),

$$E[Y] = \mu, \qquad \operatorname{Var}(Y) = \mu + \frac{\mu^2}{r}.$$

Small $r$ means strong overdispersion; as $r \to \infty$ the negative binomial tends to the Poisson.

	Waiting-time view	Count-model view
Parameters	$r$ successes needed (integer), success probability $p$	mean $\mu$, dispersion $r > 0$
Variable	trials (or failures) until the $r$-th success	count per unit
Typical use	sequences of games, sales, retries	claims, defects, website visits, gene counts
Relation	$p = r/(r + \mu)$	$\mu = r(1 - p)/p$

import numpy as np
from math import comb, exp, factorial

# 1. Waiting view: a player who scores 70% of free throws practices until the 3rd score
r, p = 3, 0.7
print(" k   P(K = k)")
for k in range(r, r + 6):
    print(f" {k}   {comb(k - 1, r - 1) * p ** r * (1 - p) ** (k - r):.4f}")
rng = np.random.default_rng(4)
trials = rng.negative_binomial(r, p, size=500_000) + r     # numpy returns failures; add the r successes
print(f"mean trials = {trials.mean():.3f} (r/p = {r / p:.3f})")

# 2. Count view: Poisson(4) vs negative binomial with mean 4 and dispersion r = 2
mu, r2 = 4.0, 2.0
rates = rng.gamma(shape=r2, scale=mu / r2, size=500_000)   # each unit gets its own rate (mean mu)
nb = rng.poisson(rates)                                    # gamma-Poisson mixture
nb_direct = rng.negative_binomial(r2, r2 / (r2 + mu), size=500_000)
pois = rng.poisson(mu, size=500_000)
for name, x in [("Poisson(4)", pois), ("gamma-Poisson mix", nb), ("negative binomial", nb_direct)]:
    print(f"{name:18s} mean = {x.mean():.3f}  variance = {x.var():6.3f}  "
          f"P(0) = {np.mean(x == 0):.3f}  P(>= 10) = {np.mean(x >= 10):.4f}")
print(f"theory: NB variance = mu + mu^2/r = {mu + mu ** 2 / r2:.0f}, NB P(0) = (r/(r+mu))^r = {(r2 / (r2 + mu)) ** r2:.3f}, "
      f"Poisson P(0) = {exp(-mu):.3f}")

 k   P(K = k)
 3   0.3430
 4   0.3087
 5   0.1852
 6   0.0926
 7   0.0417
 8   0.0175
mean trials = 4.287 (r/p = 4.286)
Poisson(4)         mean = 4.000  variance =  3.995  P(0) = 0.018  P(>= 10) = 0.0080
gamma-Poisson mix  mean = 3.998  variance = 11.993  P(0) = 0.111  P(>= 10) = 0.0754
negative binomial  mean = 4.000  variance = 11.978  P(0) = 0.111  P(>= 10) = 0.0747
theory: NB variance = mu + mu^2/r = 12, NB P(0) = (r/(r+mu))^r = 0.111, Poisson P(0) = 0.018

With the same mean, the negative binomial has many more zeros and many more large counts than the Poisson. A model that ignores overdispersion underestimates both.

Key takeaways

• The negative binomial counts the trials (or failures) until the $r$-th success; it is a sum of $r$ independent geometrics.
• $P(K = k) = \binom{k-1}{r-1} p^r (1-p)^{k-r}$, $E[K] = r/p$, $\operatorname{Var}(K) = r(1-p)/p^2$.
• As a count model it is a gamma mixture of Poissons, with $\operatorname{Var} = \mu + \mu^2/r > \mu$.
• Use it when count data are overdispersed (variance well above the mean).

🧪 Practice

1. A fair coin is tossed until the third head. Find the probability that exactly 5 tosses are needed.
2. In a best-of-7 series, a team wins each game independently with probability 0.6. Use the negative binomial to find the probability that it wins the series (4 wins before the opponent gets 4).
3. Show that the mean and variance of the waiting time for the $r$-th success follow from those of the geometric distribution.
4. Interview: The number of support tickets per customer per month has sample mean 3 and sample variance 12. A colleague fits a Poisson regression. What is wrong, and what would you use? Hint: compare the variance with the mean, solve $\mu + \mu^2/r = 12$ for the dispersion, and think about what the Poisson fit does to standard errors and p-values.

Hypergeometric

Why it matters. The binomial assumes independent trials, which in sampling means drawing with replacement. Real sampling is usually done without replacement: you do not inspect the same part twice, deal the same card twice or survey the same person twice. When the sample is a sizeable fraction of the population, the hypergeometric distribution gives the exact answer. It is used in acceptance sampling, card games and lotteries, capture-recapture estimates of population size, and Fisher's exact test (Section 9.3).

The story. A population of $N$ items contains $K$ "successes" (for example defective parts). Draw $n$ items without replacement, with all subsets of size $n$ equally likely, and let $X$ be the number of successes in the sample. By counting (Section 3.2):

$$P(X = k) = \frac{\binom{K}{k}\binom{N - K}{n - k}}{\binom{N}{n}}, \qquad \max(0, n - N + K) \le k \le \min(n, K).$$

The numerator chooses $k$ of the $K$ successes and $n - k$ of the $N - K$ failures; the denominator counts all possible samples.

  urn: N = 20 balls, K = 6 red (R), 14 blue (b)       draw n = 5 without putting back

   R R R R R R b b b b b b b b b b b b b b    ->   sample: R b b R b    X = 2

Mean and variance. With $p = K/N$, the proportion of successes in the population,

$$E[X] = np, \qquad \operatorname{Var}(X) = np(1 - p),\frac{N - n}{N - 1}.$$

The mean matches the binomial. The variance is smaller by the finite population correction $\frac{N - n}{N - 1}$ (Section 7.3): each draw removes an item, so the composition of what is left becomes more predictable. If you sample the whole population ($n = N$), the correction is 0 and $X = K$ with certainty.

	Binomial$(n, p)$	Hypergeometric$(N, K, n)$
Sampling	with replacement, or independent trials	without replacement
Draws independent?	yes	no: each draw changes the urn
Mean	$np$	$np$ with $p = K/N$
Variance	$np(1 - p)$	$np(1 - p)\frac{N - n}{N - 1}$
Use when	$p$ fixed, or population huge	finite population, $n/N$ not small

When the binomial is close enough. If the sample is small relative to the population (a common rule of thumb is $n/N < 0.05$), removing items barely changes the proportions, and Binomial$(n, K/N)$ is an excellent approximation. Polling 1,000 people out of millions is binomial in practice.

import numpy as np
from math import comb

def hypergeom_pmf(k, N, K, n):
    return comb(K, k) * comb(N - K, n - k) / comb(N, n)

# 1. Acceptance sampling: a lot of 100 parts has 5 defective. Inspect 10, accept the lot if none is defective
N, K, n = 100, 5, 10
p_accept = hypergeom_pmf(0, N, K, n)
print(f"P(accept) exact = {p_accept:.4f}, binomial approximation (0.95)^10 = {0.95 ** 10:.4f}")
rng = np.random.default_rng(6)
x = rng.hypergeometric(ngood=K, nbad=N - K, nsample=n, size=500_000)
p = K / N
print(f"mean = {x.mean():.4f} (np = {n * p}), variance = {x.var():.4f} "
      f"(hypergeometric {n * p * (1 - p) * (N - n) / (N - 1):.4f}, binomial {n * p * (1 - p):.4f})")

# 2. Lottery 6/49: 6 winning numbers among 49, a ticket picks 6
print(" matches   probability     1 in")
for k in range(7):
    pk = hypergeom_pmf(k, 49, 6, 6)
    print(f"    {k}      {pk:.3e}   {1 / pk:12,.0f}")

P(accept) exact = 0.5838, binomial approximation (0.95)^10 = 0.5987
mean = 0.5002 (np = 0.5), variance = 0.4320 (hypergeometric 0.4318, binomial 0.4750)
 matches   probability     1 in
    0      4.360e-01              2
    1      4.130e-01              2
    2      1.324e-01              8
    3      1.765e-02             57
    4      9.686e-04          1,032
    5      1.845e-05         54,201
    6      7.151e-08     13,983,816

Key takeaways

• The hypergeometric counts successes in a sample drawn without replacement from a finite population.
• $P(X = k) = \binom{K}{k}\binom{N-K}{n-k} / \binom{N}{n}$.
• Mean $nK/N$, like the binomial; variance smaller by the finite population correction $(N - n)/(N - 1)$.
• When the sampling fraction $n/N$ is small, the binomial approximation is accurate.

🧪 Practice

1. Find the probability that a 5-card poker hand contains exactly 2 aces.
2. (Capture-recapture) 50 fish are tagged and released. Later, 40 fish are caught and 8 of them are tagged. Estimate the number of fish in the lake by setting the expected proportion of tagged fish equal to the observed one, and write down the hypergeometric probability of the observation as a function of the unknown $N$.
3. Compute the finite population correction for $N = 10{,}000$ and $n = 100$, and for $N = 200$ and $n = 100$. What does each tell you about using the binomial?
4. Interview: A code review samples 10 of the 200 files changed in a release, and 15 of the 200 files contain a bug. What is the probability that the review finds no buggy file, and how many files would you need to review to catch at least one with 95% probability? Hint: $P(X = 0)$ is a ratio of two binomial coefficients; increase $n$ until it drops below 0.05, and compare with the answer the binomial would give.

Poisson

Why it matters. Many counts record how often something happens in a fixed amount of time or space, rather than how many of $n$ trials succeeded: requests per second at a server, calls per hour at a call center, typos per page, accidents per month at an intersection, radioactive decays per second. When events occur independently of each other and at a steady average rate, the count follows a Poisson distribution. It is the default model for counts and the basis of Poisson regression (Section 10.5).

Definition. $X \sim$ Poisson$(\lambda)$, with rate $\lambda > 0$, if

$$P(X = k) = e^{-\lambda} \frac{\lambda^k}{k!}, \qquad k = 0, 1, 2, \dots$$

The probabilities sum to 1 because $\sum_k \lambda^k/k! = e^{\lambda}$. The parameter $\lambda$ is the expected count in the interval considered. If events occur at rate $\nu$ per unit of time, the count in an interval of length $t$ is Poisson$(\nu t)$: 3 requests per second means Poisson$(180)$ per minute.

Mean equals variance.

$$E[X] = \lambda, \qquad \operatorname{Var}(X) = \lambda.$$

This is the Poisson's signature. The dispersion index $\operatorname{Var}/\text{mean}$ of real count data is a quick check: a value well above 1 (overdispersion) points to the negative binomial, and a value well below 1 (underdispersion) points to events that are more regular than random, like scheduled arrivals.

Analogy. Rain starting to fall on a pavement: each drop lands independently of the others, at a steady rate. The number of drops on one paving stone in one second is Poisson. It does not matter that the cloud holds an astronomically large number of drops, each with a tiny probability of hitting that stone; the next topic shows why.

Properties.

Property	Statement
Additivity	independent Poisson$(\lambda_1)$ + Poisson$(\lambda_2)$ = Poisson$(\lambda_1 + \lambda_2)$ (proof by MGFs, Section 4.3)
Thinning	if each event is kept independently with probability $q$, the kept count is Poisson$(q\lambda)$, and the kept and discarded counts are independent
Recurrence	$P(X = k + 1) = P(X = k) \cdot \frac{\lambda}{k + 1}$, which avoids huge factorials in computation
Mode	$\lfloor \lambda \rfloor$ (both $\lambda - 1$ and $\lambda$ when $\lambda$ is an integer)
Shape	right-skewed for small $\lambda$ (skewness $1/\sqrt{\lambda}$), close to normal for large $\lambda$

import numpy as np
from math import exp

lam = 3.2                                       # requests per second, on average

# PMF by the recurrence p(k+1) = p(k) * lam / (k+1), starting from p(0) = e^(-lam)
pmf = [exp(-lam)]
for k in range(30):
    pmf.append(pmf[-1] * lam / (k + 1))
print(" k   P(X = k)")
for k in range(10):
    print(f"{k:2d}   {pmf[k]:.4f}   {'#' * round(pmf[k] * 100)}")
print(f"P(X >= 8) = {1 - sum(pmf[:8]):.4f}")

# Capacity planning: smallest c with P(X > c) < 0.001
tail, c = 1.0, -1
while tail >= 0.001:
    c += 1
    tail -= pmf[c]
print(f"capacity of {c} requests/s is exceeded with probability {tail:.5f}")

# Mean = variance, and thinning: 25% of requests are writes
rng = np.random.default_rng(7)
x = rng.poisson(lam, size=500_000)
writes = rng.binomial(x, 0.25)                  # each request is a write with probability 0.25
reads = x - writes
print(f"all: mean = {x.mean():.3f}, var = {x.var():.3f};  writes: mean = {writes.mean():.3f}, "
      f"var = {writes.var():.3f};  corr(writes, reads) = {np.corrcoef(writes, reads)[0, 1]:+.4f}")

 k   P(X = k)
 0   0.0408   ####
 1   0.1304   #############
 2   0.2087   #####################
 3   0.2226   ######################
 4   0.1781   ##################
 5   0.1140   ###########
 6   0.0608   ######
 7   0.0278   ###
 8   0.0111   #
 9   0.0040
P(X >= 8) = 0.0168
capacity of 10 requests/s is exceeded with probability 0.00050
all: mean = 3.199, var = 3.196;  writes: mean = 0.801, var = 0.802;  corr(writes, reads) = -0.0024

The thinned counts are Poisson with rate $0.25 \times 3.2 = 0.8$ (mean equals variance), and the numbers of writes and reads are uncorrelated even though they come from the same stream. This is a special property of the Poisson: if the total were fixed, more writes would mean fewer reads.

Key takeaways

• Poisson$(\lambda)$ models counts of independent events occurring at a constant average rate: $P(X = k) = e^{-\lambda}\lambda^k/k!$.
• Mean and variance both equal $\lambda$; compare them in data to detect over- or underdispersion.
• The rate scales with the length of the interval, independent Poisson counts add, and random thinning keeps a count Poisson.
• Compute probabilities with the recurrence or a library, not with raw factorials.

🧪 Practice

1. A book has on average 2 typos per page. Find the probability that a given page has no typos, and that 3 consecutive pages have no typos at all.
2. A shop receives on average 10 emails per hour, and 20% of emails are spam. What is the distribution of the number of spam emails in 2 hours? What is the probability of more than 6?
3. Show directly from the PMF that $E[X] = \lambda$, and use $E[X(X - 1)] = \lambda^2$ to get the variance.
4. Interview: A service receives 1 million requests per day. How would you set an alert threshold on requests per minute that fires rarely under normal conditions? Hint: compute the Poisson rate per minute and its standard deviation, then question the assumptions: traffic varies over the day, bursts make counts overdispersed, so the threshold should depend on the time of day and on observed variance, not only on $\lambda$.

Poisson Approximation to the Binomial

Why it matters. Binomial probabilities with large $n$ and small $p$ are common (defects among millions of chips, fraud among millions of transactions, users who hit a rare bug) and awkward to compute directly. The Poisson approximation replaces them with a one-parameter formula. More importantly, it explains why the Poisson distribution is everywhere: a count of rare events among many independent opportunities is approximately Poisson. This is the law of rare events.

Statement. If $n \to \infty$ and $p \to 0$ with $np \to \lambda$, then for each fixed $k$

$$\binom{n}{k} p^k (1 - p)^{n - k} ;\longrightarrow; e^{-\lambda}\frac{\lambda^k}{k!}.$$

Derivation. Put $p = \lambda/n$ and regroup the binomial PMF:

$$\binom{n}{k}\left(\frac{\lambda}{n}\right)^k\left(1 - \frac{\lambda}{n}\right)^{n-k} = \underbrace{\frac{n(n-1)\cdots(n-k+1)}{n^k}}{\to,1};\cdot;\frac{\lambda^k}{k!};\cdot;\underbrace{\left(1 - \frac{\lambda}{n}\right)^{n}}{\to,e^{-\lambda}};\cdot;\underbrace{\left(1 - \frac{\lambda}{n}\right)^{-k}}_{\to,1}.$$

How good is it? A clean bound, due to Le Cam, measures the error with the total variation distance, the largest difference between the two distributions' probabilities of any event:

$$d_{TV} = \max_A \left|P(X \in A) - P(Y \in A)\right| = \frac12 \sum_k \left|P(X = k) - P(Y = k)\right| ;\le; \sum_{i=1}^n p_i^2.$$

The result does not require equal probabilities. For a sum of independent Bernoulli$(p_i)$ variables (the Poisson binomial) compared with Poisson$(\sum p_i)$, the error is at most $\sum p_i^2$, which is $np^2 = \lambda^2/n$ when all $p_i$ are equal. The approximation is good when every individual probability is small, however many trials there are.

Rule of thumb	Use
$n \ge 20$ and $p \le 0.05$	Poisson$(np)$ is usually accurate
$n \ge 100$ and $np \le 10$	very accurate
$p$ not small, $np(1 - p) \ge 10$	prefer the normal approximation (Section 5.3)

Note that the binomial variance $np(1 - p)$ is slightly less than the Poisson variance $np$; the factor $1 - p$ is what the approximation ignores.

import numpy as np
from math import comb, exp, factorial

def poisson_pmf(lam, kmax):
    return np.array([exp(-lam) * lam ** k / factorial(k) for k in range(kmax + 1)])

def tv_distance(p, q):
    return 0.5 * np.abs(p - q).sum()

lam, kmax = 3.0, 60
pois = poisson_pmf(lam, kmax)
print("    n       p     d_TV(Binomial, Poisson)   Le Cam bound n p^2")
for n in [10, 30, 100, 1000, 10000]:
    p = lam / n
    binom = np.array([comb(n, k) * p ** k * (1 - p) ** (n - k) if k <= n else 0.0 for k in range(kmax + 1)])
    print(f"{n:6d}  {p:.4f}   {tv_distance(binom, pois):.5f}                   {n * p * p:.5f}")

# Unequal probabilities: 2,000 users, each with their own small chance of hitting a bug today
rng = np.random.default_rng(8)
p_i = rng.uniform(0, 0.003, size=2000)
pmf = np.zeros(kmax + 1); pmf[0] = 1.0
for q in p_i:                                   # exact PMF of the sum, one Bernoulli at a time
    pmf[1:] = pmf[1:] * (1 - q) + pmf[:-1] * q
    pmf[0] *= 1 - q
print(f"Poisson binomial vs Poisson({p_i.sum():.3f}): d_TV = {tv_distance(pmf, poisson_pmf(p_i.sum(), kmax)):.5f}, "
      f"bound sum p_i^2 = {np.sum(p_i ** 2):.5f}")

    n       p     d_TV(Binomial, Poisson)   Le Cam bound n p^2
    10  0.3000   0.08641                   0.90000
    30  0.1000   0.02618                   0.30000
   100  0.0300   0.00761                   0.09000
  1000  0.0030   0.00075                   0.00900
 10000  0.0003   0.00008                   0.00090
Poisson binomial vs Poisson(2.999): d_TV = 0.00049, bound sum p_i^2 = 0.00588

Key takeaways

• For large $n$ and small $p$, Binomial$(n, p) \approx$ Poisson$(np)$.
• The error in total variation is at most $\sum p_i^2$ ($= np^2$ for equal $p$), so what matters is that every $p_i$ is small.
• The approximation also covers trials with unequal probabilities, which is why counts of rare events are Poisson in practice.
• When $p$ is not small, use the exact binomial or the normal approximation instead.

🧪 Practice

1. A 500-page book contains 300 typos, each independently on a uniformly random page. Find the probability that a given page has at least 2 typos, exactly (binomial) and approximately (Poisson).
2. Compute the total variation distance between Binomial$(20, 0.3)$ and Poisson$(6)$ by adapting the code. Is the approximation acceptable? Compare the two variances.
3. Prove that $\left(1 - \frac{\lambda}{n}\right)^n \to e^{-\lambda}$ by taking logarithms and using $\ln(1 - x) = -x + O(x^2)$.
4. Interview: Website visitors are not identical: each has a different probability of buying, and the number of potential visitors is finite. Why is it still reasonable to model daily purchases as Poisson, and when would it fail? Hint: use Le Cam's bound for unequal $p_i$; then consider what happens if a few visitors have large $p_i$ (wholesale buyers) or if purchases are not independent (a viral post, a sale).

<a id="52-continuous-distributions"></a>

5.2 Continuous Distributions

Continuous distributions model measurements such as durations, sizes, amounts and proportions. As in Section 5.1, each family comes from a mechanism, and the mechanism tells you which data it fits and how its tails behave.

Continuous Uniform

Why it matters. The continuous uniform distribution says "any value in this interval, none preferred". It describes the arrival time within a cycle when you show up at a random moment, the error made by rounding a measurement, and a random angle. Above all, Uniform$(0, 1)$ is what random number generators produce: inverse transform sampling (Section 4.4) turns it into any other distribution, so every simulation in this guide starts from it.

Definition. $X \sim$ Uniform$(a, b)$, with $a < b$, if its density is constant on $[a, b]$:

$$f(x) = \frac{1}{b - a} ;; (a \le x \le b), \qquad F(x) = \frac{x - a}{b - a} ;; (a \le x \le b), \qquad Q(p) = a + p,(b - a).$$

Probability is proportional to length: $P(c \le X \le d) = \frac{d - c}{b - a}$ for any $[c, d] \subseteq [a, b]$.

  f(x)
  1/(b-a) |      +---------------------+
          |      |                     |
          |      |   area = 1          |
        0 +------+---------------------+------> x
                 a                     b

Moments.

$$E[X] = \frac{a + b}{2}, \qquad \operatorname{Var}(X) = \frac{(b - a)^2}{12}, \qquad \operatorname{SD}(X) = \frac{b - a}{\sqrt{12}} \approx 0.289,(b - a).$$

Uniform$(a, b)$ is a location-scale family (Section 4.4): $X = a + (b - a)U$ with $U \sim$ Uniform$(0, 1)$, so all results follow from those for $U$.

Rounding errors. If a measurement is rounded to the nearest unit and the true value is spread over many units, the rounding error is approximately Uniform$(-0.5, 0.5)$, with standard deviation $1/\sqrt{12} \approx 0.29$ units. This is the standard way to account for the resolution of an instrument or the quantization error of a digital signal.

Sums of uniforms are not uniform. The sum of two independent Uniform$(0, 1)$ variables has a triangular density on $[0, 2]$, peaked at 1. The sum of $n$ of them (the Irwin-Hall distribution) quickly looks normal. Old software generated approximate standard normals as $U_1 + \dots + U_{12} - 6$, which has mean 0 and variance $12 \cdot \frac{1}{12} = 1$, but which can never exceed 6 in absolute value.

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(9)

# 1. A bus comes every 10 minutes; you arrive at a uniformly random moment
wait = rng.uniform(0, 10, size=1_000_000)
print(f"wait: mean = {wait.mean():.3f} (5), variance = {wait.var():.3f} (100/12 = {100 / 12:.3f}), "
      f"P(wait > 7) = {np.mean(wait > 7):.4f} (0.3)")

# 2. Rounding error of heights measured to the nearest cm
true = rng.normal(170, 8, size=1_000_000)
err = np.round(true) - true
print(f"rounding error: range [{err.min():.3f}, {err.max():.3f}], mean = {err.mean():+.4f}, "
      f"SD = {err.std():.4f} (1/sqrt(12) = {1 / np.sqrt(12):.4f})")

# 3. Twelve uniforms minus 6: close to N(0, 1) in the middle, but the tails are cut off
s = rng.random((1_000_000, 12)).sum(axis=1) - 6
print(f"sum of 12 uniforms - 6: mean = {s.mean():+.4f}, variance = {s.var():.4f}")
for z in [1, 2, 3, 4]:
    print(f"  P(S > {z}) = {np.mean(s > z):.2e}   normal: {1 - NormalDist().cdf(z):.2e}")

wait: mean = 4.996 (5), variance = 8.343 (100/12 = 8.333), P(wait > 7) = 0.2992 (0.3)
rounding error: range [-0.500, 0.500], mean = +0.0000, SD = 0.2885 (1/sqrt(12) = 0.2887)
sum of 12 uniforms - 6: mean = -0.0001, variance = 1.0002
  P(S > 1) = 1.61e-01   normal: 1.59e-01
  P(S > 2) = 2.23e-02   normal: 2.28e-02
  P(S > 3) = 9.58e-04   normal: 1.35e-03
  P(S > 4) = 9.00e-06   normal: 3.17e-05

The approximation is excellent within 2 standard deviations and poor in the far tail, a pattern that recurs whenever a distribution is approximated by the normal.

Key takeaways

• Uniform$(a, b)$ has constant density $1/(b - a)$; probabilities are proportional to lengths.
• Mean $(a + b)/2$, variance $(b - a)^2/12$.
• Uniform$(0, 1)$ is the raw material of all simulation; rounding errors are approximately uniform.
• Sums of uniforms are not uniform: they become bell-shaped, with bounded tails.

🧪 Practice

1. $X \sim$ Uniform$(2, 8)$. Find $P(|X - 5| < 1)$, the mean, the variance and the 90th percentile.
2. Derive $\operatorname{Var}(X) = (b - a)^2/12$ by computing $E[X^2]$ for $X \sim$ Uniform$(0, 1)$ and using the location-scale relation.
3. $U_1, U_2$ are independent Uniform$(0, 1)$. Use the CDF method (Section 4.4) to find the distribution of $\max(U_1, U_2)$ and of $\min(U_1, U_2)$.
4. Interview: Break a stick of length 1 at two points chosen independently and uniformly. What is the probability that the three pieces can form a triangle? Hint: the pieces form a triangle exactly when every piece is shorter than 1/2; draw the region in the unit square and compute its area, then check by simulation.

Exponential and Memorylessness

Why it matters. The exponential distribution is the standard model for waiting times when events occur at random at a constant rate: the time until the next request arrives, until a radioactive atom decays, until a component without wear fails. It is the continuous counterpart of the geometric distribution, and its defining feature, memorylessness, is either exactly right or badly wrong for a given application. Knowing which is essential in reliability, queueing and survival analysis (Section 12.3).

Definition. $X \sim$ Exponential$(\lambda)$, with rate $\lambda > 0$, if

$$f(x) = \lambda e^{-\lambda x}, \qquad F(x) = 1 - e^{-\lambda x}, \qquad P(X > x) = e^{-\lambda x}, \qquad x \ge 0.$$

Quantity	Value
Mean	$1/\lambda$
Variance	$1/\lambda^2$ (so SD = mean, coefficient of variation 1)
Median	$\ln 2/\lambda \approx 0.693/\lambda$, below the mean
Quantile	$Q(p) = -\ln(1 - p)/\lambda$ (Section 4.2)
$P(X \le \text{mean})$	$1 - e^{-1} \approx 0.632$

Rate or scale? The same distribution is often parameterized by the mean $\beta = 1/\lambda$ (the scale). numpy's rng.exponential(scale) and Python's random.expovariate(lambd) use opposite conventions: always check.

Memorylessness. For all $s, t \ge 0$,

$$P(X > s + t \mid X > s) = \frac{e^{-\lambda(s + t)}}{e^{-\lambda s}} = e^{-\lambda t} = P(X > t).$$

An item that has survived $s$ hours is, statistically, as good as new. Equivalently, its hazard rate, the instantaneous failure rate among survivors $h(t) = f(t)/(1 - F(t))$, is the constant $\lambda$. The exponential is the only continuous distribution with this property.

Analogy. A radioactive atom does not age. Whether it was created a second ago or a billion years ago, the chance that it decays in the next minute is the same. That is why decay is described by a half-life: the time by which half the atoms have decayed, the median $\ln 2/\lambda$.

Memorylessness is right for events triggered by external shocks at a constant rate (requests, calls, decays, accidents) and wrong for anything that wears out, such as tires, bearings or human lifetimes, where the Weibull (later in this section) is a better model.

Competing exponentials. If $X_1, \dots, X_n$ are independent with $X_i \sim$ Exponential$(\lambda_i)$, then

$$\min_i X_i \sim \text{Exponential}(\lambda_1 + \dots + \lambda_n), \qquad P(X_j \text{ is the smallest}) = \frac{\lambda_j}{\lambda_1 + \dots + \lambda_n}.$$

The proof is one line: $P(\min > t) = \prod_i e^{-\lambda_i t} = e^{-(\sum \lambda_i)t}$. A system that fails as soon as any of its parts fails has a failure rate equal to the sum of the parts' rates.

import numpy as np

rng = np.random.default_rng(10)

# Component lifetimes with mean 100 hours (rate 0.01 per hour); numpy uses the scale = mean
life = rng.exponential(scale=100, size=1_000_000)
print(f"mean = {life.mean():.2f}, SD = {life.std():.2f}, median = {np.median(life):.2f} "
      f"(100 ln 2 = {100 * np.log(2):.2f}), P(T < mean) = {np.mean(life < 100):.4f}")

# Memorylessness: components that survived 150 hours
rest = life[life > 150] - 150
print(f"exponential: remaining life after 150 h: mean = {rest.mean():.2f}, "
      f"P(> 50 more) = {np.mean(rest > 50):.4f} vs P(T > 50) = {np.mean(life > 50):.4f}")

# Contrast: a wearing-out lifetime, Uniform(0, 200), with the same mean of 100 hours
worn = rng.uniform(0, 200, size=1_000_000)
rest_w = worn[worn > 150] - 150
print(f"uniform:     remaining life after 150 h: mean = {rest_w.mean():.2f}, "
      f"P(> 50 more) = {np.mean(rest_w > 50):.4f} vs P(T > 50) = {np.mean(worn > 50):.4f}")

# Competing risks: disk (mean 100 h) and fan (mean 50 h); the machine stops at the first failure
disk, fan = rng.exponential(100, 1_000_000), rng.exponential(50, 1_000_000)
first = np.minimum(disk, fan)
print(f"time to first failure: mean = {first.mean():.2f} (1/(0.01 + 0.02) = {1 / 0.03:.2f}), "
      f"P(disk fails first) = {np.mean(disk < fan):.4f} (1/3)")

mean = 100.01, SD = 100.03, median = 69.43 (100 ln 2 = 69.31), P(T < mean) = 0.6319
exponential: remaining life after 150 h: mean = 99.93, P(> 50 more) = 0.6049 vs P(T > 50) = 0.6070
uniform:     remaining life after 150 h: mean = 25.03, P(> 50 more) = 0.0000 vs P(T > 50) = 0.7497
time to first failure: mean = 33.35 (1/(0.01 + 0.02) = 33.33), P(disk fails first) = 0.3330 (1/3)

Key takeaways

• Exponential$(\lambda)$: $P(X > x) = e^{-\lambda x}$, mean $1/\lambda$, SD equal to the mean.
• It is memoryless and has a constant hazard rate: the past does not change the future waiting time.
• The minimum of independent exponentials is exponential with the summed rate.
• It fits random shocks at a constant rate, not processes with wear-out; check whether software uses the rate or the scale.

🧪 Practice

1. A radioactive isotope has a half-life of 8 days. Find the decay rate $\lambda$, the mean lifetime of an atom, and the fraction of atoms remaining after 30 days.
2. Show that $P(X > 2E[X]) = e^{-2}$ for every exponential distribution. Why does the answer not depend on $\lambda$?
3. A server cluster has 20 independent machines, each failing after an exponential time with mean 1,000 hours. What is the distribution of the time until the first failure, and its mean?
4. Interview: Buses arrive at a stop as a Poisson process, with exponential gaps averaging 10 minutes. You arrive at a random time. What is your expected wait? Most people say 5 minutes; explain why it is 10. Hint: use memorylessness; then explain the inspection paradox: a random arrival is more likely to fall into a long gap than a short one, and compute the expected length of the gap you land in.

Gamma

Why it matters. Where the exponential is the waiting time for one event, the gamma distribution is the waiting time for several: the time until the 10th customer arrives, or until a component has absorbed 5 shocks. More generally, it is a flexible model for positive, right-skewed quantities such as rainfall amounts, insurance claims and service times. It contains the exponential and the chi-square (Section 5.4) as special cases, and it is the conjugate prior for a Poisson rate (Section 11.1).

The gamma function. The normalizing constant involves

$$\Gamma(\alpha) = \int_0^\infty x^{\alpha - 1} e^{-x},dx, \qquad \alpha > 0,$$

which extends the factorial to non-integers: $\Gamma(\alpha + 1) = \alpha,\Gamma(\alpha)$, $\Gamma(n) = (n - 1)!$ for positive integers, and $\Gamma(1/2) = \sqrt{\pi}$.

Definition. $X \sim$ Gamma$(\alpha, \beta)$, with shape $\alpha > 0$ and rate $\beta > 0$, if

$$f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)},x^{\alpha - 1} e^{-\beta x}, \qquad x > 0.$$

$$E[X] = \frac{\alpha}{\beta}, \qquad \operatorname{Var}(X) = \frac{\alpha}{\beta^2}, \qquad \text{mode} = \frac{\alpha - 1}{\beta} ;; (\alpha \ge 1), \qquad \text{skewness} = \frac{2}{\sqrt{\alpha}}.$$

The shape parameter controls the form of the density, and the rate only rescales it ($cX \sim$ Gamma$(\alpha, \beta/c)$). As with the exponential, many sources use the scale $\theta = 1/\beta$ instead; numpy's rng.gamma(shape, scale) does.

Shape $\alpha$	Form of the density	Note
$\alpha < 1$	infinite at 0, decreasing	many tiny values, a few large ones
$\alpha = 1$	$\beta e^{-\beta x}$	Exponential$(\beta)$
$\alpha = n$ integer	hump, right-skewed	Erlang: sum of $n$ exponentials
$\alpha = \nu/2$, $\beta = 1/2$	hump, right-skewed	chi-square with $\nu$ degrees of freedom
$\alpha$ large	nearly symmetric, nearly normal	skewness $2/\sqrt{\alpha} \to 0$

Additivity. If $X_1 \sim$ Gamma$(\alpha_1, \beta)$ and $X_2 \sim$ Gamma$(\alpha_2, \beta)$ are independent, with the same rate, then $X_1 + X_2 \sim$ Gamma$(\alpha_1 + \alpha_2, \beta)$. The MGF $\left(\frac{\beta}{\beta - t}\right)^\alpha$ makes this immediate (Section 4.3). In particular, the sum of $n$ independent Exponential$(\lambda)$ variables is Gamma$(n, \lambda)$.

Gamma and Poisson are two views of the same process. If events occur at random at rate $\lambda$, the time $T_n$ of the $n$-th event is Gamma$(n, \lambda)$, and the count $N(t)$ of events by time $t$ is Poisson$(\lambda t)$. "The $n$-th event happened by time $t$" is the same as "at least $n$ events happened by time $t$":

$$P(T_n \le t) = P(N(t) \ge n) = 1 - \sum_{k=0}^{n-1} e^{-\lambda t}\frac{(\lambda t)^k}{k!}.$$

Section 5.5 develops this connection.

import numpy as np
from math import exp, factorial

rng = np.random.default_rng(11)

# Support tickets arrive at 6 per hour. T = time until the 10th ticket ~ Gamma(10, rate 6)
n, lam, t = 10, 6.0, 1.5
via_poisson = 1 - sum(exp(-lam * t) * (lam * t) ** k / factorial(k) for k in range(n))
sum_of_exp = rng.exponential(1 / lam, size=(500_000, n)).sum(axis=1)   # 10 exponential gaps
direct = rng.gamma(shape=n, scale=1 / lam, size=500_000)               # numpy: shape and SCALE
print(f"P(10th ticket within 1.5 h): Poisson identity = {via_poisson:.4f}, "
      f"sum of exponentials = {np.mean(sum_of_exp <= t):.4f}, rng.gamma = {np.mean(direct <= t):.4f}")
print(f"mean = {direct.mean():.4f} h (alpha/beta = {n / lam:.4f}), variance = {direct.var():.4f} (alpha/beta^2 = {n / lam ** 2:.4f})")

# The shape controls skewness: 2 / sqrt(alpha)
def skewness(x):
    z = (x - x.mean()) / x.std()
    return np.mean(z ** 3)

print(" alpha   skewness (sample)   2/sqrt(alpha)")
for alpha in [0.5, 1, 2, 5, 20, 100]:
    x = rng.gamma(alpha, 1.0, size=500_000)
    print(f"{alpha:6}   {skewness(x):8.3f}            {2 / np.sqrt(alpha):.3f}")

P(10th ticket within 1.5 h): Poisson identity = 0.4126, sum of exponentials = 0.4128, rng.gamma = 0.4121
mean = 1.6676 h (alpha/beta = 1.6667), variance = 0.2789 (alpha/beta^2 = 0.2778)
 alpha   skewness (sample)   2/sqrt(alpha)
   0.5      2.815            2.828
     1      1.994            2.000
     2      1.414            1.414
     5      0.892            0.894
    20      0.444            0.447
   100      0.197            0.200

Key takeaways

• Gamma$(\alpha, \beta)$ has density proportional to $x^{\alpha - 1}e^{-\beta x}$, mean $\alpha/\beta$ and variance $\alpha/\beta^2$.
• With integer shape $n$, it is the waiting time for the $n$-th event of a Poisson process (sum of $n$ exponentials).
• Gammas with the same rate add; the exponential ($\alpha = 1$) and chi-square ($\beta = 1/2$) are special cases.
• Skewness $2/\sqrt{\alpha}$: large shapes look normal. Check whether software expects a rate or a scale.

🧪 Practice

1. $X \sim$ Gamma$(4, 2)$ (rate 2). Find the mean, the standard deviation and the mode.
2. Use the MGF $M(t) = (\beta/(\beta - t))^\alpha$ to prove additivity, and to derive the mean and variance.
3. Compute $\Gamma(5/2)$ from $\Gamma(1/2) = \sqrt{\pi}$ and the recurrence. Check with math.gamma.
4. Interview: A batch of 50 independent jobs runs sequentially, each taking an exponential time with mean 2 seconds. Estimate the 99th percentile of the batch duration. Hint: the total is Gamma$(50, 0.5)$; a normal approximation with the same mean and variance gives a quick answer, and the remaining right skew ($2/\sqrt{50}$) tells you whether the true quantile is above or below it. Check by simulation.

Beta

Why it matters. Many quantities are proportions: a conversion rate, a batting average, the fraction of a day a server is busy, the share of votes for a candidate. The beta distribution lives on $[0, 1]$ and is the standard model for an uncertain probability. It is the conjugate prior of the binomial, which makes Bayesian A/B testing a matter of adding counts (Section 11.1), and it describes the order statistics of uniform samples (Section 6.4).

Definition. $X \sim$ Beta$(\alpha, \beta)$, with $\alpha, \beta > 0$, if

$$f(x) = \frac{x^{\alpha - 1}(1 - x)^{\beta - 1}}{B(\alpha, \beta)}, \quad 0 < x < 1, \qquad B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}.$$

$$E[X] = \frac{\alpha}{\alpha + \beta}, \qquad \operatorname{Var}(X) = \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)}, \qquad \text{mode} = \frac{\alpha - 1}{\alpha + \beta - 2} ;; (\alpha, \beta > 1).$$

Shapes. The two parameters pull mass toward 1 ($\alpha$) and toward 0 ($\beta$):

$(\alpha, \beta)$	Shape	Mean
$(1, 1)$	flat: Uniform$(0, 1)$	0.5
$(2, 2)$	symmetric hump, $6x(1 - x)$ (Section 4.2)	0.5
$(0.5, 0.5)$	U-shaped: mass near 0 and 1	0.5
$(2, 5)$	hump near 0.2, right-skewed	0.286
$(5, 2)$	hump near 0.8, left-skewed	0.714
$(50, 50)$	narrow hump at 0.5	0.5

Pseudo-counts. The most useful way to read the parameters: Beta$(\alpha, \beta)$ is roughly what you know about a success probability after seeing $\alpha - 1$ successes and $\beta - 1$ failures, starting from a uniform prior. The mean is the success fraction, and the sum $\alpha + \beta$ measures how concentrated the distribution is: the variance shrinks like $1/(\alpha + \beta)$.

Bayesian updating. If the prior for $p$ is Beta$(\alpha, \beta)$ and you observe $k$ successes in $n$ Bernoulli$(p)$ trials, the posterior is

$$p \mid \text{data} \sim \text{Beta}(\alpha + k, ; \beta + n - k).$$

The successes are added to $\alpha$ and the failures to $\beta$. Section 11.1 derives this; here it is a convenient fact.

Other connections.

• If $X \sim$ Gamma$(\alpha, 1)$ and $Y \sim$ Gamma$(\beta, 1)$ are independent, $X/(X + Y) \sim$ Beta$(\alpha, \beta)$: the share of one gamma in a total.
• The $k$-th smallest of $n$ independent Uniform$(0, 1)$ values is Beta$(k, n - k + 1)$, with mean $k/(n + 1)$.

import numpy as np

rng = np.random.default_rng(12)

# Bayesian A/B test with a uniform Beta(1, 1) prior on each conversion rate
data = {"A": (120, 1000), "B": (140, 1000)}            # (conversions, visitors)
post = {}
for name, (k, n) in data.items():
    post[name] = rng.beta(1 + k, 1 + n - k, size=1_000_000)   # draws from Beta(1 + k, 1 + n - k)
    lo, hi = np.quantile(post[name], [0.025, 0.975])
    print(f"{name}: posterior Beta({1 + k}, {1 + n - k}), mean = {post[name].mean():.4f}, "
          f"95% interval = ({lo:.4f}, {hi:.4f})")
print(f"P(rate B > rate A) = {np.mean(post['B'] > post['A']):.3f}")
print(f"expected relative lift of B over A = {np.mean(post['B'] / post['A'] - 1):+.1%}")

# Order statistics: the 3rd smallest of 10 uniforms is Beta(3, 8)
u3 = np.sort(rng.random((200_000, 10)), axis=1)[:, 2]
a, b = 3, 8
print(f"3rd smallest of 10 uniforms: mean = {u3.mean():.4f} (3/11 = {3 / 11:.4f}), "
      f"variance = {u3.var():.5f} (Beta(3, 8): {a * b / ((a + b) ** 2 * (a + b + 1)):.5f})")

A: posterior Beta(121, 881), mean = 0.1208, 95% interval = (0.1013, 0.1417)
B: posterior Beta(141, 861), mean = 0.1407, 95% interval = (0.1199, 0.1630)
P(rate B > rate A) = 0.908
expected relative lift of B over A = +17.4%
3rd smallest of 10 uniforms: mean = 0.2724 (3/11 = 0.2727), variance = 0.01658 (Beta(3, 8): 0.01653)

The posterior answers the question a product manager actually asks ("how likely is B better than A, and by how much?") directly, by simulation from two beta distributions.

Key takeaways

• Beta$(\alpha, \beta)$ is a flexible distribution on $[0, 1]$ for proportions and probabilities.
• Mean $\alpha/(\alpha + \beta)$; the sum $\alpha + \beta$ controls the concentration.
• Read the parameters as pseudo-counts of successes and failures; observing data adds the counts (Beta-binomial conjugacy).
• Beta$(1, 1)$ is uniform, and order statistics of uniforms are beta.

🧪 Practice

1. Compute the mean, variance and mode of Beta$(2, 5)$.
2. Choose $\alpha$ and $\beta$ for a prior on a click-through rate with mean 0.02 that is "worth" 500 observations ($\alpha + \beta = 500$). After 10 clicks in 300 impressions, what is the posterior mean?
3. Verify by simulation that $X/(X + Y)$ is Beta$(2, 3)$ when $X \sim$ Gamma$(2, 1)$ and $Y \sim$ Gamma$(3, 1)$ are independent, by comparing means and variances.
4. Interview: You need to rank products by customer rating, where some products have 2 reviews (both positive) and others have 2,000 reviews (95% positive). How do you rank them? Hint: the raw proportion puts the 2-review product first; use a beta posterior, then rank by the posterior mean under a sensible prior or by a lower credible bound, and explain how the prior acts as a set of pseudo-reviews.

Weibull

Why it matters. The exponential distribution assumes that items never age, which is false for most physical things. The Weibull distribution adds one shape parameter that lets the failure rate fall, stay constant or rise over time. It is the workhorse of reliability engineering (bearing lives, battery cycles, fatigue failures), is used for wind speeds and material strengths, and arises naturally as the distribution of the weakest link in a chain of many components.

Definition. $X \sim$ Weibull$(k, \lambda)$, with shape $k > 0$ and scale $\lambda > 0$, if

$$F(x) = 1 - e^{-(x/\lambda)^k}, \qquad f(x) = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k - 1} e^{-(x/\lambda)^k}, \qquad x \ge 0.$$

For $k = 1$ it is Exponential$(1/\lambda)$. In general, $X = \lambda E^{1/k}$ with $E \sim$ Exponential$(1)$.

The hazard rate tells the story.

$$h(x) = \frac{f(x)}{1 - F(x)} = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k - 1}.$$

Shape	Hazard	Failure mode	Examples
$k < 1$	decreasing	infant mortality: weak units fail early, survivors are robust	manufacturing defects, early software crashes
$k = 1$	constant	random failures (exponential)	external shocks
$k > 1$	increasing	wear-out: older units fail more often	bearings, tires, batteries ($k \approx 2$ to 4)

Over a product's life, the three phases often appear in sequence, giving the bathtub curve:

  hazard
    |\
    | \                                              /
    |  \                                           /
    |   '-.__________________________________..--'
    |   infant         useful life (random         wear-out
    |   mortality      failures, k ~ 1)            (k > 1)
    +-----------------------------------------------------> age

Summaries.

Quantity	Value
Mean	$\lambda,\Gamma(1 + 1/k)$
Variance	$\lambda^2\left[\Gamma(1 + 2/k) - \Gamma(1 + 1/k)^2\right]$
Quantile	$Q(p) = \lambda,(-\ln(1 - p))^{1/k}$
Characteristic life	$F(\lambda) = 1 - e^{-1} \approx 63.2%$ for every $k$
B10 life	$Q(0.10)$: the age by which 10% of units have failed, a standard warranty metric

Weibull plot. Taking logarithms twice linearizes the CDF:

$$\ln\left(-\ln(1 - F(x))\right) = k \ln x - k \ln \lambda.$$

Plotting $\ln(-\ln(1 - \hat F))$ against $\ln x$ for observed failure times gives a straight line if the data are Weibull, with slope $k$. The slope immediately tells the engineer whether failures are infant mortality, random, or wear-out.

import numpy as np
from math import exp, gamma, log

rng = np.random.default_rng(13)

# Bearings: Weibull with shape k = 2.5 (wear-out) and scale 1000 hours
k, lam = 2.5, 1000.0
S = lambda x: exp(-(x / lam) ** k)                       # survival function
print(f"mean = {lam * gamma(1 + 1 / k):.1f} h, B10 life = {lam * (-log(0.9)) ** (1 / k):.1f} h, "
      f"P(survive 500 h) = {S(500):.3f}")
print(f"P(survive 500 more h | survived 500 h) = {S(1000) / S(500):.3f}   (an aging unit is worse than a new one)")

# Hazard rates per 1,000 hours at several ages, for three shapes with the same scale
print("  age    k = 0.5   k = 1   k = 2.5")
for x in [100, 500, 1000, 2000]:
    h = [kk / lam * (x / lam) ** (kk - 1) * 1000 for kk in (0.5, 1, 2.5)]
    print(f"{x:5d}    {h[0]:6.2f}   {h[1]:5.2f}   {h[2]:6.2f}")

# Weibull plot: estimate the shape from 40 observed failure times
t = np.sort(lam * rng.weibull(k, size=40))               # numpy's weibull has scale 1
i = np.arange(1, len(t) + 1)
F_hat = (i - 0.3) / (len(t) + 0.4)                       # median-rank plotting positions
slope, intercept = np.polyfit(np.log(t), np.log(-np.log(1 - F_hat)), 1)
print(f"Weibull plot fit: shape = {slope:.2f} (true {k}), scale = {np.exp(-intercept / slope):.0f} (true {lam:.0f})")

mean = 887.3 h, B10 life = 406.5 h, P(survive 500 h) = 0.838
P(survive 500 more h | survived 500 h) = 0.439   (an aging unit is worse than a new one)
  age    k = 0.5   k = 1   k = 2.5
  100      1.58    1.00     0.08
  500      0.71    1.00     0.88
 1000      0.50    1.00     2.50
 2000      0.35    1.00     7.07
Weibull plot fit: shape = 2.23 (true 2.5), scale = 947 (true 1000)

With only 40 failures the estimated shape is rough, but it clearly signals wear-out ($k > 1$), which is the decision the plot is usually made for.

Key takeaways

• Weibull$(k, \lambda)$ has survival $e^{-(x/\lambda)^k}$ and hazard proportional to $x^{k-1}$.
• $k < 1$: decreasing hazard (infant mortality); $k = 1$: exponential; $k > 1$: wear-out.
• 63.2% of units fail by the scale $\lambda$, whatever the shape; B10 life is the 10% quantile.
• A Weibull plot of $\ln(-\ln(1 - \hat F))$ against $\ln x$ is a straight line with slope $k$.

🧪 Practice

1. Show that the Weibull with $k = 1$ is the exponential distribution, and that $\lambda E^{1/k}$ is Weibull when $E \sim$ Exponential$(1)$.
2. For $k = 0.5$ and $\lambda = 100$, compute the hazard at ages 1, 10 and 100. Would a burn-in period (running every unit for a while before shipping) help? Would it help for $k = 3$?
3. Derive the B10 life formula from the CDF, and compute it for $k = 1.5$, $\lambda = 5$ years.
4. Interview: A disk vendor quotes an MTBF (mean time between failures) of 1,000,000 hours. Does that mean a typical disk lasts over 100 years? Hint: MTBF is derived from the failure rate of a large fleet during its useful life, where failures are roughly exponential; convert it to an annual failure rate for a fleet, and contrast it with wear-out after a few years, which a Weibull with $k > 1$ would capture.

Pareto and Heavy Tails

Why it matters. Wealth, city populations, file sizes, website traffic, insurance losses and the sizes of software projects share a pattern: most values are modest, and a few are enormous. A small fraction of items accounts for a large fraction of the total (the "80/20 rule"). The Pareto distribution is the basic model for this pattern. Its heavy tail breaks the intuitions built on the normal distribution: averages converge slowly or not at all, the largest observation can dominate the sum, and variances can be infinite.

Definition. $X \sim$ Pareto$(x_m, \alpha)$, with minimum $x_m > 0$ and tail index $\alpha > 0$, if

$$P(X > x) = \left(\frac{x_m}{x}\right)^{\alpha}, \qquad f(x) = \frac{\alpha,x_m^\alpha}{x^{\alpha + 1}}, \qquad x \ge x_m.$$

The survival function decays as a power of $x$, not exponentially. On a log-log plot it is a straight line with slope $-\alpha$:

  log P(X > x)                               log P(X > x)
    |'.                                        |'-._
    |   '.     Pareto: a straight line         |    '-.     exponential: bends down,
    |     '.   with slope -alpha               |       '.   steeper and steeper
    |       '.                                 |         \
    |         '.                               |          |
    +---------------------> log x              +---------------------> log x

Moments exist only below the tail index. $E[X^k] < \infty$ only for $k < \alpha$:

Tail index	Mean	Variance	Consequence
$\alpha \le 1$	infinite	infinite	the sample mean grows without bound
$1 < \alpha \le 2$	$\frac{\alpha x_m}{\alpha - 1}$	infinite	the sample mean converges, but very slowly and erratically
$\alpha > 2$	finite	$\frac{\alpha x_m^2}{(\alpha - 1)^2(\alpha - 2)}$	the usual statistics work, but extremes still matter

Many real datasets have $1 < \alpha < 2$ (city sizes, file sizes, wealth), so their variance is effectively infinite.

Scale invariance. For every $x \ge x_m$ and $c > 1$, $P(X > cx \mid X > x) = c^{-\alpha}$. Among cities above 100,000 people, the fraction above 200,000 is the same as the fraction above 2 million among cities above 1 million. Conditioning on exceeding a threshold gives another Pareto, and the mean excess grows with the threshold: $E[X \mid X > t] = \frac{\alpha}{\alpha - 1}t$. For an exponential, by contrast, the mean excess is constant (memorylessness). "The bigger it already is, the bigger it is likely to get" is the signature of a heavy tail.

The 80/20 rule. The top fraction $q$ of a Pareto population holds a share $q^{1 - 1/\alpha}$ of the total. The share of the top 20% is 80% when $\alpha = \log_4 5 \approx 1.16$.

Heavy tails in general. A distribution is heavy-tailed if its tail decays more slowly than any exponential.

Tail class	Examples	Behavior of sample averages
Bounded	uniform, beta	converge fast
Light (exponential or faster)	normal, exponential, gamma, Poisson	converge at rate $1/\sqrt{n}$; CLT works well
Heavy, all moments finite	lognormal, Weibull with $k < 1$	converge, but slowly; CLT needs large $n$
Power law	Pareto, Student's $t$, Cauchy	may not converge at all; one value can dominate

import numpy as np

rng = np.random.default_rng(14)

def pareto(alpha, xm, size):
    # numpy's rng.pareto samples the Lomax (Pareto II) distribution; add 1 and scale to get Pareto(xm, alpha)
    return xm * (1 + rng.pareto(alpha, size))

# Running means: finite variance (alpha = 3), infinite variance (1.5), infinite mean (0.9)
print(" alpha   theory    mean of first 10^3   10^4      10^5      10^6    largest / total")
for alpha in [3.0, 1.5, 0.9]:
    x = pareto(alpha, 1.0, 1_000_000)
    means = [x[:n].mean() for n in (10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6)]
    theory = alpha / (alpha - 1) if alpha > 1 else float("inf")
    print(f"{alpha:5.1f}   {theory:6.2f}   " + "   ".join(f"{m:8.2f}" for m in means) + f"   {x.max() / x.sum():8.3f}")

# 80/20: share of the total held by the top 20% when alpha = 1.16, in five samples of 10^6
shares = []
for _ in range(5):
    x = np.sort(pareto(1.16, 1.0, 1_000_000))[::-1]
    shares.append(x[:200_000].sum() / x.sum())
print("top 20% share in five samples:", " ".join(f"{s:.3f}" for s in shares),
      f"  theory 0.2^(1 - 1/alpha) = {0.2 ** (1 - 1 / 1.16):.3f}")

# Estimating alpha from the log-log slope of the empirical survival function
x = np.sort(pareto(1.5, 1.0, 100_000))
surv = 1 - np.arange(1, len(x) + 1) / (len(x) + 1)     # empirical P(X > x)
keep = (x > 2) & (x < 200)                             # avoid the noisy extreme tail
slope = np.polyfit(np.log(x[keep]), np.log(surv[keep]), 1)[0]
print(f"log-log slope = {slope:.3f}  ->  alpha estimate = {-slope:.3f} (true 1.5)")

 alpha   theory    mean of first 10^3   10^4      10^5      10^6    largest / total
  3.0     1.50       1.52       1.49       1.50       1.50      0.000
  1.5     3.00       2.90       3.23       2.95       2.96      0.001
  0.9      inf      20.78      22.80      40.18      49.28      0.218
top 20% share in five samples: 0.980 0.791 0.783 0.783 0.761   theory 0.2^(1 - 1/alpha) = 0.801
log-log slope = -1.503  ->  alpha estimate = 1.503 (true 1.5)

With $\alpha = 3$ the running mean settles quickly. With $\alpha = 1.5$ it converges, but erratically. With $\alpha = 0.9$ it keeps growing, and a single observation can be a large part of the entire total. Even the top-20% share, a summary of a million observations, changes noticeably from one sample to the next, because the denominator is a sum with infinite variance.

Key takeaways

• Pareto$(x_m, \alpha)$: $P(X > x) = (x_m/x)^\alpha$, a straight line on log-log axes.
• Moments of order $k \ge \alpha$ are infinite; $1 < \alpha \le 2$ gives a finite mean but infinite variance.
• Power laws are scale-invariant: the mean excess grows with the threshold, and a few items hold most of the total.
• With heavy tails, report medians and quantiles, look at the largest values, and distrust sample means and standard deviations.

🧪 Practice

1. For Pareto$(1, 3)$, compute the mean, the variance and the median.
2. Show that if $X \sim$ Pareto$(x_m, \alpha)$ and $t \ge x_m$, then $X$ conditioned on $X > t$ is Pareto$(t, \alpha)$. Use this to derive $E[X \mid X > t]$.
3. Derive an inverse transform sampler for the Pareto (Section 4.4) and use it instead of rng.pareto in the code. Then estimate $\alpha$ with the Hill estimator $\hat\alpha = m / \sum_{i=1}^{m} \ln(x_{(i)}/x_{(m+1)})$, where $x_{(1)} \ge x_{(2)} \ge \dots$ are the data sorted from largest to smallest and $m$ is, say, 1,000.
4. Interview: Average revenue per user jumped 30% this week and the team wants to celebrate. What do you check first? Hint: revenue per user is heavy-tailed; look at the median, a trimmed mean, the share of revenue from the top few users, and whether one large customer (a "whale") explains the change.

Lognormal

Why it matters. When a quantity is the product of many independent positive factors, its logarithm is a sum, and sums tend to be normal. The quantity itself is then lognormal. This is the mechanism behind incomes (raises are percentages), stock prices (returns compound), file sizes, response times, particle sizes and the durations of tasks. The lognormal is right-skewed and positive, and it is often confused with the Pareto because both have long right tails. Section 4.4 derived its density.

Definition. $X$ is lognormal, $X \sim$ Lognormal$(\mu, \sigma^2)$, if $\ln X \sim$ Normal$(\mu, \sigma^2)$, that is, $X = e^Y$ with $Y$ normal:

$$f(x) = \frac{1}{x,\sigma\sqrt{2\pi}} \exp\left(-\frac{(\ln x - \mu)^2}{2\sigma^2}\right), \qquad x > 0.$$

The parameters $\mu$ and $\sigma$ are the mean and SD of $\ln X$, not of $X$, a frequent source of errors in software (numpy's rng.lognormal(mean, sigma) takes the parameters of the log).

Summaries. Quantiles pass through the increasing function $e^y$ (Section 4.2), but means do not:

Quantity	Value	Note
Median	$e^{\mu}$	the geometric mean of $X$
Mode	$e^{\mu - \sigma^2}$	below the median
Mean	$e^{\mu + \sigma^2/2}$	above the median: mode < median < mean
Variance	$(e^{\sigma^2} - 1),e^{2\mu + \sigma^2}$
Coefficient of variation	$\sqrt{e^{\sigma^2} - 1}$	depends only on $\sigma$
Quantile	$e^{\mu + \sigma z_p}$	$z_p$ = standard normal quantile (Section 5.3)

Multiplicative thinking. On the log scale, the familiar normal intervals become "times or divide" intervals. About 68% of values lie between $e^\mu / e^\sigma$ and $e^\mu \cdot e^\sigma$, and $e^\sigma$ is called the geometric standard deviation. A latency distribution with median 100 ms and geometric SD 1.8 has 68% of requests between $100/1.8 \approx 56$ ms and $100 \times 1.8 = 180$ ms.

Working on the log scale. To analyze lognormal data, take logs, use normal-theory methods on $\ln X$, and transform back. Be careful with the last step: $\exp(\text{mean of the logs})$ estimates the median of $X$, not its mean. To estimate the mean, multiply by $e^{\hat\sigma^2/2}$.

Lognormal vs. Pareto. All moments of the lognormal are finite, but it has no MGF, and for large $\sigma$ its tail looks like a straight line on a log-log plot over several orders of magnitude. Distinguishing the two from data requires a lot of data in the tail, and the choice matters for extreme quantiles.

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(15)

def skewness(x):
    z = (x - x.mean()) / x.std()
    return np.mean(z ** 3)

# 1. Multiplicative growth: 100,000 firms, 50 years of random yearly growth factors in [0.8, 1.25]
size = rng.uniform(0.8, 1.25, size=(100_000, 50)).prod(axis=1)
print(f"firm size: skewness = {skewness(size):.2f};  log(size): skewness = {skewness(np.log(size)):+.3f}")

# 2. Latencies: lognormal with median 100 ms and sigma = 0.6 (geometric SD e^0.6 = 1.82)
mu, sigma = np.log(100), 0.6
z99 = NormalDist().inv_cdf(0.99)
print(f"theory: mode = {np.exp(mu - sigma ** 2):.1f}, median = {np.exp(mu):.1f}, "
      f"mean = {np.exp(mu + sigma ** 2 / 2):.1f}, p99 = {np.exp(mu + sigma * z99):.1f} ms")
lat = rng.lognormal(mean=mu, sigma=sigma, size=1_000_000)      # numpy takes the log-scale parameters
print(f"sample: median = {np.median(lat):.1f}, mean = {lat.mean():.1f}, p99 = {np.quantile(lat, 0.99):.1f} ms, "
      f"within [median / 1.82, median * 1.82]: {np.mean((lat > 100 / np.exp(sigma)) & (lat < 100 * np.exp(sigma))):.3f}")

# 3. Back-transforming: exp(mean of logs) estimates the median, not the mean
sample = lat[:200]
logs = np.log(sample)
naive = np.exp(logs.mean())
corrected = naive * np.exp(logs.var(ddof=1) / 2)
print(f"200 latencies: sample mean = {sample.mean():.1f}, exp(mean log) = {naive:.1f}, "
      f"with e^(s^2/2) correction = {corrected:.1f}  (true mean {np.exp(mu + sigma ** 2 / 2):.1f})")

firm size: skewness = 4.36;  log(size): skewness = -0.019
theory: mode = 69.8, median = 100.0, mean = 119.7, p99 = 403.8 ms
sample: median = 99.9, mean = 119.5, p99 = 400.9 ms, within [median / 1.82, median * 1.82]: 0.683
200 latencies: sample mean = 120.4, exp(mean log) = 99.4, with e^(s^2/2) correction = 120.7  (true mean 119.7)

Key takeaways

• $X$ is lognormal when $\ln X$ is normal; products of many positive factors are approximately lognormal.
• $\mu, \sigma$ describe $\ln X$. Median $e^\mu$, mean $e^{\mu + \sigma^2/2}$, and mode < median < mean.
• Think multiplicatively: geometric mean and geometric SD give "times or divide" intervals.
• Analyze on the log scale, but remember that back-transformed means estimate medians unless corrected.

🧪 Practice

1. Household incomes are lognormal with median 50,000 and mean 65,000. Find $\mu$ and $\sigma$, and the fraction of households earning more than 100,000.
2. Show that the product of two independent lognormal variables is lognormal, and find its parameters. What happens to $1/X$?
3. Simulate a quantity that grows by a random percentage each year and check that its logarithm, not the quantity itself, has a symmetric distribution. Repeat with additive random increments and compare.
4. Interview: You fit a linear regression to log(revenue) and forecast revenue as exp(prediction). Summed over all stores, your forecasts are consistently below actual revenue. Why? Hint: exp of the expected log is the conditional median, not the conditional mean; think about what factor $e^{\sigma^2/2}$ does and how a "smearing" correction estimates it from residuals.

<a id="53-the-normal-distribution"></a>

5.3 The Normal Distribution

The normal distribution is the most important distribution in statistics, because sums and averages of many small independent effects are approximately normal (the central limit theorem, Section 6.5). This subchapter covers its properties, how to compute normal probabilities, when it approximates discrete distributions, and how to check whether data are normal.

Normal Density and Properties

Why it matters. Measurement errors, heights, test scores and the average of many observations all tend to follow the same bell-shaped curve. The reason is the central limit theorem: whenever a quantity is the sum of many small, independent contributions, its distribution is close to normal, whatever the contributions look like. Add mathematical convenience (the normal is fully described by its mean and variance, and linear combinations of normals stay normal), and you get the default model behind confidence intervals, $t$-tests, regression and control charts.

Definition. $X \sim$ Normal$(\mu, \sigma^2)$, also written $N(\mu, \sigma^2)$, if

$$f(x) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right), \qquad -\infty < x < \infty.$$

$\mu$ is the mean (and the median and mode), and $\sigma$ is the standard deviation. Note the convention: the second parameter is the variance $\sigma^2$ in mathematical notation, but most software (numpy's rng.normal(loc, scale), Python's NormalDist(mu, sigma)) takes the standard deviation.

Shape.

                         _.-'''-._
                      .-'    |    '-.            inflection points at mu - sigma and mu + sigma
                    .'       |       '.
                  .'|        |        |'.
               _.'  |        |        |  '._
         __..-'     |        |        |     '-..__
  ---''''-----+-----+--------+--------+-----+-----''''---
           mu-2s  mu-s      mu      mu+s  mu+2s
                    |<---- 68.3% ---->|
              |<---------- 95.4% ---------->|

Property	Value
Symmetry	$f(\mu + t) = f(\mu - t)$; mean = median = mode = $\mu$
Inflection points	$\mu \pm \sigma$
Skewness, kurtosis	0 and 3 (excess kurtosis 0, the reference value of Section 2.4)
Within $1, 2, 3\sigma$	68.27%, 95.45%, 99.73% (the empirical rule of Section 2.3)
Tails	decay like $e^{-x^2/2}$: faster than exponential, so extreme values are extremely rare

Why the constant is $\sqrt{2\pi}$. The integral $I = \int_{-\infty}^{\infty} e^{-x^2/2},dx$ has no elementary antiderivative, but its square can be computed in polar coordinates: $I^2 = \iint e^{-(x^2 + y^2)/2},dx,dy = \int_0^{2\pi}\int_0^\infty e^{-r^2/2},r,dr,d\theta = 2\pi$. So $I = \sqrt{2\pi}$.

Closure properties. These make the normal uniquely convenient.

Operation	Result
Linear transformation	$aX + b \sim N(a\mu + b, ; a^2\sigma^2)$
Standardization	$Z = (X - \mu)/\sigma \sim N(0, 1)$
Sum of independent normals	$X_1 + X_2 \sim N(\mu_1 + \mu_2, ; \sigma_1^2 + \sigma_2^2)$
Difference of independent normals	$X_1 - X_2 \sim N(\mu_1 - \mu_2, ; \sigma_1^2 + \sigma_2^2)$ (variances still add)
Mean of $n$ independent $N(\mu, \sigma^2)$	$\bar X \sim N(\mu, \sigma^2/n)$ exactly, for every $n$

Among all distributions with a given mean and variance, the normal has the largest entropy (Section 12.5): it is the "least informative" choice when only those two quantities are known.

import numpy as np
from statistics import NormalDist

# 68-95-99.7, and how quickly the tails vanish
Z = NormalDist(0, 1)
print(" k   P(|Z| <= k)    P(|Z| > k)    about 1 in")
for k in range(1, 7):
    outside = 2 * (1 - Z.cdf(k))
    print(f" {k}   {1 - outside:.7f}      {outside:.2e}     {1 / outside:,.0f}")

# Differences of independent normals: men's heights N(178, 7^2), women's heights N(165, 6^2)
M, W = NormalDist(178, 7), NormalDist(165, 6)
D = M - W                                         # NormalDist arithmetic assumes independence
print(f"D = M - W: mean = {D.mean}, SD = {D.stdev:.3f} (sqrt(49 + 36));  P(random woman taller) = {D.cdf(0):.4f}")

rng = np.random.default_rng(17)
m, w = rng.normal(178, 7, 1_000_000), rng.normal(165, 6, 1_000_000)
print(f"simulation: SD of M - W = {np.std(m - w):.3f},  P(W > M) = {np.mean(w > m):.4f}")

# Linear transformation: temperatures N(20, 5^2) in Celsius, converted to Fahrenheit
F = NormalDist(20, 5) * 9 / 5 + 32
print(f"Fahrenheit: N({F.mean:.0f}, {F.stdev:.0f}^2)")

 k   P(|Z| <= k)    P(|Z| > k)    about 1 in
 1   0.6826895      3.17e-01     3
 2   0.9544997      4.55e-02     22
 3   0.9973002      2.70e-03     370
 4   0.9999367      6.33e-05     15,787
 5   0.9999994      5.73e-07     1,744,278
 6   1.0000000      1.97e-09     506,797,317
D = M - W: mean = 13.0, SD = 9.220 (sqrt(49 + 36));  P(random woman taller) = 0.0793
simulation: SD of M - W = 9.225,  P(W > M) = 0.0789
Fahrenheit: N(68, 9^2)

Key takeaways

• The normal density is a symmetric bell with mean $\mu$ and SD $\sigma$; 68%, 95% and 99.7% of the mass lie within 1, 2 and 3 SDs.
• Linear transformations and sums of independent normals are normal; means add and variances add (also for differences).
• The tails are extremely thin: a 6-SD event has probability about $2 \times 10^{-9}$ under normality.
• Its importance comes from the central limit theorem and from these closure properties.

🧪 Practice

1. Show by differentiating the density twice that the inflection points are at $\mu \pm \sigma$.
2. A machine fills bags with $N(500, 4^2)$ grams, and the empty bag weighs $N(10, 1^2)$ grams, independently. Find the distribution of the total weight and the probability that it exceeds 520 grams.
3. $X_1, \dots, X_{25}$ are independent $N(100, 15^2)$. Find the distribution of $\bar X$ and $P(\bar X > 105)$. Compare with $P(X_1 > 105)$.
4. Interview: "Six Sigma" quality programs claim that a process with specification limits 6 SDs from the mean produces 3.4 defects per million, but the normal tail beyond 6 SDs is about 1 per billion. Where does 3.4 come from? Hint: the method assumes the process mean drifts by 1.5 SDs in the long run; compute the one-sided tail beyond 4.5 SDs, and discuss whether the normal tail is trustworthy that far out.

Standard Normal and z-Tables

Why it matters. There are infinitely many normal distributions, but they are all rescaled copies of one: the standard normal $Z \sim N(0, 1)$. Its CDF $\Phi$ has no formula in elementary functions, so for a century people used printed tables of $\Phi$; today software computes it. Either way, any normal probability is computed the same way: standardize, then look up $\Phi$. The same few critical values (1.645, 1.96, 2.576) appear in every confidence interval and test in Chapters 8 and 9.

Definition. The standard normal has density and CDF

$$\varphi(z) = \frac{1}{\sqrt{2\pi}}e^{-z^2/2}, \qquad \Phi(z) = P(Z \le z) = \int_{-\infty}^{z}\varphi(t),dt = \frac{1}{2}\left[1 + \operatorname{erf}\left(\frac{z}{\sqrt 2}\right)\right].$$

Two directions. If $X \sim N(\mu, \sigma^2)$:

Question	Method	Formula
Value $\to$ probability: "what fraction scores below $x$?"	standardize, then $\Phi$	$P(X \le x) = \Phi\left(\frac{x - \mu}{\sigma}\right)$
Probability $\to$ value: "which score is the 90th percentile?"	$\Phi^{-1}$, then un-standardize	$x_p = \mu + \sigma,\Phi^{-1}(p)$

Symmetry rules. Tables list only $\Phi(z)$ for $z \ge 0$; the rest follows from symmetry:

$$\Phi(-z) = 1 - \Phi(z), \qquad P(Z > z) = 1 - \Phi(z), \qquad P(|Z| \le z) = 2\Phi(z) - 1.$$

Reading a z-table. Rows give $z$ to one decimal, columns give the second decimal. To find $\Phi(1.96)$, go to row 1.9 and column 0.06, giving 0.9750. To find the 95th percentile, search the body of the table for 0.95: it lies between 0.9495 (1.64) and 0.9505 (1.65), so $z \approx 1.645$.

Critical values to know.

Central probability	Two-sided $z$	One-sided $z$ (same tail area)
80%	1.282	0.842
90%	1.645	1.282
95%	1.960	1.645
99%	2.576	2.326
99.9%	3.291	3.090

from math import erf, sqrt
from statistics import NormalDist

Phi = lambda z: 0.5 * (1 + erf(z / sqrt(2)))            # standard normal CDF via the error function

# A slice of the classic z-table: P(Z <= z) for z = row + column
cols = [i / 100 for i in range(10)]
print("  z  " + "".join(f"  {c:.2f} " for c in cols))
for row in [0.0, 0.5, 1.0, 1.2, 1.6, 1.9, 2.3, 2.5, 3.0]:
    print(f" {row:.1f} " + "".join(f" {Phi(row + c):.4f}" for c in cols))

# Forward and inverse problems: exam scores ~ N(1050, 200^2)
X = NormalDist(1050, 200)
z = (1300 - 1050) / 200
print(f"P(score > 1300): z = {z:.2f}, 1 - Phi(z) = {1 - Phi(z):.4f}  (NormalDist: {1 - X.cdf(1300):.4f})")
print(f"P(900 < score < 1200) = {X.cdf(1200) - X.cdf(900):.4f}")
print(f"90th percentile = 1050 + 200 * {NormalDist().inv_cdf(0.90):.4f} = {X.inv_cdf(0.90):.1f}")

  z    0.00   0.01   0.02   0.03   0.04   0.05   0.06   0.07   0.08   0.09
 0.0  0.5000 0.5040 0.5080 0.5120 0.5160 0.5199 0.5239 0.5279 0.5319 0.5359
 0.5  0.6915 0.6950 0.6985 0.7019 0.7054 0.7088 0.7123 0.7157 0.7190 0.7224
 1.0  0.8413 0.8438 0.8461 0.8485 0.8508 0.8531 0.8554 0.8577 0.8599 0.8621
 1.2  0.8849 0.8869 0.8888 0.8907 0.8925 0.8944 0.8962 0.8980 0.8997 0.9015
 1.6  0.9452 0.9463 0.9474 0.9484 0.9495 0.9505 0.9515 0.9525 0.9535 0.9545
 1.9  0.9713 0.9719 0.9726 0.9732 0.9738 0.9744 0.9750 0.9756 0.9761 0.9767
 2.3  0.9893 0.9896 0.9898 0.9901 0.9904 0.9906 0.9909 0.9911 0.9913 0.9916
 2.5  0.9938 0.9940 0.9941 0.9943 0.9945 0.9946 0.9948 0.9949 0.9951 0.9952
 3.0  0.9987 0.9987 0.9987 0.9988 0.9988 0.9989 0.9989 0.9989 0.9990 0.9990
P(score > 1300): z = 1.25, 1 - Phi(z) = 0.1056  (NormalDist: 0.1056)
P(900 < score < 1200) = 0.5467
90th percentile = 1050 + 200 * 1.2816 = 1306.3

Key takeaways

• Every normal probability reduces to the standard normal CDF $\Phi$ through $z = (x - \mu)/\sigma$.
• Percentiles go the other way: $x_p = \mu + \sigma,\Phi^{-1}(p)$.
• Symmetry, $\Phi(-z) = 1 - \Phi(z)$, handles negative $z$ and two-sided probabilities.
• Memorize 1.645, 1.96 and 2.576; in code use statistics.NormalDist, math.erf or scipy.stats.norm.

🧪 Practice

1. Use the table to find $P(-1.2 < Z < 0.8)$, $P(Z > 2.35)$ and the value $z$ with $P(|Z| > z) = 0.01$.
2. A bottling machine fills $N(\mu, 4^2)$ ml. Where should $\mu$ be set so that only 1% of bottles contain less than the 495 ml printed on the label?
3. IQ scores are $N(100, 15^2)$. What fraction of people score above 130? What score marks the top 2%?
4. Interview: Implement $\Phi(z)$ without any statistics library, accurate to 6 decimal places. Hint: options include numerical integration of $\varphi$ from 0 to $z$ (Simpson's rule) plus symmetry, the Taylor series of $e^{-t^2/2}$ integrated term by term for small $|z|$, and published rational approximations for the tails; test against math.erf and watch the accuracy for large $|z|$.

Normal Approximation to the Binomial

Why it matters. A Binomial$(n, p)$ count is a sum of $n$ independent Bernoulli variables, so by the central limit theorem it is approximately normal when $n$ is large. Historically this result (de Moivre, 1733; Laplace, 1812) was the first central limit theorem, and it made binomial probabilities with large $n$ computable by hand. Today exact binomial probabilities are cheap, but the approximation is still how people reason about counts: a proportion from a sample of $n$ is roughly normal with standard error $\sqrt{p(1-p)/n}$, which is the basis of the margin of error in polls and of the $z$-test for proportions (Section 9.2).

Statement (de Moivre-Laplace). If $X \sim$ Binomial$(n, p)$ with $0 < p < 1$ fixed, then as $n \to \infty$

$$\frac{X - np}{\sqrt{np(1 - p)}} ;\longrightarrow; N(0, 1) \text{ in distribution}, \qquad \text{so} \qquad X \approx N\big(np, ; np(1 - p)\big).$$

Equivalently, the sample proportion $\hat p = X/n$ is approximately $N\big(p, ; p(1 - p)/n\big)$.

When is it good? The binomial is skewed when $p$ is close to 0 or 1, and the normal is symmetric, so the approximation needs enough expected successes and failures. The common rule of thumb is

$$np \ge 10 \quad \text{and} \quad n(1 - p) \ge 10$$

(some textbooks use 5). The quantity that really matters is the variance $np(1 - p)$: the larger it is, the better the approximation.

Situation	Best approximation
$n$ large, $p$ small, $np$ moderate (say $\le 10$)	Poisson$(np)$ (Section 5.1)
$np$ and $n(1-p)$ both $\ge 10$	Normal$(np, np(1-p))$, with continuity correction (next topic)
Small $n$, or exact tail probabilities needed	exact binomial (always available in software)

The same reasoning applies to the Poisson: Poisson$(\lambda) \approx N(\lambda, \lambda)$ for large $\lambda$ (say $\lambda \ge 20$), since a Poisson$(\lambda)$ count is a sum of $\lambda$ independent Poisson$(1)$ counts when $\lambda$ is an integer.

from math import comb, floor, sqrt
from statistics import NormalDist

def binom_cdf(k, n, p):
    return sum(comb(n, j) * p ** j * (1 - p) ** (n - j) for j in range(k + 1))

# P(X <= k) at one SD below the mean: exact vs normal approximation (no continuity correction yet)
print("     n      p   np(1-p)    k    exact    normal   abs. error")
for n, p in [(20, 0.5), (100, 0.5), (1000, 0.5), (20, 0.1), (100, 0.1), (100, 0.02), (1000, 0.02), (5000, 0.02)]:
    mu, sd = n * p, sqrt(n * p * (1 - p))
    k = floor(mu - sd)
    exact, approx = binom_cdf(k, n, p), NormalDist(mu, sd).cdf(k)
    print(f"{n:6d}  {p:5.2f}  {sd ** 2:7.2f}  {k:4d}   {exact:.4f}   {approx:.4f}   {abs(exact - approx):.4f}")

     n      p   np(1-p)    k    exact    normal   abs. error
    20   0.50     5.00     7   0.1316   0.0899   0.0417
   100   0.50    25.00    45   0.1841   0.1587   0.0254
  1000   0.50   250.00   484   0.1635   0.1558   0.0077
    20   0.10     1.80     0   0.1216   0.0680   0.0536
   100   0.10     9.00     7   0.2061   0.1587   0.0474
   100   0.02     1.96     0   0.1326   0.0766   0.0561
  1000   0.02    19.60    15   0.1539   0.1294   0.0245
  5000   0.02    98.00    90   0.1689   0.1562   0.0126

Two patterns stand out. The error shrinks as $np(1 - p)$ grows, but slowly, roughly like $1/\sqrt{np(1-p)}$. And skewed cases do worse for a similar variance: Binomial$(100, 0.1)$, with variance 9, has a larger error than Binomial$(20, 0.5)$, with variance 5. Much of the remaining error, especially for small $n$, comes from treating a discrete count as continuous, which the continuity correction fixes.

Key takeaways

• For large $n$, Binomial$(n, p) \approx N(np, np(1 - p))$, and $\hat p \approx N(p, p(1 - p)/n)$.
• Use it when $np \ge 10$ and $n(1 - p) \ge 10$; for rare events use the Poisson instead.
• Accuracy depends on the variance $np(1 - p)$ and on how skewed the binomial is.
• Exact binomial probabilities are cheap today; the approximation's value is in quick reasoning and in the theory of tests and intervals.

🧪 Practice

1. A fair coin is tossed 400 times. Approximate the probability that the number of heads is between 190 and 210, and compare with the exact value.
2. In an election with 52% support for candidate A, a poll asks 1,000 random voters. Approximate the probability that the poll shows A with less than 50%.
3. For which of these is the normal approximation reasonable: Binomial$(30, 0.5)$, Binomial$(30, 0.05)$, Binomial$(500, 0.05)$, Binomial$(10^6, 10^{-5})$? Which would you approximate by a Poisson instead?
4. Interview: A flight has 180 seats and each ticket holder independently fails to show up with probability 10%. How many tickets can the airline sell while keeping the probability of more passengers than seats below 5%? Hint: the number of passengers who show up is Binomial$(n, 0.9)$; write the condition with the normal approximation as an inequality in $n$, solve it, then check neighbouring values of $n$ with the exact binomial.

Continuity Correction

Why it matters. A binomial count only takes integer values, but the normal distribution is continuous. Under the normal approximation, $P(X = 50)$ would be 0, and $P(X \le 50)$ and $P(X < 51)$, which are equal for a count, would get different values. The continuity correction fixes this with a simple rule, and it often reduces the approximation error several-fold for moderate $n$.

The idea. Draw the PMF of the count as a histogram whose bars have width 1, centered on the integers. The bar for $k$ covers $[k - 0.5, k + 0.5]$, and its area is $P(X = k)$. The normal curve approximates the histogram, so the probability of a set of integers is the normal area over their bars, which extends half a unit beyond the outermost integers.

            +-----+
       +----+     +----+                  P(X <= 17) covers the bars of 0, 1, ..., 17,
  +----+    |     |    +----+             so it ends at the right edge of bar 17:
  |    |    |     |    |    |                  P(X <= 17)  ~  P(Y <= 17.5)
  | 15 | 16 | 17  | 18 | 19 |
  +----+----+-----+----+----+
 14.5 15.5 16.5  17.5 18.5 19.5
                  ^
                  cut here, not at 17

Rules. Let $X$ be integer-valued and $Y$ the approximating normal variable.

Event about $X$	Same event for integers	Normal approximation
$X \le k$	$X \le k$	$P(Y \le k + 0.5)$
$X < k$	$X \le k - 1$	$P(Y \le k - 0.5)$
$X \ge k$	$X \ge k$	$P(Y \ge k - 0.5)$
$X > k$	$X \ge k + 1$	$P(Y \ge k + 0.5)$
$X = k$	$X = k$	$P(k - 0.5 \le Y \le k + 0.5)$
$a \le X \le b$	$a \le X \le b$	$P(a - 0.5 \le Y \le b + 0.5)$

The safe procedure: first rewrite the event with $\le$ and $\ge$ on integers, then widen it by 0.5 on each open end.

from math import comb, exp, factorial, sqrt
from statistics import NormalDist

n, p = 40, 0.5
pmf = [comb(n, k) * p ** k * (1 - p) ** (n - k) for k in range(n + 1)]
Y = NormalDist(n * p, sqrt(n * p * (1 - p)))                 # N(20, 10)

events = [                                                     # (label, exact, without, with correction)
    ("P(X <= 17)",      sum(pmf[:18]),    Y.cdf(17),               Y.cdf(17.5)),
    ("P(X < 17)",       sum(pmf[:17]),    Y.cdf(17),               Y.cdf(16.5)),
    ("P(X >= 25)",      sum(pmf[25:]),    1 - Y.cdf(25),           1 - Y.cdf(24.5)),
    ("P(X = 20)",       pmf[20],          0.0,                     Y.cdf(20.5) - Y.cdf(19.5)),
    ("P(18 <= X <= 22)", sum(pmf[18:23]), Y.cdf(22) - Y.cdf(18),   Y.cdf(22.5) - Y.cdf(17.5)),
]
print("Binomial(40, 0.5)     exact    without   with correction")
for label, exact, raw, corr in events:
    print(f"{label:18s}   {exact:.4f}   {raw:.4f}    {corr:.4f}")

# The same correction works for other integer-valued variables: Poisson(20), P(X <= 15)
lam = 20
exact = sum(exp(-lam) * lam ** j / factorial(j) for j in range(16))
W = NormalDist(lam, sqrt(lam))
print(f"Poisson(20) P(X <= 15): exact = {exact:.4f}, without = {W.cdf(15):.4f}, with = {W.cdf(15.5):.4f}")

Binomial(40, 0.5)     exact    without   with correction
P(X <= 17)           0.2148   0.1714    0.2146
P(X < 17)            0.1341   0.1714    0.1342
P(X >= 25)           0.0769   0.0569    0.0774
P(X = 20)            0.1254   0.0000    0.1256
P(18 <= X <= 22)     0.5704   0.4729    0.5708
Poisson(20) P(X <= 15): exact = 0.1565, without = 0.1318, with = 0.1572

Without the correction, $P(X \le 17)$ and $P(X < 17)$ get the same value, which cannot be right. With it, every error falls below 0.001. For the skewed Poisson, the correction helps but cannot remove the error caused by skewness.

When to skip it. For very large $n$, half a unit is negligible relative to the standard deviation. The correction is specific to integer-valued (lattice) variables. In chi-square tests for $2 \times 2$ tables, the analogous Yates correction (Section 9.3) is known to be conservative and is often omitted.

Key takeaways

• Each integer $k$ corresponds to the interval $[k - 0.5, k + 0.5]$ under the approximating normal curve.
• Rewrite the event using $\le$ and $\ge$ on integers, then extend each end by 0.5.
• The correction makes $P(X = k)$ positive and distinguishes $<$ from $\le$; it typically cuts the error several-fold for moderate $n$.
• It matters most for small to moderate $n$ and becomes negligible as the standard deviation grows.

🧪 Practice

1. For $X \sim$ Binomial$(100, 0.5)$, approximate $P(X = 50)$ with the continuity correction and compare with the exact value $\binom{100}{50}/2^{100}$.
2. Write the corrected normal approximation for each event: $10 < X < 20$, $X \ne 15$, "at most 3 below the mean of 20".
3. For $X \sim$ Poisson$(30)$, approximate $P(X \ge 40)$ with and without the correction and compare with the exact value.
4. Interview: Your variable takes only even values (for example, twice a binomial count). What continuity correction would you use, and why does the standard $\pm 0.5$ fail? Hint: the correction is half the spacing between possible values; draw the bars of the histogram for the even values and see where their edges fall.

Assessing Normality and Q-Q Plots

Why it matters. Many classical methods assume normal data or normal errors: $t$-intervals for small samples, prediction intervals in regression, control charts, and tests for variances, which are especially sensitive. Real data are often skewed or heavy-tailed. Before relying on normality you need to check it, and the most informative tool for doing so is the quantile-quantile (Q-Q) plot.

Tools compared.

Tool	What it shows	Weakness
Histogram (Section 2.1)	overall shape	depends on bin choices; tails are hard to see
Boxplot (Section 2.4)	skew and outliers	hides multimodality; little detail in the tails
Skewness and kurtosis	one number each for asymmetry and tails	very noisy, driven by a few points (Section 2.4)
Q-Q plot	the whole distribution, especially the tails	needs some practice to read
Formal tests (Shapiro-Wilk, Anderson-Darling, Kolmogorov-Smirnov, Section 9.5)	a p-value	small samples: little power; large samples: reject trivial deviations

Constructing a normal Q-Q plot. If the data come from $N(\mu, \sigma^2)$, the $i$-th smallest of $n$ values should be close to the corresponding normal quantile.

1. Sort the data: $x_{(1)} \le x_{(2)} \le \dots \le x_{(n)}$.
2. Compute the theoretical quantiles $z_i = \Phi^{-1}\left(\frac{i - 0.5}{n}\right)$ (other plotting positions, such as $\frac{i - 3/8}{n + 1/4}$, differ slightly).
3. Plot the points $(z_i, x_{(i)})$.

If the data are normal, the points lie close to a straight line with intercept $\approx \mu$ and slope $\approx \sigma$. Because the line can have any intercept and slope, you do not need to know $\mu$ and $\sigma$ in advance. The correlation $r_{QQ}$ between $z_i$ and $x_{(i)}$ summarizes straightness, and the Shapiro-Wilk test is based on a similar idea.

Reading the patterns. Deviations show up at the ends, where the tails are:

  normal                right-skewed          heavy tails           light tails
  (straight line)       (convex: bends up)    (S-shape)             (inverted S)

  x |           .       x |           .       x |           .       x |       . . .
    |         .           |          .          |           .         |      .
    |       .             |                     |        .            |     .
    |     .               |       .             |   . . .             |    .
    |   .                 |    . .              |  .                  | . . .
    | .                   | . .                 |  .                  |
    +----------- z        +----------- z        +----------- z        +----------- z

Pattern	Meaning	Typical data
Straight line	consistent with normal	measurement errors, heights
Curve bending upward (convex)	right skew: the upper tail is stretched	incomes, latencies, counts
Curve bending downward (concave)	left skew	exam scores near a ceiling
S-shape: low end below the line, high end above	heavy tails on both sides	financial returns, $t$-distributed data
Inverted S	light tails	uniform or bounded data
A few points far off at the ends	outliers	data errors, rare events
Steps or horizontal runs	discrete or rounded data	integer counts, ratings

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(18)
n = 500
samples = {
    "normal":              rng.normal(50, 10, n),
    "lognormal (skewed)":  rng.lognormal(0, 0.6, n),
    "t, 3 df (heavy)":     rng.standard_t(3, n),
    "uniform (light)":     rng.uniform(0, 1, n),
}
z = np.array([NormalDist().inv_cdf((i - 0.5) / n) for i in range(1, n + 1)])   # theoretical quantiles

print("standardized sample quantile vs normal quantile at p:")
probs = [0.01, 0.05, 0.25, 0.50, 0.75, 0.95, 0.99]
print(f"{'':20s} r_QQ   " + "  ".join(f"{p:6.2f}" for p in probs))
print(f"{'normal quantile':20s}        " + "  ".join(f"{NormalDist().inv_cdf(p):+6.2f}" for p in probs))
for name, x in samples.items():
    xs = np.sort(x)
    r_qq = np.corrcoef(z, xs)[0, 1]                      # straightness of the Q-Q plot
    std = (xs - xs.mean()) / xs.std(ddof=1)
    qs = [std[int(p * n)] for p in probs]                # the point plotted at each probability
    print(f"{name:20s} {r_qq:.3f}  " + "  ".join(f"{q:+6.2f}" for q in qs))

standardized sample quantile vs normal quantile at p:
                     r_QQ     0.01    0.05    0.25    0.50    0.75    0.95    0.99
normal quantile              -2.33   -1.64   -0.67   +0.00   +0.67   +1.64   +2.33
normal               0.998   -2.12   -1.64   -0.66   -0.04   +0.65   +1.75   +2.68
lognormal (skewed)   0.912   -1.20   -1.03   -0.69   -0.23   +0.42   +1.84   +3.32
t, 3 df (heavy)      0.961   -3.10   -1.49   -0.45   +0.03   +0.45   +1.55   +2.75
uniform (light)      0.977   -1.63   -1.55   -0.85   -0.04   +0.88   +1.60   +1.74

Each row is a summary of a Q-Q plot: the normal sample matches the normal quantiles; the lognormal sample has a short lower tail and a long upper tail (convex); the $t$ sample is too extreme at both ends (S-shape); the uniform sample is not extreme enough (inverted S). To draw the actual plot, use scipy.stats.probplot or statsmodels.api.qqplot, or plot z against the sorted data with matplotlib.

Practical advice. Judge normality by its consequences. Mild skew is harmless for $t$-intervals on means with moderate $n$ (the central limit theorem helps), while heavy tails or outliers can ruin variance-based methods regardless of $n$. With small samples, Q-Q plots of truly normal data wiggle noticeably; comparing your plot with a few plots of simulated normal samples of the same size calibrates your eye.

Key takeaways

• A normal Q-Q plot plots sorted data against standard normal quantiles; normal data give a straight line with slope $\approx \sigma$.
• Convex curves mean right skew, S-shapes mean heavy tails, and isolated end points are outliers.
• Q-Q plots show how data deviate from normality, which histograms and tests do not.
• Formal normality tests answer a different question than "is normality good enough for my method?"; use them with care, especially for very small or very large samples.

🧪 Practice

1. Build a Q-Q plot by hand for the data 3.1, 4.7, 5.0, 5.9, 8.8: compute the plotting positions $(i - 0.5)/5$, their normal quantiles, and decide whether any point looks unusual.
2. What does a normal Q-Q plot of exponential data look like? Predict the shape, then check by simulation.
3. A Q-Q plot can compare any two distributions. Write code that compares a sample with the exponential distribution (using $Q(p) = -\ln(1 - p)$) and apply it to simulated interarrival times.
4. Interview: A Shapiro-Wilk test on 1 million latency measurements gives $p < 10^{-10}$, and a colleague concludes that no method assuming normality can be used. How do you respond? Hint: with huge $n$, any tiny deviation is significant; ask which quantity the method needs to be normal (the data, the errors, or the sampling distribution of a mean), what the Q-Q plot shows, and whether a log transform or a method designed for skewed data answers the real question better.

<a id="54-sampling-distributions-families"></a>

5.4 Sampling Distributions Families

When data are normal, the statistics computed from them (sample variances, standardized means, ratios of variances) follow three families built from the normal: chi-square, $t$ and $F$. They are the reference distributions for most classical confidence intervals and tests in Chapters 8 to 10.

Chi-Square Distribution

Why it matters. Squared quantities are everywhere in statistics: variances are averages of squared deviations, and many test statistics are sums of squared standardized differences. The chi-square distribution describes a sum of squared independent standard normals. It gives the sampling distribution of the sample variance, which leads to confidence intervals for variances (Section 8.3), and it is the reference distribution for goodness-of-fit and independence tests (Section 9.3) and for likelihood ratio tests (Section 9.6).

Definition. If $Z_1, \dots, Z_k$ are independent $N(0, 1)$, then

$$Q = Z_1^2 + Z_2^2 + \dots + Z_k^2 \sim \chi^2_k,$$

the chi-square distribution with $k$ degrees of freedom. Its density is

$$f(x) = \frac{x^{k/2 - 1}e^{-x/2}}{2^{k/2},\Gamma(k/2)}, \qquad x > 0,$$

which is Gamma$(k/2, 1/2)$ (rate 1/2, Section 5.2). Section 4.4 derived the case $k = 1$ with the CDF method.

Properties.

Property	Value
Mean	$k$ (each $Z_i^2$ has mean 1)
Variance	$2k$ (each $Z_i^2$ has variance $E[Z^4] - 1 = 3 - 1 = 2$)
Additivity	independent $\chi^2_{k_1} + \chi^2_{k_2} = \chi^2_{k_1 + k_2}$
Special case	$\chi^2_2$ is Exponential with mean 2
Shape	right-skewed; skewness $\sqrt{8/k}$, so approximately $N(k, 2k)$ for large $k$
Better approximation (Wilson-Hilferty)	$(Q/k)^{1/3} \approx N\left(1 - \frac{2}{9k}, \frac{2}{9k}\right)$

Degrees of freedom. The parameter counts the independent squared normal pieces. This becomes concrete with the sample variance: for $X_1, \dots, X_n$ independent $N(\mu, \sigma^2)$,

$$\frac{(n - 1)S^2}{\sigma^2} = \sum_{i=1}^n \left(\frac{X_i - \bar X}{\sigma}\right)^2 \sim \chi^2_{n - 1}.$$

There are $n$ squared terms but only $n - 1$ degrees of freedom, because the deviations $X_i - \bar X$ always sum to zero: once $n - 1$ of them are known, the last is determined. Estimating the mean "uses up" one degree of freedom, which is also why $S^2$ divides by $n - 1$ (Section 2.3).

Geometric view. $Q$ is the squared length of a random vector $(Z_1, \dots, Z_k)$. Its mean $k$ and standard deviation $\sqrt{2k}$ show that in high dimensions, the length is concentrated near $\sqrt{k}$: almost all of the probability of a high-dimensional standard normal lies in a thin shell, far from the origin where the density is highest.

import numpy as np
from math import exp, factorial, sqrt
from statistics import NormalDist

def chi2_sf_even(x, k):
    """P(Q > x) for even k, via the Gamma-Poisson identity of Section 5.2."""
    m = k // 2
    return sum(exp(-x / 2) * (x / 2) ** j / factorial(j) for j in range(m))

def chi2_quantile_even(p, k):
    lo, hi = 0.0, 10.0 * k + 100
    for _ in range(100):                                    # bisection on the CDF
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if 1 - chi2_sf_even(mid, k) < p else (lo, mid)
    return lo

def wilson_hilferty(p, k):
    z = NormalDist().inv_cdf(p)
    return k * (1 - 2 / (9 * k) + z * sqrt(2 / (9 * k))) ** 3

rng = np.random.default_rng(19)
print(" k   95% quantile: exact   Wilson-Hilferty   N(k, 2k)   simulated sum of k Z^2")
for k in [2, 4, 10, 20, 50]:
    q = (rng.standard_normal((200_000, k)) ** 2).sum(axis=1)
    print(f"{k:2d}               {chi2_quantile_even(0.95, k):7.3f}   {wilson_hilferty(0.95, k):7.3f}"
          f"        {k + 1.645 * sqrt(2 * k):7.3f}    {np.quantile(q, 0.95):7.3f}")

# Sample variance: n = 5 observations from N(10, 2^2); (n - 1) S^2 / sigma^2 should be chi-square with 4 df
n, sigma = 5, 2.0
s2 = rng.normal(10, sigma, size=(200_000, n)).var(axis=1, ddof=1)
stat = (n - 1) * s2 / sigma ** 2
print(f"(n-1)S^2/sigma^2: mean = {stat.mean():.3f} (4), variance = {stat.var():.3f} (8), "
      f"P(> 9.488) = {np.mean(stat > chi2_quantile_even(0.95, 4)):.4f} (0.05)")
print(f"S^2 itself: mean = {s2.mean():.3f} (sigma^2 = 4), but P(S^2 < 2) = {np.mean(s2 < 2):.3f}: "
      f"with 5 points, S^2 is often far from sigma^2")

 k   95% quantile: exact   Wilson-Hilferty   N(k, 2k)   simulated sum of k Z^2
 2                 5.991     5.937          5.290      5.983
 4                 9.488     9.456          8.653      9.485
10                18.307    18.292         17.357     18.293
20                31.410    31.402         30.404     31.392
50                67.505    67.501         66.450     67.504
(n-1)S^2/sigma^2: mean = 3.992 (4), variance = 7.952 (8), P(> 9.488) = 0.0491 (0.05)
S^2 itself: mean = 3.992 (sigma^2 = 4), but P(S^2 < 2) = 0.265: with 5 points, S^2 is often far from sigma^2

The normal approximation $N(k, 2k)$ is poor even at $k = 50$ because of the right skew, while the Wilson-Hilferty cube-root approximation is within 0.06 of the exact quantile even at $k = 2$, and much closer for larger $k$. The last line is a warning: variances estimated from small samples are very imprecise.

Key takeaways

• $\chi^2_k$ is the distribution of a sum of $k$ independent squared standard normals: mean $k$, variance $2k$, right-skewed.
• It is Gamma$(k/2, 1/2)$, and chi-squares add their degrees of freedom.
• For normal data, $(n - 1)S^2/\sigma^2 \sim \chi^2_{n-1}$; one degree of freedom is lost by estimating the mean.
• It is the reference distribution for variance intervals, chi-square tests and likelihood ratio tests.

🧪 Practice

1. Use the definition to find the mean and variance of $\chi^2_{10}$, and the probability that it exceeds 18.31 (use the even-$k$ formula in the code).
2. Show that $\chi^2_2$ is the exponential distribution with mean 2, both from the gamma form and from chi2_sf_even.
3. A sample of $n = 10$ normal observations has $s^2 = 9.5$. If the true variance were $\sigma^2 = 4$, would this value of $(n - 1)s^2/\sigma^2$ be unusual? Use the simulation to estimate the probability of a value at least this large.
4. Interview: Where does most of the probability mass of a 1,000-dimensional standard normal vector lie? Why does this matter for nearest-neighbor methods and for sampling algorithms in high dimensions? Hint: the squared length is $\chi^2_{1000}$; compute its mean and SD to find the typical length and how much it varies, and compare with the location of the density's peak.

Student's t-Distribution

Why it matters. To judge how far a sample mean is from a hypothesized value, we divide the difference by its standard error $\sigma/\sqrt n$. But $\sigma$ is almost never known, and replacing it by the sample SD $S$ adds extra randomness: the ratio is no longer normal. William Gosset, working on small samples at the Guinness brewery and publishing as "Student" (1908), found its exact distribution. The $t$-distribution is behind the $t$-intervals and $t$-tests of Chapters 8 and 9 and the coefficient tests of regression (Chapter 10). It is also a useful heavy-tailed model in its own right.

Definition. If $Z \sim N(0, 1)$ and $V \sim \chi^2_\nu$ are independent, then

$$T = \frac{Z}{\sqrt{V/\nu}} \sim t_\nu,$$

the $t$-distribution with $\nu$ degrees of freedom, with density

$$f(t) = \frac{\Gamma\left(\frac{\nu + 1}{2}\right)}{\sqrt{\nu\pi},\Gamma\left(\frac{\nu}{2}\right)}\left(1 + \frac{t^2}{\nu}\right)^{-(\nu + 1)/2}.$$

The key application. For $X_1, \dots, X_n$ independent $N(\mu, \sigma^2)$,

$$T = \frac{\bar X - \mu}{S/\sqrt{n}} \sim t_{n - 1}.$$

This works because, for normal samples, $\bar X$ and $S^2$ are independent (a special property of the normal, see the last topic of this subchapter). Then $\frac{\bar X - \mu}{\sigma/\sqrt n}$ is the $Z$, and $\frac{(n-1)S^2}{\sigma^2}$ is the $\chi^2_{n-1}$ variable $V$.

Properties.

Property	Value
Symmetry	symmetric about 0, bell-shaped
Mean	0 (for $\nu > 1$)
Variance	$\frac{\nu}{\nu - 2}$ (for $\nu > 2$), always above 1
Tails	polynomial, like $
$\nu = 1$	the Cauchy distribution: no mean
$\nu \to \infty$	converges to $N(0, 1)$

Why the tails are heavier. In a small sample, $S$ is noisy. When it happens to underestimate $\sigma$, the ratio becomes too large. The $t$-distribution accounts for this extra uncertainty by putting more probability in the tails, so its critical values are larger than the normal's:

$\nu$	1	2	4	5	10	20	30	100	$\infty$
$t_{0.975,,\nu}$	12.706	4.303	2.776	2.571	2.228	2.086	2.042	1.984	1.960

Beyond about 30 degrees of freedom, the difference from 1.96 is small, which is where the old "use $z$ for $n > 30$" rule comes from.

import numpy as np

rng = np.random.default_rng(20)
N = 1_000_000

# T = Z / sqrt(V / nu): tail probability beyond the normal's 1.96
print(" nu   P(|T| > 1.96)   via rng.standard_t   variance (theory)")
for nu in [1, 2, 5, 10, 30]:
    t = rng.standard_normal(N) / np.sqrt(rng.chisquare(nu, N) / nu)
    var_theory = f"{nu / (nu - 2):.3f}" if nu > 2 else "infinite"
    print(f"{nu:3d}       {np.mean(np.abs(t) > 1.96):.4f}           {np.mean(np.abs(rng.standard_t(nu, N)) > 1.96):.4f}"
          f"            {var_theory}")

# Coverage of 95% intervals for a mean from n = 5 normal observations (true mean 100, sigma 15)
n, mu, sigma = 5, 100.0, 15.0
x = rng.normal(mu, sigma, size=(200_000, n))
xbar, s = x.mean(axis=1), x.std(axis=1, ddof=1)
cover = lambda half: np.mean(np.abs(xbar - mu) <= half)
print(f"known sigma, xbar +- 1.96 sigma/sqrt(n):   {cover(1.96 * sigma / np.sqrt(n)):.4f}")
print(f"estimated s, xbar +- 1.96 s/sqrt(n):       {cover(1.96 * s / np.sqrt(n)):.4f}   <- too short")
print(f"estimated s, xbar +- 2.776 s/sqrt(n) (t_4): {cover(2.776 * s / np.sqrt(n)):.4f}")

 nu   P(|T| > 1.96)   via rng.standard_t   variance (theory)
  1       0.3001           0.3004            infinite
  2       0.1891           0.1888            infinite
  5       0.1073           0.1069            1.667
 10       0.0786           0.0786            1.250
 30       0.0595           0.0592            1.071
known sigma, xbar +- 1.96 sigma/sqrt(n):   0.9502
estimated s, xbar +- 1.96 s/sqrt(n):       0.8785   <- too short
estimated s, xbar +- 2.776 s/sqrt(n) (t_4): 0.9503

Using 1.96 with an estimated SD on 5 observations gives an interval that is advertised as 95% but covers the true mean only about 88% of the time. The $t$ critical value restores the promised coverage.

Key takeaways

• $t_\nu = Z/\sqrt{\chi^2_\nu/\nu}$; for normal data, $(\bar X - \mu)/(S/\sqrt n) \sim t_{n-1}$.
• It is symmetric with heavier tails than the normal; its variance $\nu/(\nu - 2)$ exceeds 1.
• Use $t$ critical values whenever $\sigma$ is estimated; they approach 1.96 as $\nu$ grows.
• $t_1$ is the Cauchy distribution; small $\nu$ gives a useful heavy-tailed model for robust methods.

🧪 Practice

1. Compute the variance of $t_5$ and of $t_{30}$, and explain why $t_2$ has infinite variance.
2. Show that the $t_1$ density is $\frac{1}{\pi(1 + t^2)}$, the Cauchy density.
3. A sample of $n = 12$ has $\bar x = 52.3$ and $s = 4.1$. Compute the $t$ statistic for testing $\mu = 50$ and compare it with the critical value $t_{0.975, 11} = 2.201$.
4. Interview: With $n = 200$ observations, does it matter whether you use the $t$ or the normal distribution for a confidence interval for the mean? What matters more at that sample size? Hint: compare $t_{0.975, 199} \approx 1.972$ with 1.96; then think about what the derivation of the $t$-distribution assumed about the data, and when the central limit theorem rescues the interval and when it does not (heavy tails, dependence).

F-Distribution

Why it matters. Many questions compare two sources of variability: is one machine less consistent than another, are differences between group means larger than the noise within groups, does adding predictors to a regression explain significantly more variance? Each is answered with a ratio of two variance estimates, and for normal data such a ratio follows an $F$-distribution (named after R. A. Fisher). It is the reference distribution of the ANOVA $F$-test (Section 9.4), the overall regression $F$-test (Section 10.2) and tests comparing nested models.

Definition. If $U \sim \chi^2_{d_1}$ and $V \sim \chi^2_{d_2}$ are independent, then

$$F = \frac{U/d_1}{V/d_2} \sim F(d_1, d_2),$$

with $d_1$ numerator and $d_2$ denominator degrees of freedom. Each chi-square is divided by its degrees of freedom, so numerator and denominator both have mean 1, and $F$ is centered near 1.

The key application. For independent normal samples of sizes $n_1$ and $n_2$,

$$\frac{S_1^2/\sigma_1^2}{S_2^2/\sigma_2^2} \sim F(n_1 - 1, ; n_2 - 1).$$

If the population variances are equal, the ratio of sample variances $S_1^2/S_2^2$ is itself $F$-distributed, which gives a test for equal variances.

Properties.

Property	Value
Support	$(0, \infty)$, right-skewed
Mean	$\frac{d_2}{d_2 - 2}$ for $d_2 > 2$ (slightly above 1)
Reciprocal	$1/F \sim F(d_2, d_1)$, so lower quantiles come from upper ones: $F_{p}(d_1, d_2) = 1/F_{1-p}(d_2, d_1)$
Relation to $t$	$T^2 \sim F(1, \nu)$ when $T \sim t_\nu$
Large $d_2$	$d_1 F \to \chi^2_{d_1}$ (the denominator becomes the constant 1)

A warning. The $F$-test for equal variances is very sensitive to non-normality: with heavy-tailed data it rejects far too often, even when the variances are equal. This is because, unlike for means, there is no central limit theorem rescue for variance ratios; their distribution depends on the kurtosis of the data. Levene's or the Brown-Forsythe test (Section 9.2) is preferred in practice. The ANOVA $F$-test for means is much more robust.

import numpy as np

rng = np.random.default_rng(21)
N = 400_000
n1, n2 = 10, 15                                           # two machines, 10 and 15 bottles measured
d1, d2 = n1 - 1, n2 - 1

# Ratio of sample variances from two normal populations with the same variance
F = rng.normal(0, 3, (N, n1)).var(axis=1, ddof=1) / rng.normal(0, 3, (N, n2)).var(axis=1, ddof=1)
crit_hi = np.quantile(F, 0.95)
print(f"S1^2/S2^2: mean = {F.mean():.4f} (d2/(d2-2) = {d2 / (d2 - 2):.4f}), 95th percentile = {crit_hi:.3f}")

# Reciprocal property: the 5th percentile of F(9, 14) is 1 / the 95th percentile of F(14, 9)
F_swapped = (rng.chisquare(d2, N) / d2) / (rng.chisquare(d1, N) / d1)
print(f"5th percentile = {np.quantile(F, 0.05):.4f}, 1 / 95th percentile of F(14, 9) = {1 / np.quantile(F_swapped, 0.95):.4f}")

# T^2 ~ F(1, nu)
nu = 8
t2 = rng.standard_t(nu, N) ** 2
f1 = rng.chisquare(1, N) / (rng.chisquare(nu, N) / nu)
print(f"95th percentiles: T^2 with T ~ t_8 = {np.quantile(t2, 0.95):.3f}, F(1, 8) = {np.quantile(f1, 0.95):.3f}, "
      f"t_0.975,8 squared = {2.306 ** 2:.3f}")

# Sensitivity to non-normality: equal variances, but heavy-tailed or skewed data
def false_alarm_rate(sampler):
    ratio = sampler((N, n1)).var(axis=1, ddof=1) / sampler((N, n2)).var(axis=1, ddof=1)
    return np.mean(ratio > crit_hi)

print("P(ratio > normal-theory 95th percentile) when the variances are EQUAL:")
print(f"  normal data       {false_alarm_rate(lambda size: rng.normal(0, 1, size)):.3f}")
print(f"  t with 5 df       {false_alarm_rate(lambda size: rng.standard_t(5, size)):.3f}")
print(f"  exponential       {false_alarm_rate(lambda size: rng.exponential(1, size)):.3f}")

S1^2/S2^2: mean = 1.1673 (d2/(d2-2) = 1.1667), 95th percentile = 2.648
5th percentile = 0.3312, 1 / 95th percentile of F(14, 9) = 0.3310
95th percentiles: T^2 with T ~ t_8 = 5.289, F(1, 8) = 5.343, t_0.975,8 squared = 5.318
P(ratio > normal-theory 95th percentile) when the variances are EQUAL:
  normal data       0.050
  t with 5 df       0.104
  exponential       0.156

Key takeaways

• $F(d_1, d_2)$ is a ratio of two independent chi-squares, each divided by its degrees of freedom.
• For normal samples, the ratio of sample variances (scaled by the true variances) is $F(n_1 - 1, n_2 - 1)$.
• $1/F \sim F(d_2, d_1)$ and $t_\nu^2 \sim F(1, \nu)$.
• It underlies ANOVA and regression $F$-tests; the $F$-test for equal variances is unreliable for non-normal data.

🧪 Practice

1. Compute the mean of $F(5, 10)$ and of $F(5, 100)$. Why does the mean not depend on $d_1$?
2. Prove that $T^2 \sim F(1, \nu)$ from the definitions of $t_\nu$ and $F$.
3. Two production lines give $s_1^2 = 12.4$ from $n_1 = 10$ items and $s_2^2 = 4.1$ from $n_2 = 15$ items. Compute the variance ratio and compare it with the simulated 95th percentile of $F(9, 14)$ from the code. What do you conclude, and what assumption must you check?
4. Interview: In one-way ANOVA the alternative hypothesis "the group means are not all equal" has no direction, yet only large values of $F$ lead to rejection. Why is the test one-sided? Hint: write $F$ as (variance estimate based on the spread of group means) / (variance estimate within groups); ask what differences between means do to the numerator, and what a very small $F$ would suggest instead.

Relationships Among Normal-Derived Distributions

Why it matters. The normal, chi-square, $t$ and $F$ distributions, together with the gamma, exponential and beta from Section 5.2, form a single family tree. Knowing the links lets you derive one distribution from another instead of memorizing each, check software output, understand the approximations used for large samples, and see why the classical tests depend so much on normality.

The family tree.

                               Z ~ N(0, 1)
                                    |
                 square and add k   |   Z_1^2 + ... + Z_k^2   (k = 1: Z^2 ~ chi^2_1)
                                    v
   Gamma(k/2, rate 1/2)   =      chi^2_k      (k = 2: Exponential with mean 2)
                                    |
              +---------------------+----------------------+
              |                                            |
   t_nu = Z / sqrt(chi^2_nu / nu)          F(d1, d2) = (chi^2_d1 / d1) / (chi^2_d2 / d2)
              |                                            |
              |-- square -->  t_nu^2 = F(1, nu)            |-- d2 -> infinity -->  d1 F ~ chi^2_d1
              |                                            |
              v  nu -> infinity                            v  d1 F / (d1 F + d2)
           N(0, 1)                                      Beta(d1/2, d2/2)

Relationships at a glance.

Relationship	Statement	Where it is used
Normal $\to$ chi-square	$Z_1^2 + \dots + Z_k^2 \sim \chi^2_k$	variance estimates, chi-square tests
Chi-square = gamma	$\chi^2_k = $ Gamma$(k/2, 1/2)$; $\chi^2_2 = $ Exponential(mean 2)	computing chi-square probabilities
Normal and chi-square $\to t$	$Z/\sqrt{\chi^2_\nu/\nu} \sim t_\nu$	$t$-tests and intervals
Chi-squares $\to F$	$(\chi^2_{d_1}/d_1)/(\chi^2_{d_2}/d_2) \sim F(d_1, d_2)$	ANOVA, regression $F$-tests
$t \to F$	$t_\nu^2 \sim F(1, \nu)$	a two-group ANOVA equals a two-sided $t$-test
$F \to$ beta	$\frac{d_1F}{d_1F + d_2} \sim$ Beta$(d_1/2, d_2/2)$	$R^2$ in regression is beta under the null
Large-sample limits	$t_\nu \to N(0, 1)$; $d_1F(d_1, d_2) \to \chi^2_{d_1}$ as $d_2 \to \infty$	$z$-tests and chi-square tests as approximations

The special role of normality. For a normal sample, $\bar X$ and $S^2$ are independent, and $(n - 1)S^2/\sigma^2 \sim \chi^2_{n-1}$ exactly (Cochran's theorem, proved with the multivariate normal in Section 6.4). This independence is what turns $(\bar X - \mu)/(S/\sqrt n)$ into a $t$ variable. It is a characterization: the normal is the only distribution for which the sample mean and sample variance are independent. For skewed data they are correlated (a large mean often comes with a large variance), and the exact $t$ and $F$ results no longer hold.

import numpy as np

rng = np.random.default_rng(22)
N = 1_000_000
Z = lambda size: rng.standard_normal(size)
chi2 = lambda k: rng.chisquare(k, N)
qs = [0.50, 0.90, 0.99]

def compare(label, left, right):
    l, r = np.quantile(left, qs), np.quantile(right, qs)
    print(f"{label:38s} " + "  ".join(f"{a:7.3f} vs {b:7.3f}" for a, b in zip(l, r)))

print(f"{'relationship (quantiles 0.50, 0.90, 0.99)':38s}")
compare("Z1^2 + Z2^2 + Z3^2  vs  chi^2_3", (Z((N, 3)) ** 2).sum(axis=1), chi2(3))
compare("chi^2_2 / 2  vs  Exponential(1)", chi2(2) / 2, rng.exponential(1, N))
compare("t_8 squared  vs  F(1, 8)", rng.standard_t(8, N) ** 2, chi2(1) / (chi2(8) / 8))
compare("5 * F(5, 10000)  vs  chi^2_5", 5 * (chi2(5) / 5) / (chi2(10_000) / 10_000), chi2(5))
F = (chi2(4) / 4) / (chi2(6) / 6)
compare("4F / (4F + 6)  vs  Beta(2, 3)", 4 * F / (4 * F + 6), rng.beta(2, 3, N))
compare("t_200  vs  N(0, 1)", rng.standard_t(200, N), Z(N))

# Independence of the sample mean and variance: only for normal data
for name, sampler in [("normal", lambda s: rng.normal(0, 1, s)), ("exponential", lambda s: rng.exponential(1, s))]:
    x = sampler((200_000, 5))
    r = np.corrcoef(x.mean(axis=1), x.var(axis=1, ddof=1))[0, 1]
    print(f"corr(sample mean, sample variance), n = 5, {name:11s} data: {r:+.3f}")

relationship (quantiles 0.50, 0.90, 0.99)
Z1^2 + Z2^2 + Z3^2  vs  chi^2_3          2.365 vs   2.363    6.254 vs   6.240   11.360 vs  11.321
chi^2_2 / 2  vs  Exponential(1)          0.693 vs   0.693    2.307 vs   2.299    4.611 vs   4.588
t_8 squared  vs  F(1, 8)                 0.497 vs   0.497    3.461 vs   3.451   11.334 vs  11.212
5 * F(5, 10000)  vs  chi^2_5             4.348 vs   4.358    9.230 vs   9.233   15.083 vs  15.090
4F / (4F + 6)  vs  Beta(2, 3)            0.386 vs   0.386    0.679 vs   0.679    0.859 vs   0.859
t_200  vs  N(0, 1)                      -0.001 vs   0.000    1.288 vs   1.282    2.342 vs   2.331
corr(sample mean, sample variance), n = 5, normal      data: -0.003
corr(sample mean, sample variance), n = 5, exponential data: +0.685

Every pair of quantiles agrees to within simulation error, and the sample mean and variance are uncorrelated for normal data but strongly correlated for exponential data.

Key takeaways

• Chi-square, $t$ and $F$ are all built from independent standard normals: squares and sums give $\chi^2$, a normal over a scaled $\chi^2$ gives $t$, a ratio of scaled $\chi^2$s gives $F$.
• Useful identities: $\chi^2_k = $ Gamma$(k/2, 1/2)$, $t_\nu^2 = F(1, \nu)$, $1/F(d_1, d_2) = F(d_2, d_1)$.
• In the large-sample limits, $t \to N(0, 1)$ and $d_1F \to \chi^2_{d_1}$, which is why $z$- and chi-square approximations work for large samples.
• The exact results rely on the independence of $\bar X$ and $S^2$, which holds only for normal data.

🧪 Practice

1. Starting from $E[Z^2] = 1$ and $E[Z^4] = 3$, derive the mean and variance of $\chi^2_k$.
2. Show that if $X \sim F(d_1, d_2)$ then $1/X \sim F(d_2, d_1)$, directly from the definition.
3. Express an Exponential$(\lambda)$ variable as a multiple of a $\chi^2_2$ variable, and use this to write a 95% interval for the mean of an exponential distribution based on one observation.
4. Interview: Your environment only provides a standard normal random number generator. How would you generate chi-square, $t$ and $F$ random variables with integer degrees of freedom? What is the cost, and how would you handle 10,000 degrees of freedom efficiently? Hint: follow the definitions using sums of squares and ratios; for large degrees of freedom, use the gamma representation or a normal approximation of the chi-square instead of summing 10,000 squares.

<a id="55-stochastic-processes-preview"></a>

5.5 Stochastic Processes Preview

A stochastic process is a collection of random variables ${X_t}$ indexed by time, describing how randomness unfolds step by step. This preview shows that many distributions of this chapter are snapshots of a few basic processes (Bernoulli and Poisson processes for arrivals, Markov chains for systems with memory of their current state, random walks for accumulated noise), which reappear in time series (Section 12.1), MCMC (Section 11.3) and survival analysis (Section 12.3).

Bernoulli Processes

Why it matters. A Bernoulli process is the simplest model of events in discrete time: at each step something happens or not, independently, with the same probability. Website sessions that convert or not, daily server failures, coin tosses: all are Bernoulli processes. Its value is unification: the Bernoulli, binomial, geometric and negative binomial distributions of Section 5.1 are all questions about one sequence.

Definition. A Bernoulli process with parameter $p$ is a sequence $X_1, X_2, X_3, \dots$ of independent Bernoulli$(p)$ random variables. Think of $X_t = 1$ as "an arrival at time $t$".

  time      1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
  X_t       0  0  1  0  0  0  0  1  1  0  0  0  0  0  1  0
                  ^              ^  ^                 ^
  arrival         1st            2nd 3rd              4th

  successes in the first 10 steps:   3          Binomial(10, p)
  time of the 1st arrival:           3          Geometric(p)
  gaps between arrivals:             3, 5, 1, 6  i.i.d. Geometric(p)
  time of the 3rd arrival:           9          Negative binomial(3, p)

One process, four distributions.

Question about the sequence	Distribution
Is there an arrival at time $t$?	Bernoulli$(p)$
How many arrivals in $n$ steps?	Binomial$(n, p)$
When is the first arrival? How long is each gap?	Geometric$(p)$, with independent gaps
When is the $r$-th arrival?	Negative binomial$(r, p)$ (sum of $r$ gaps)

Fresh-start property. Because the trials are independent, the process "restarts" at any fixed time, and also at any time determined by the past alone, such as the moment of the 5th arrival: the future after that point is again a Bernoulli process, independent of the past. This is memorylessness (Section 5.1) applied to the whole process.

Merging and splitting.

• Splitting: if each arrival is independently sent to stream A with probability $q$, stream A is a Bernoulli process with parameter $pq$.
• Merging: if two independent Bernoulli processes with parameters $p_1, p_2$ are combined by recording an arrival when either has one, the result is a Bernoulli process with parameter $1 - (1 - p_1)(1 - p_2) = p_1 + p_2 - p_1p_2$.

import numpy as np

rng = np.random.default_rng(23)
p = 0.2
x = rng.random(1_000_000) < p                          # one long Bernoulli process

print("first 40 steps:", "".join("1" if v else "0" for v in x[:40]))

counts = x.reshape(-1, 20).sum(axis=1)                 # arrivals in consecutive windows of 20 steps
print(f"arrivals per 20 steps: mean = {counts.mean():.3f} (np = 4), variance = {counts.var():.3f} (np(1-p) = 3.2)")

times = np.flatnonzero(x) + 1                          # time index of each arrival
gaps = np.diff(np.concatenate([[0], times]))
print(f"gaps: mean = {gaps.mean():.3f} (1/p = 5), variance = {gaps.var():.2f} ((1-p)/p^2 = 20), "
      f"corr(gap, next gap) = {np.corrcoef(gaps[:-1], gaps[1:])[0, 1]:+.4f}")

g3 = gaps[: len(gaps) // 3 * 3].reshape(-1, 3).sum(axis=1)   # sums of 3 consecutive gaps
print(f"waiting time for 3 arrivals: mean = {g3.mean():.3f} (r/p = 15), variance = {g3.var():.2f} (r(1-p)/p^2 = 60)")

# Fresh start: after a run of 4 failures, the next step is still a success with probability p
after_run = x[4:][~x[:-4] & ~x[1:-3] & ~x[2:-2] & ~x[3:-1]]
print(f"P(arrival | previous 4 steps empty) = {after_run.mean():.4f} (p = {p})")

# Splitting: each arrival goes to stream A with probability 0.3
stream_a = x & (rng.random(x.size) < 0.3)
print(f"stream A: rate = {stream_a.mean():.4f} (pq = {p * 0.3:.2f})")

first 40 steps: 0011000000001000000010000100000000001000
arrivals per 20 steps: mean = 4.004 (np = 4), variance = 3.196 (np(1-p) = 3.2)
gaps: mean = 4.994 (1/p = 5), variance = 20.03 ((1-p)/p^2 = 20), corr(gap, next gap) = +0.0015
waiting time for 3 arrivals: mean = 14.983 (r/p = 15), variance = 60.75 (r(1-p)/p^2 = 60)
P(arrival | previous 4 steps empty) = 0.1999 (p = 0.2)
stream A: rate = 0.0600 (pq = 0.06)

Key takeaways

• A Bernoulli process is a sequence of independent Bernoulli$(p)$ trials, the discrete-time model of random arrivals.
• Counts in fixed windows are binomial, gaps are independent geometrics, and the time of the $r$-th arrival is negative binomial.
• The process starts afresh after any fixed time or arrival time: the past does not predict the future.
• Splitting multiplies the parameter by the routing probability; merging gives $p_1 + p_2 - p_1p_2$.

🧪 Practice

1. With $p = 0.3$, find the probability that the first arrival is at time 4 and the second at time 9.
2. Given that there were exactly 3 arrivals in the first 10 steps, show that every set of 3 positions is equally likely. What is the probability that the first step is one of them?
3. A website gets a purchase in each minute with probability 0.1, and 40% of purchases are from returning customers. Describe the process of purchases by new customers, and the expected waiting time between them.
4. Interview: In 100 tosses of a fair coin, how long is the longest run of heads, typically? People asked to write down "random" sequences almost never include runs that long; how could you use that to tell real sequences from invented ones? Hint: simulate it, then relate the answer to $\log_2 n$; a run of length $L$ starting at a given position has probability $2^{-L}$, and there are about $n$ starting positions.

Poisson Processes

Why it matters. Most real arrivals happen in continuous time: requests hit a server at any instant, not at discrete steps. The Poisson process is the continuous-time version of the Bernoulli process and the standard model for random arrivals: customers, calls, network packets, insurance claims, earthquakes, photons. It is the foundation of queueing theory, and it explains why counts are Poisson (Section 5.1) and waiting times exponential (Section 5.2).

From Bernoulli to Poisson. Divide time into tiny slots of length $\delta$, and let each slot independently contain an arrival with probability $\lambda\delta$. In a time interval of length $t$ there are $t/\delta$ slots, so the count is Binomial$(t/\delta, \lambda\delta)$, which tends to Poisson$(\lambda t)$ as $\delta \to 0$ (the Poisson approximation of Section 5.1). The geometric gaps, measured in time, tend to exponential gaps.

Definition. A Poisson process with rate $\lambda > 0$ is a counting process ${N(t), t \ge 0}$, where $N(t)$ is the number of arrivals in $[0, t]$, such that

1. $N(0) = 0$;
2. independent increments: counts in disjoint time intervals are independent;
3. stationary increments: the count in any interval of length $t$ is Poisson$(\lambda t)$.

An equivalent description: in a short interval of length $h$, one arrival occurs with probability $\approx \lambda h$, and two or more with probability negligible compared with $h$.

Bernoulli process (discrete time)	Poisson process (continuous time)
parameter $p$ per step	rate $\lambda$ per unit of time
arrivals in $n$ steps: Binomial$(n, p)$	arrivals in time $t$: Poisson$(\lambda t)$
gap: Geometric$(p)$, mean $1/p$	gap: Exponential$(\lambda)$, mean $1/\lambda$
time of $r$-th arrival: negative binomial	time of $r$-th arrival: Gamma$(r, \lambda)$
memoryless, fresh start	memoryless, fresh start

Properties.

Property	Statement
Superposition	merging independent Poisson processes with rates $\lambda_1, \lambda_2$ gives a Poisson process with rate $\lambda_1 + \lambda_2$
Thinning	keeping each arrival independently with probability $q$ gives a Poisson process with rate $q\lambda$, independent of the discarded arrivals
Conditional uniformity	given $N(t) = n$, the arrival times are distributed like $n$ independent Uniform$(0, t)$ points, sorted

Conditional uniformity gives the simplest way to simulate the process on $[0, T]$: draw $N \sim$ Poisson$(\lambda T)$, then $N$ uniform times, and sort them.

Non-homogeneous Poisson processes. Real traffic varies over the day. If the rate is a function $\lambda(t)$, the count in $[a, b]$ is Poisson with mean $\int_a^b \lambda(t),dt$, and counts in disjoint intervals are still independent. This is the usual model for arrivals with daily or weekly cycles.

import numpy as np

rng = np.random.default_rng(24)
lam, T = 4.0, 50_000.0                             # 4 calls per minute, observed for 50,000 minutes

# Simulation by conditional uniformity: a Poisson number of arrivals, placed uniformly at random
n = rng.poisson(lam * T)
times = np.sort(rng.uniform(0, T, n))

per_minute = np.bincount(times.astype(int), minlength=int(T))      # counts in each minute
print(f"calls per minute: mean = {per_minute.mean():.3f}, variance = {per_minute.var():.3f} (both 4)")
print(f"corr(count in minute i, count in minute i+1) = {np.corrcoef(per_minute[:-1], per_minute[1:])[0, 1]:+.4f}")

half_minute = np.bincount((times * 2).astype(int), minlength=int(2 * T))
print(f"P(no call in 30 s) = {np.mean(half_minute == 0):.4f} (e^-2 = {np.exp(-2):.4f})")

# Superposition: a second, independent stream of 1 call per minute
other = np.sort(rng.uniform(0, T, rng.poisson(1.0 * T)))
merged = np.bincount(np.concatenate([times, other]).astype(int), minlength=int(T))
print(f"merged: mean = {merged.mean():.3f}, variance = {merged.var():.3f} (both 5)")

# Thinning: 30% of calls are sales calls
sales = times[rng.random(times.size) < 0.3]
s_per_minute = np.bincount(sales.astype(int), minlength=int(T))
print(f"sales calls per minute: mean = {s_per_minute.mean():.3f}, variance = {s_per_minute.var():.3f} (both 1.2)")

calls per minute: mean = 3.996, variance = 4.000 (both 4)
corr(count in minute i, count in minute i+1) = -0.0083
P(no call in 30 s) = 0.1341 (e^-2 = 0.1353)
merged: mean = 4.997, variance = 5.039 (both 5)
sales calls per minute: mean = 1.198, variance = 1.205 (both 1.2)

Key takeaways

• A Poisson process has independent, stationary increments with Poisson$(\lambda t)$ counts; it is the continuous-time limit of the Bernoulli process.
• Merging independent Poisson processes adds the rates; random thinning multiplies the rate by the keep probability.
• Given the number of arrivals in an interval, their times are independent uniform points.
• Time-varying rates are handled by the non-homogeneous version, with mean count $\int \lambda(t),dt$.

🧪 Practice

1. Calls arrive at 4 per minute. Find the probability of no call in 30 seconds, and of exactly 2 calls in the first minute and 3 in the second.
2. Customers arrive at a shop at 20 per hour, and each buys something with probability 0.3, independently. What is the distribution of the number of buyers in 2 hours? Are the numbers of buyers and non-buyers independent?
3. Given that 5 arrivals occurred in an hour, what is the probability that all of them happened in the first 30 minutes?
4. Interview: How would you simulate a non-homogeneous Poisson process whose rate $\lambda(t)$ follows a daily cycle and never exceeds $\lambda_{\max}$? Hint: simulate a homogeneous process with rate $\lambda_{\max}$, then keep each point at time $t$ with probability $\lambda(t)/\lambda_{\max}$; explain why thinning makes this correct.

Interarrival Times

Why it matters. Counts per interval are one view of an arrival process; the other is the time between consecutive arrivals. Interarrival times determine how long a customer waits, how often a server is idle, how bursty traffic is, and whether the Poisson model fits real data at all. They also lead to the inspection paradox, one of the most common traps in reasoning about waiting times.

The result. In a Poisson process with rate $\lambda$, the interarrival times $T_1, T_2, \dots$ (the time to the first arrival, then between the first and the second, and so on) are independent Exponential$(\lambda)$. The first one is a one-line calculation:

$$P(T_1 > t) = P(\text{no arrival in } [0, t]) = P(N(t) = 0) = e^{-\lambda t}.$$

By the fresh-start property, the same holds after each arrival. Conversely, adding up independent Exponential$(\lambda)$ gaps generates a Poisson process. The time of the $n$-th arrival, $S_n = T_1 + \dots + T_n$, is Gamma$(n, \lambda)$, and

$${S_n \le t} = {N(t) \ge n},$$

which is the Gamma-Poisson identity of Section 5.2.

  N(t)
    4 |                                       +--------
    3 |                         +-------------+
    2 |              +----------+
    1 |      +-------+
    0 +------+-------+----------+-------------+-----------> t
             S_1     S_2        S_3           S_4
      |<-T1->|<--T2->|<---T3--->|<-----T4---->|
          interarrival times T_i ~ i.i.d. Exponential(lambda)
          arrival times     S_n ~ Gamma(n, lambda)

The inspection paradox. Arrive at a bus stop at a random time, with buses arriving as a Poisson process of rate $\lambda$. By memorylessness, your wait is Exponential$(\lambda)$ with mean $1/\lambda$, and the time since the last bus is also about $1/\lambda$. The gap you landed in therefore has mean length $2/\lambda$, twice the average gap. There is no contradiction: a random moment is more likely to fall into a long gap than a short one (length-biased sampling). The same effect makes your gym classes feel crowded (you are present when it is busy) and makes sampled hospital stays longer than the average stay.

Checking the Poisson assumption with gaps. The exponential distribution has coefficient of variation (SD/mean) exactly 1, which gives a quick diagnostic for real arrival data:

CV of interarrival times	Pattern	Example
$\approx 1$	random, consistent with Poisson	independent users hitting a website
$> 1$	bursty: clusters of arrivals separated by quiet spells	retries after an error, batch jobs, viral traffic
$< 1$	more regular than random	scheduled appointments, a production line

import numpy as np

rng = np.random.default_rng(25)
lam = 4.0                                              # arrivals per minute

# Build a Poisson process from exponential gaps
gaps = rng.exponential(1 / lam, size=1_000_000)
times = np.cumsum(gaps)
per_minute = np.bincount(times.astype(int))[:-1]       # drop the last, incomplete minute
print(f"counts per minute: mean = {per_minute.mean():.3f}, variance = {per_minute.var():.3f} (Poisson: both 4)")

s10 = gaps[:999_990].reshape(-1, 10).sum(axis=1)       # time of every 10th arrival, measured from the previous one
print(f"time for 10 arrivals: mean = {s10.mean():.4f} (n/lambda = 2.5), variance = {s10.var():.4f} (n/lambda^2 = 0.625)")

# Inspection paradox: the gap that contains a random moment is twice as long on average
moments = rng.uniform(0, times[-1], size=200_000)
idx = np.searchsorted(times, moments)                  # index of the first arrival after each moment
covering = times[idx] - np.where(idx > 0, times[idx - 1], 0.0)
wait = times[idx] - moments
print(f"average gap = {gaps.mean():.4f}, average gap containing a random moment = {covering.mean():.4f}, "
      f"average wait = {wait.mean():.4f}")

# The coefficient of variation of gaps as a diagnostic
bursty = np.where(rng.random(100_000) < 0.9, rng.exponential(0.05, 100_000), rng.exponential(2.0, 100_000))
regular = 0.25 + rng.uniform(-0.05, 0.05, 100_000)     # scheduled every 15 s, with jitter
for name, g in [("Poisson", gaps), ("bursty", bursty), ("scheduled", regular)]:
    print(f"CV of gaps, {name:9s}: {g.std() / g.mean():.3f}")

counts per minute: mean = 4.003, variance = 3.984 (Poisson: both 4)
time for 10 arrivals: mean = 2.4983 (n/lambda = 2.5), variance = 0.6249 (n/lambda^2 = 0.625)
average gap = 0.2498, average gap containing a random moment = 0.4994, average wait = 0.2494
CV of gaps, Poisson  : 1.000
CV of gaps, bursty   : 3.523
CV of gaps, scheduled: 0.116

Key takeaways

• In a Poisson process, the gaps between arrivals are independent Exponential$(\lambda)$, and the $n$-th arrival time is Gamma$(n, \lambda)$.
• "The $n$-th arrival by time $t$" equals "at least $n$ arrivals by time $t$", which links gamma and Poisson probabilities.
• The gap containing a random moment is length-biased, with mean $2/\lambda$ (inspection paradox).
• A coefficient of variation of the gaps near 1 is consistent with Poisson arrivals; above 1 means bursty, below 1 means regular.

🧪 Practice

1. Arrivals occur at rate 2 per hour. Find the probability that the 3rd arrival occurs within the first hour, using the Poisson count and then the gamma distribution.
2. A server receives requests as a Poisson process with rate 10 per second. What is the expected time until the 50th request, and its standard deviation?
3. You have already waited 5 minutes at a stop served by a Poisson bus process with mean gap 10 minutes. What is the distribution of your remaining wait? How would the answer change if buses ran on a fixed 10-minute schedule?
4. Interview: You have the timestamps of all API requests for one week. How would you test whether arrivals follow a Poisson process? Hint: check the gaps (exponential Q-Q plot, Section 5.3, and a CV near 1), the counts per interval (variance close to the mean), independence (autocorrelation of gaps or counts) and stationarity (the rate over the day); a daily cycle calls for a non-homogeneous model rather than a rejection of the Poisson idea.

Markov Chains and Transition Matrices

Why it matters. Bernoulli and Poisson processes have no memory at all. Many systems have a state that shapes what happens next: the weather today makes tomorrow's weather more predictable, a user on a product page is more likely to add to the cart than a user on the home page, a queue with 10 people will not be empty a second later. A Markov chain is the simplest model with memory: the future depends on the present state, but not on how the system got there. Markov chains power PageRank, text generation models, credit rating migrations, queueing models, and the MCMC algorithms of Bayesian statistics (Section 11.3).

Definition. A sequence $X_0, X_1, X_2, \dots$ with values in a set of states $S$ is a (time-homogeneous) Markov chain if, for all states and times,

$$P(X_{n+1} = j \mid X_n = i, X_{n-1}, \dots, X_0) = P(X_{n+1} = j \mid X_n = i) = P_{ij}.$$

This is the Markov property: given the present, the past is irrelevant for the future. The probabilities $P_{ij}$ form the transition matrix $P$, in which row $i$ is the distribution of the next state given the current state $i$. Every entry is non-negative and every row sums to 1 (a stochastic matrix).

Example: weather.

  transition matrix P                        state diagram, arrow by arrow

             Sunny  Cloudy  Rainy            Sunny  --0.2-->  Cloudy  --0.3-->  Rainy
  Sunny   [   0.7    0.2    0.1  ]           Sunny  <--0.3--  Cloudy  <--0.3--  Rainy
  Cloudy  [   0.3    0.4    0.3  ]           Sunny  --0.1---------------------> Rainy
  Rainy   [   0.2    0.3    0.5  ]           Sunny  <--0.2--------------------- Rainy
                                             stay: Sunny 0.7, Cloudy 0.4, Rainy 0.5

Multi-step transitions. The probability of going from $i$ to $j$ in $n$ steps is the $(i, j)$ entry of the matrix power $P^n$ (the Chapman-Kolmogorov equations): to go from $i$ to $j$ in two steps, sum over the intermediate state, $P^{(2)}{ij} = \sum_k P{ik}P_{kj}$, which is matrix multiplication. If the initial state has distribution given by the row vector $\pi_0$, the distribution at time $n$ is

$$\pi_n = \pi_0 P^n.$$

Classifying states.

Term	Meaning
Communicating states	each can be reached from the other
Irreducible chain	every state can be reached from every other state
Absorbing state	once entered, never left ($P_{ii} = 1$)
Transient state	the chain eventually leaves it for good
Recurrent state	the chain returns to it with probability 1
Period	gcd of the possible return times; aperiodic if 1

Absorbing chains. In a chain with absorbing states, such as a purchase funnel that ends in "purchase" or "leave", write the transition matrix among the transient states as $Q$ and from transient to absorbing states as $R$. The fundamental matrix $N = (I - Q)^{-1}$ gives the expected number of visits to each transient state; $NR$ gives the absorption probabilities, and $N\mathbf{1}$ the expected number of steps before absorption. These formulas are first-step analysis (Section 3.4) in matrix form.

Estimating a chain from data. Given an observed sequence of states, the maximum likelihood estimate of $P_{ij}$ is the fraction of transitions out of $i$ that went to $j$: count the transitions, then normalize each row.

import numpy as np

states = ["Sunny", "Cloudy", "Rainy"]
P = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])
assert np.allclose(P.sum(axis=1), 1)                   # each row is a probability distribution

P2 = P @ P
print(f"P(Rainy in 2 days | Sunny today) = {P2[0, 2]:.3f}")
pi0 = np.array([1.0, 0.0, 0.0])                       # it is sunny today
for n in [1, 2, 5, 10, 20]:
    print(f"day {n:2d}: " + "  ".join(f"{s} {p:.4f}" for s, p in zip(states, pi0 @ np.linalg.matrix_power(P, n))))

# Simulate 100,000 days, then estimate P back from the transition counts
rng = np.random.default_rng(26)
cum = P.cumsum(axis=1)
path = np.empty(100_000, dtype=int); path[0] = 0
u = rng.random(path.size)
for t in range(1, path.size):
    path[t] = np.searchsorted(cum[path[t - 1]], u[t])  # inverse transform on the current row
counts = np.zeros((3, 3))
np.add.at(counts, (path[:-1], path[1:]), 1)
print("estimated P:\n", np.round(counts / counts.sum(axis=1, keepdims=True), 3))

# Absorbing chain: a purchase funnel. Transient states: Home, Product, Cart. Absorbing: Purchase, Leave
Q = np.array([[0.1, 0.5, 0.0],                         # from Home
              [0.2, 0.1, 0.3],                         # from Product
              [0.0, 0.2, 0.1]])                        # from Cart
R = np.array([[0.0, 0.4],                              # columns: Purchase, Leave
              [0.0, 0.4],
              [0.5, 0.2]])
N = np.linalg.inv(np.eye(3) - Q)                       # fundamental matrix
for i, s in enumerate(["Home", "Product", "Cart"]):
    print(f"start at {s:7s}: P(purchase) = {(N @ R)[i, 0]:.3f}, expected page views before the end = {N.sum(axis=1)[i]:.2f}")

P(Rainy in 2 days | Sunny today) = 0.180
day  1: Sunny 0.7000  Cloudy 0.2000  Rainy 0.1000
day  2: Sunny 0.5700  Cloudy 0.2500  Rainy 0.1800
day  5: Sunny 0.4685  Cloudy 0.2794  Rainy 0.2521
day 10: Sunny 0.4568  Cloudy 0.2825  Rainy 0.2607
day 20: Sunny 0.4565  Cloudy 0.2826  Rainy 0.2609
estimated P:
 [[0.701 0.201 0.098]
 [0.3   0.397 0.303]
 [0.205 0.293 0.502]]
start at Home   : P(purchase) = 0.128, expected page views before the end = 2.31
start at Product: P(purchase) = 0.231, expected page views before the end = 2.15
start at Cart   : P(purchase) = 0.607, expected page views before the end = 1.59

By day 20 the distribution has stopped changing, and starting from a cloudy or rainy day leads to the same limit (change pi0 to check): the chain forgets its starting point. That limiting distribution is the subject of the next topic.

Key takeaways

• A Markov chain's next state depends only on its current state, through a transition matrix whose rows sum to 1.
• $n$-step transition probabilities are the entries of $P^n$, and the distribution at time $n$ is $\pi_0 P^n$.
• In absorbing chains, $N = (I - Q)^{-1}$ gives expected visits, absorption probabilities $NR$ and expected time to absorption.
• Estimate a transition matrix by counting observed transitions and normalizing each row.

🧪 Practice

1. For the weather chain, compute by hand the probability that it is sunny two days after a rainy day, and check it with the code.
2. Let $X_n$ be the outcome of the $n$-th toss of a fair coin and $Y_n = X_n + X_{n-1}$. Is $Y_n$ a Markov chain? Is the pair $(X_{n-1}, X_n)$ a Markov chain?
3. Modify the funnel so that 10% of users on the Cart page go back to Home. How do the purchase probabilities change?
4. Interview: What is the expected number of fair coin tosses until you see two heads in a row? Hint: define the states "start", "last toss was H" and "HH" (absorbing), write the transition matrix, and either solve the first-step equations for the expected time from each state or use the fundamental matrix.

Stationary Distributions

Why it matters. The weather example showed the distribution of a Markov chain settling down to a limit that does not depend on the starting state. That limit, the stationary distribution, describes the long-run behavior of the system: the fraction of days that are rainy, the fraction of time a server is overloaded, the importance of a web page in PageRank. MCMC algorithms (Section 11.3) work by building a Markov chain whose stationary distribution is the distribution we want to sample from.

Definition. A probability vector $\pi$ (non-negative entries summing to 1) is a stationary distribution of a chain with transition matrix $P$ if

$$\pi P = \pi, \qquad \text{that is,} \qquad \pi_j = \sum_i \pi_i P_{ij} \text{ for every state } j.$$

If the chain starts in $\pi$, it stays in $\pi$ forever: the probability flowing into each state equals the probability already there. In linear algebra terms, $\pi$ is a left eigenvector of $P$ with eigenvalue 1.

Existence, uniqueness and convergence. For a chain with finitely many states:

Condition	Consequence
Irreducible	a unique stationary distribution $\pi$ exists, with all $\pi_i > 0$
Irreducible and aperiodic	in addition, $\pi_0 P^n \to \pi$ from every starting distribution
Irreducible (ergodic theorem)	the long-run fraction of time spent in state $i$ is $\pi_i$, and the mean return time to $i$ is $1/\pi_i$

Periodicity matters for convergence: the chain that alternates deterministically between two states has stationary distribution $(1/2, 1/2)$, but $\pi_n$ oscillates forever and never converges, although the long-run time fractions are still $1/2$.

Computing $\pi$.

1. Linear system: solve $\pi(P - I) = 0$ together with $\sum_i \pi_i = 1$.
2. Eigenvector: take the eigenvector of $P^\top$ for eigenvalue 1 and normalize it to sum 1.
3. Power iteration: repeat $\pi \leftarrow \pi P$ from any start until it stops changing. This is how PageRank is computed for billions of pages.
4. Simulation: run the chain for a long time and record the fraction of time in each state.

Detailed balance. If $\pi_i P_{ij} = \pi_j P_{ji}$ for all pairs $i, j$, then $\pi$ is stationary (sum both sides over $i$). Such a chain is reversible: in equilibrium, the flow from $i$ to $j$ equals the flow back. Detailed balance is much easier to check than $\pi P = \pi$, and it is the design principle behind the Metropolis-Hastings algorithm.

PageRank. Model a web surfer who, at each step, follows a random link on the current page with probability $d$ (typically 0.85) or jumps to a uniformly random page with probability $1 - d$. A page's PageRank is its probability in the stationary distribution of this chain. The random jumps make the chain irreducible and aperiodic, so the stationary distribution exists, is unique, and power iteration converges to it.

import numpy as np

P = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])

# 1. Linear system: pi (P - I) = 0 plus the normalization sum(pi) = 1
A = np.vstack([(P - np.eye(3)).T, np.ones(3)])
pi_lin = np.linalg.lstsq(A, np.array([0, 0, 0, 1.0]), rcond=None)[0]

# 2. Left eigenvector for eigenvalue 1
w, v = np.linalg.eig(P.T)
pi_eig = np.real(v[:, np.argmin(np.abs(w - 1))])
pi_eig /= pi_eig.sum()

# 3. Power iteration from the uniform distribution
pi_pow = np.ones(3) / 3
for _ in range(100):
    pi_pow = pi_pow @ P

# 4. Long-run time fractions in a simulation of 200,000 days
rng = np.random.default_rng(27)
cum, state, visits, u = P.cumsum(axis=1), 0, np.zeros(3), rng.random(200_000)
for x in u:
    state = np.searchsorted(cum[state], x)
    visits[state] += 1
for name, pi in [("linear system", pi_lin), ("eigenvector", pi_eig), ("power iteration", pi_pow), ("simulation", visits / visits.sum())]:
    print(f"{name:16s} Sunny {pi[0]:.4f}  Cloudy {pi[1]:.4f}  Rainy {pi[2]:.4f}")
print(f"mean time between rainy days = 1/pi_Rainy = {1 / pi_lin[2]:.2f} days")

# PageRank on a tiny web: A -> B, A -> C, B -> C, C -> A, D -> C
links = {0: [1, 2], 1: [2], 2: [0], 3: [2]}
n, d = 4, 0.85
M = np.zeros((n, n))
for i, outs in links.items():
    M[i, outs] = 1 / len(outs)                       # follow a random outgoing link
G = d * M + (1 - d) / n                              # or teleport to a random page
rank = np.ones(n) / n
for _ in range(100):
    rank = rank @ G
print("PageRank:", "  ".join(f"{page} {r:.3f}" for page, r in zip("ABCD", rank)))

linear system    Sunny 0.4565  Cloudy 0.2826  Rainy 0.2609
eigenvector      Sunny 0.4565  Cloudy 0.2826  Rainy 0.2609
power iteration  Sunny 0.4565  Cloudy 0.2826  Rainy 0.2609
simulation       Sunny 0.4546  Cloudy 0.2836  Rainy 0.2618
mean time between rainy days = 1/pi_Rainy = 3.83 days
PageRank: A 0.373  B 0.196  C 0.394  D 0.038

The four methods agree. In the PageRank example, page C, which receives links from three pages, ranks highest, and page D, which no page links to, only gets the teleportation share $(1 - d)/n = 0.0375$.

Key takeaways

• A stationary distribution satisfies $\pi P = \pi$: a left eigenvector of $P$ with eigenvalue 1, normalized to sum 1.
• A finite irreducible chain has a unique $\pi$; if it is also aperiodic, the chain converges to $\pi$ from any start.
• $\pi_i$ is the long-run fraction of time in state $i$, and $1/\pi_i$ is the mean return time.
• Detailed balance $\pi_iP_{ij} = \pi_jP_{ji}$ implies stationarity; it is the basis of MCMC. PageRank is a stationary distribution.

🧪 Practice

1. For the two-state chain $P = \begin{pmatrix} 1 - a & a \ b & 1 - b \end{pmatrix}$, show that $\pi = \left(\frac{b}{a + b}, \frac{a}{a + b}\right)$. What happens when $a = b = 1$?
2. A random walk on an undirected graph moves from a node to a uniformly random neighbor. Show that $\pi_i \propto \deg(i)$ satisfies detailed balance.
3. Remove the random jumps from the PageRank example ($d = 1$). Does power iteration still converge? What happens if page C had no outgoing links?
4. Interview: Explain PageRank as a Markov chain to a non-specialist, and explain why the damping factor is needed. Hint: the random surfer; pages without outgoing links and groups of pages that link only to each other break irreducibility, and periodic structures break convergence; teleportation fixes both and also determines how fast power iteration converges.

Random Walks

Why it matters. Add up independent random steps and you get a random walk: the position of a diffusing particle, the fortune of a gambler, the log-price of a stock in the simplest model, the running difference between two A/B test variants. Random walks show how randomness accumulates. The typical distance grows like $\sqrt{n}$, not $n$, and the walk produces long excursions that look like trends although none exists. They are the discrete ancestors of Brownian motion and the engine of random-walk MCMC.

Definition. Let $X_1, X_2, \dots$ be independent and identically distributed steps. The random walk is

$$S_0 = 0, \qquad S_n = X_1 + X_2 + \dots + X_n.$$

The simple random walk has steps $\pm 1$ with probabilities $p$ and $q = 1 - p$; it is symmetric if $p = 1/2$.

Spread grows like $\sqrt n$. For the symmetric simple walk, $E[S_n] = 0$ and $\operatorname{Var}(S_n) = n$, so the typical distance from the start after $n$ steps is $\sqrt n$. By the central limit theorem, $S_n/\sqrt n \approx N(0, 1)$. After 10,000 fair coin tosses, the difference between heads and tails is typically around 100, even though the proportion of heads is very close to 1/2.

  S_n
   +20 |                                   .            upper envelope ~ +2 sqrt(n)
       |                       .'.  ..   .' '.  .
   +10 |            ..   .''..'   ''  '.'     ''.''.
       |   . .'.  .'  '.'                           '.
     0 +-.'-'---''-----------------------------------'-----> n
       |
   -20 |                                                    lower envelope ~ -2 sqrt(n)

Recurrence. Does the symmetric walk return to its starting point? On a line ($\mathbb{Z}$) and in the plane ($\mathbb{Z}^2$), it returns with probability 1, infinitely often. In three dimensions it returns with probability only about 0.34 (Pólya's theorem). In the words of Shizuo Kakutani: "A drunk man will find his way home, but a drunk bird may get lost forever." The expected time to return in one dimension, however, is infinite.

Gambler's ruin. A gambler starts with $k$ units and bets 1 unit at a time, winning each bet with probability $p$, and stops at 0 (ruin) or at $N$ (target). The probability of reaching $N$ is

$$P(\text{reach } N) = \begin{cases} \dfrac{k}{N} & p = \tfrac12, \[2ex] \dfrac{1 - (q/p)^k}{1 - (q/p)^N} & p \ne \tfrac12, \end{cases}$$

and for the fair game the expected duration is $k(N - k)$ bets. A small house edge has a large effect over many bets: with $p$ slightly below 1/2, the chance of doubling a large bankroll in unit bets is close to 0.

The arcsine law. In a fair game, one player is usually ahead for most of the time. The fraction of time the walk spends above 0 is not concentrated near 1/2: it follows the arcsine distribution, with density $\frac{1}{\pi\sqrt{x(1 - x)}}$ on $(0, 1)$, which is U-shaped. Long leads are the rule, not the exception, and that is why random noise is so easily mistaken for a trend.

import numpy as np

rng = np.random.default_rng(28)

# 10,000 symmetric walks of 1,000 steps
steps = rng.choice(np.array([-1, 1], dtype=np.int8), size=(10_000, 1000))
S = steps.cumsum(axis=1, dtype=np.int32)
for n in [10, 100, 1000]:
    print(f"n = {n:4d}: mean S_n = {S[:, n - 1].mean():+.2f}, SD = {S[:, n - 1].std():6.2f} (sqrt(n) = {np.sqrt(n):5.2f})")
print(f"fraction of walks that returned to 0 within 1,000 steps = {np.mean((S == 0).any(axis=1)):.3f}")

# Gambler's ruin: start with 10, stop at 0 or 20
def ruin_sim(p, k=10, N=20, players=20_000):
    x = np.full(players, k)
    active = np.ones(players, dtype=bool)
    bets = 0
    while active.any():
        step = np.where(rng.random(active.sum()) < p, 1, -1)
        x[active] += step
        bets += active.sum()
        active &= (x > 0) & (x < N)
    return np.mean(x == N), bets / players

for p in [0.5, 18 / 37]:                                  # fair game, and red/black at European roulette
    r = (1 - p) / p
    exact = 10 / 20 if p == 0.5 else (1 - r ** 10) / (1 - r ** 20)
    sim, duration = ruin_sim(p)
    print(f"p = {p:.4f}: P(reach 20) exact = {exact:.4f}, simulated = {sim:.4f}, mean number of bets = {duration:.1f}")

# Arcsine law: the fraction of time spent ahead is U-shaped, not centered at 1/2
frac = (S > 0).mean(axis=1)
hist, edges = np.histogram(frac, bins=10, range=(0, 1))
arcsine = np.diff(2 / np.pi * np.arcsin(np.sqrt(edges)))
print("time ahead   simulated   arcsine law")
for a, b, h, e in zip(edges[:-1], edges[1:], hist / len(frac), arcsine):
    print(f"{a:.1f} - {b:.1f}     {h:.3f}       {e:.3f}   {'#' * round(h * 100)}")

n =   10: mean S_n = -0.01, SD =   3.19 (sqrt(n) =  3.16)
n =  100: mean S_n = -0.10, SD =  10.10 (sqrt(n) = 10.00)
n = 1000: mean S_n = +0.30, SD =  31.60 (sqrt(n) = 31.62)
fraction of walks that returned to 0 within 1,000 steps = 0.973
p = 0.5000: P(reach 20) exact = 0.5000, simulated = 0.5003, mean number of bets = 100.5
p = 0.4865: P(reach 20) exact = 0.3680, simulated = 0.3643, mean number of bets = 97.1
time ahead   simulated   arcsine law
0.0 - 0.1     0.214       0.205   #####################
0.1 - 0.2     0.090       0.090   #########
0.2 - 0.3     0.074       0.074   #######
0.3 - 0.4     0.064       0.067   ######
0.4 - 0.5     0.062       0.064   ######
0.5 - 0.6     0.065       0.064   ######
0.6 - 0.7     0.066       0.067   #######
0.7 - 0.8     0.078       0.074   ########
0.8 - 0.9     0.095       0.090   #########
0.9 - 1.0     0.192       0.205   ###################

The extreme bins, where one side leads more than 90% of the time, are the most likely outcomes, while a near-even split of the lead is the least likely.

Key takeaways

• A random walk is a sum of independent steps; for the symmetric simple walk, $E[S_n] = 0$ and $\operatorname{SD}(S_n) = \sqrt n$.
• One- and two-dimensional symmetric walks always return to the start; three-dimensional walks may not.
• Gambler's ruin: $P(\text{reach } N) = k/N$ in a fair game, and a small negative edge makes reaching a distant target nearly impossible.
• The arcsine law: long leads are typical, so random walks produce convincing fake trends.

🧪 Practice

1. For a symmetric simple walk, find $E[S_{100}]$, $\operatorname{SD}(S_{100})$, and approximate $P(S_{100} \ge 20)$ with the normal distribution and a continuity correction (note that $S_{100}$ is always even).
2. A gambler with 50 units bets 1 unit at a time on a game with $p = 0.47$, aiming to reach 100. Compute the probability of reaching 100, and compare it with a single bet of all 50 units at the same $p$.
3. Prove that the expected duration of the fair gambler's ruin is $k(N - k)$ by first-step analysis: show that $D_k = 1 + \frac12 D_{k-1} + \frac12 D_{k+1}$ with $D_0 = D_N = 0$.
4. Interview: An A/B test is analyzed every day, and the team stops as soon as the p-value drops below 0.05. Even when the two variants are identical, they "find" a winner far more often than 5% of the time. Explain why. Hint: under the null, the test statistic over time behaves like a scaled random walk; checking it repeatedly gives it many chances to cross a fixed boundary, and such a walk eventually crosses any fixed level; think about sequential tests or a fixed sample size decided in advance (Section 9.6).

<a id="6-multivariate-probability"></a>