6. Multivariate Probability

Real questions almost always involve more than one random quantity at a time: height and weight, the return of two stocks, the number of visits and the number of purchases, the features and the label of a machine learning example. This chapter extends the single-variable tools of Chapter 4 to several random variables at once: joint, marginal and conditional distributions, covariance and correlation as measures of dependence, the distribution of sums and the most important multivariate models, and finally the inequalities and limit theorems (the laws of large numbers and the central limit theorem) that explain why averages of many random variables are so predictable. These results are the theoretical foundation of every method of estimation, testing and regression in the chapters that follow.

6.1 Joint Distributions

A joint distribution describes how several random variables behave together, which is information that the separate distributions of each variable do not contain.

Joint PMF and PDF

Why it matters. Knowing that a customer spends on average 40 dollars, and that a customer visits on average 3 times a month, says nothing about whether frequent visitors spend more. To answer questions that involve two variables at once, such as "what is the probability that a user visits often and spends little?", we need a single object that assigns probabilities to pairs of values. That object is the joint distribution. It is the multivariate version of the PMF and PDF of Sections 4.1 and 4.2.

Discrete case. For discrete random variables $X$ and $Y$ defined on the same sample space, the joint probability mass function is

$$p(x, y) = P(X = x, Y = y),$$

where the comma means "and". A joint PMF must satisfy $p(x, y) \ge 0$ and $\sum_x \sum_y p(x, y) = 1$, and the probability of any set $A$ of pairs is the sum of $p$ over the pairs in $A$:

$$P\big((X, Y) \in A\big) = \sum_{(x, y) \in A} p(x, y).$$

The natural picture of a joint PMF is a table, with one row per value of $X$, one column per value of $Y$, and probabilities in the cells.

Example: minimum and maximum of two dice. Roll two fair dice and let $X$ be the smaller value and $Y$ the larger. The pair $(X, Y) = (2, 5)$ happens for the rolls (2, 5) and (5, 2), so it has probability $2/36$; the pair $(3, 3)$ only for the roll (3, 3), so $1/36$; and $X > Y$ is impossible.

  p(x, y) x 36            Y = max
                   1    2    3    4    5    6
             1  [  1    2    2    2    2    2  ]
             2  [  0    1    2    2    2    2  ]
  X = min    3  [  0    0    1    2    2    2  ]
             4  [  0    0    0    1    2    2  ]
             5  [  0    0    0    0    1    2  ]
             6  [  0    0    0    0    0    1  ]

  P(Y - X >= 3) = cells on the 3rd diagonal and beyond = (2 * 6) / 36 = 1/3

Continuous case. Continuous random variables $X$ and $Y$ are jointly continuous if there is a function $f(x, y) \ge 0$, the joint probability density function, such that for every (reasonable) region $A$ of the plane

$$P\big((X, Y) \in A\big) = \iint_A f(x, y), dx, dy, \qquad \iint_{\mathbb{R}^2} f(x, y), dx, dy = 1.$$

Probability is now volume under the surface $z = f(x, y)$, just as in one dimension it was area under a curve (Section 4.2). As in one dimension, $f(x, y)$ is not a probability and can exceed 1; for a small rectangle around $(x, y)$,

$$P(x < X \le x + dx,; y < Y \le y + dy) \approx f(x, y), dx, dy.$$

A consequence is that any curve, including the line $X = Y$, has probability 0 for a jointly continuous pair.

Example. Let $f(x, y) = x + y$ on the unit square $0 \le x, y \le 1$ (and 0 elsewhere). It integrates to $\frac12 + \frac12 = 1$. The probability that $X + Y \le 1$ is the volume over the triangle below the anti-diagonal:

$$P(X + Y \le 1) = \int_0^1 \int_0^{1 - x} (x + y), dy, dx = \int_0^1 \frac{(1 - x)(1 + x)}{2}, dx = \frac13.$$

  y
  1 +-----------+          density f(x, y) = x + y grows toward the corner (1, 1)
    |\          |
    | \  x+y>1  |          region x + y <= 1 (lower-left triangle):
    |  \        |          half of the area, but only 1/3 of the probability,
    |   \       |          because the density is small near the origin
    |    \      |
    |x+y<=1\    |
    |       \   |
  0 +--------\--+-- x
    0           1

Analogy. Think of a joint PDF as a sheet of sand spread over a map, with one unit of sand in total. The probability of a region is the amount of sand lying over it. Thick sand means likely outcomes; the height of the sand at one point (the density) is not itself an amount of sand.

The code computes probabilities from the dice table and from the density, the latter both by numerical integration on a grid and by simulation. To simulate from $f(x, y) = x + y$ we use rejection sampling: propose a uniform point in the square and accept it with probability $f(x, y)/2$, where 2 is the maximum of the density.

import numpy as np

rng = np.random.default_rng(61)

# Discrete: joint PMF of (min, max) of two dice, built by enumerating the 36 rolls
p = np.zeros((6, 6))
for a in range(1, 7):
    for b in range(1, 7):
        p[min(a, b) - 1, max(a, b) - 1] += 1 / 36      # row = min - 1, column = max - 1
print(f"total probability = {p.sum():.4f}")
x, y = np.meshgrid(np.arange(1, 7), np.arange(1, 7), indexing="ij")
print(f"P(Y - X >= 3) = {p[y - x >= 3].sum():.4f}")    # sum the cells of the event
print(f"P(X = 2, Y = 5) = {p[1, 4]:.4f} (2/36 = {2 / 36:.4f})")

# Continuous: f(x, y) = x + y on the unit square
m = 2000
g = (np.arange(m) + 0.5) / m                          # midpoints of a fine grid
gx, gy = np.meshgrid(g, g, indexing="ij")
f = gx + gy
cell = (1 / m) ** 2                                   # area of one grid cell
print(f"integral of f       = {(f * cell).sum():.4f}")
print(f"P(X + Y <= 1), grid = {(f * cell)[gx + gy <= 1].sum():.4f} (exact 1/3)")

# Rejection sampling: accept a uniform proposal with probability f / max(f)
u = rng.random((400_000, 2))
keep = rng.random(len(u)) < u.sum(axis=1) / 2
s = u[keep]
print(f"accepted {len(s)} points; P(X + Y <= 1), simulation = {np.mean(s.sum(axis=1) <= 1):.4f}")
print(f"P(X > 0.5, Y > 0.5), simulation = {np.mean((s > 0.5).all(axis=1)):.4f} (exact 0.375)")

total probability = 1.0000
P(Y - X >= 3) = 0.3333
P(X = 2, Y = 5) = 0.0556 (2/36 = 0.0556)
integral of f       = 1.0000
P(X + Y <= 1), grid = 0.3336 (exact 1/3)
accepted 199788 points; P(X + Y <= 1), simulation = 0.3316
P(X > 0.5, Y > 0.5), simulation = 0.3750 (exact 0.375)

The last line is worth a second look: the upper-right quarter of the square holds 37.5% of the probability on 25% of the area, because that is where the density is largest. The grid estimate of 1/3 is slightly off because the grid cells that straddle the line $x + y = 1$ are counted whole.

Key takeaways

The joint PMF $p(x, y) = P(X = x, Y = y)$ is a table of probabilities that are non-negative and sum to 1; probabilities of events are sums over cells.
The joint PDF $f(x, y)$ gives probabilities as volumes, $P((X, Y) \in A) = \iint_A f$; it is a density, not a probability, and can exceed 1.
For jointly continuous variables, lines and curves (such as $X = Y$) have probability 0.
Simulation by rejection sampling turns any bounded joint density into samples.

🧪 Practice

From the dice table above, compute $P(X = Y)$, $P(Y = 2X)$ and $P(X + Y = 7)$.

Let $f(x, y) = c,xy$ on $0 \le x \le 1$, $0 \le y \le 2$. Find $c$, then compute $P(X \le 1/2, Y \ge 1)$ and $P(Y \le X)$.

Let $f(x, y) = 2$ on the triangle $0 < x < y < 1$. Sketch the support and compute $P(Y > 2X)$.

Interview: You have a table of 10 million (user, item) click events. How would you estimate the joint distribution of (number of clicks per user, number of distinct categories per user), and what would you do about the cells that have very few users? Hint: the empirical joint PMF is the normalized 2D histogram (a crosstab); sparse cells can be pooled into coarser bins, smoothed, or modeled parametrically; think about what goes wrong when you estimate a probability from 3 observations.

Joint CDF

Why it matters. The CDF was the one description that worked for every single random variable, discrete, continuous or mixed (Sections 4.1 and 4.2). The joint CDF plays the same role for several variables: it always exists, it determines the joint distribution completely, and it is what multivariate tools such as copulas (Section 6.3) are built on.

Definition. The joint cumulative distribution function of $X$ and $Y$ is

$$F(x, y) = P(X \le x,; Y \le y).$$

It is the probability of the lower-left quadrant with corner $(x, y)$:

          y
          ^
          |
  ////////+ (x, y)         F(x, y) = probability of the shaded quadrant
  ////////|                         { X <= x and Y <= y }
  ////////|
  --------+---------> x
  ////////|

Properties.

Property	Statement
Monotone	non-decreasing in each argument
Limits at $-\infty$	$F(x, y) \to 0$ if $x \to -\infty$ or $y \to -\infty$
Limit at $+\infty$	$F(x, y) \to 1$ if $x \to \infty$ and $y \to \infty$
Marginals	$F_X(x) = \lim_{y \to \infty} F(x, y)$, $F_Y(y) = \lim_{x \to \infty} F(x, y)$
Right-continuous	in each argument

Rectangle probabilities. Being non-decreasing in each argument is not enough for a function of two variables to be a joint CDF: every rectangle must also get non-negative probability. By inclusion-exclusion (Section 3.2), the probability of the rectangle $(a, b] \times (c, d]$ is

$$P(a < X \le b,; c < Y \le d) = F(b, d) - F(a, d) - F(b, c) + F(a, c).$$

  d +-----------+----------+       F(b, d) counts everything below and left of (b, d);
    |           |##########|       subtract the strip left of x = a  : F(a, d)
    |           |##########|       subtract the strip below y = c    : F(b, c)
  c +-----------+----------+       the corner below-left of (a, c) was subtracted twice,
    |           |          |       so add it back                    : F(a, c)
    +-----------+----------+
                a          b

From the CDF to the density. For jointly continuous variables, the density is the mixed partial derivative of the CDF, and the CDF is the double integral of the density:

$$f(x, y) = \frac{\partial^2 F(x, y)}{\partial x, \partial y}, \qquad F(x, y) = \int_{-\infty}^{x} \int_{-\infty}^{y} f(s, t), dt, ds.$$

Example. For $f(x, y) = x + y$ on the unit square, for $0 \le x, y \le 1$:

$$F(x, y) = \int_0^x \int_0^y (s + t), dt, ds = \frac{x^2 y + x y^2}{2} = \frac{xy(x + y)}{2}.$$

Check: $\partial F/\partial y = (x^2 + 2xy)/2$, and then $\partial/\partial x$ gives $x + y$. Setting $y = 1$ gives the marginal CDF of $X$, $F_X(x) = (x^2 + x)/2$.

The empirical joint CDF. From a sample of pairs $(x_i, y_i)$, the empirical joint CDF $\hat F(x, y)$ is the fraction of points with $x_i \le x$ and $y_i \le y$. It is the bivariate version of the ECDF, and it converges to $F$ as the sample grows.

import numpy as np

rng = np.random.default_rng(62)

def F(x, y):
    """Joint CDF of f(x, y) = x + y on the unit square (arguments clipped to [0, 1])."""
    x, y = np.clip(x, 0, 1), np.clip(y, 0, 1)
    return x * y * (x + y) / 2

# Rectangle (0.25, 0.5] x (0.5, 1] by inclusion-exclusion
a, b, c, d = 0.25, 0.5, 0.5, 1.0
rect = F(b, d) - F(a, d) - F(b, c) + F(a, c)
print(f"P(0.25 < X <= 0.5, 0.5 < Y <= 1) = {rect:.6f} (direct integral: 0.140625)")

# Mixed partial derivative of F recovers the density (finite differences)
h, x0, y0 = 1e-4, 0.3, 0.6
mixed = (F(x0 + h, y0 + h) - F(x0 + h, y0) - F(x0, y0 + h) + F(x0, y0)) / h**2
print(f"d2F/dxdy at (0.3, 0.6) = {mixed:.4f} (f = {x0 + y0:.4f})")

# Empirical joint CDF from rejection samples
u = rng.random((200_000, 2))
s = u[rng.random(len(u)) < u.sum(axis=1) / 2]
for (x, y) in [(0.5, 0.5), (0.3, 0.9), (1.0, 0.4)]:
    emp = np.mean((s[:, 0] <= x) & (s[:, 1] <= y))
    print(f"F({x}, {y}): exact {F(x, y):.4f}, empirical {emp:.4f}")

P(0.25 < X <= 0.5, 0.5 < Y <= 1) = 0.140625 (direct integral: 0.140625)
d2F/dxdy at (0.3, 0.6) = 0.9001 (f = 0.9000)
F(0.5, 0.5): exact 0.1250, empirical 0.1256
F(0.3, 0.9): exact 0.1620, empirical 0.1630
F(1.0, 0.4): exact 0.2800, empirical 0.2806

Key takeaways

$F(x, y) = P(X \le x, Y \le y)$ exists for every pair of random variables and determines the joint distribution.
Rectangle probabilities come from four CDF values by inclusion-exclusion: $F(b, d) - F(a, d) - F(b, c) + F(a, c)$.
Marginal CDFs are limits of the joint CDF: let the other variable go to $+\infty$.
For jointly continuous variables, $f = \partial^2 F / \partial x, \partial y$.

🧪 Practice

Using the dice table of the previous topic, compute $F(2, 4) = P(\min \le 2, \max \le 4)$.

$X$ and $Y$ have joint CDF $F(x, y) = (1 - e^{-x})(1 - e^{-2y})$ for $x, y \ge 0$. Find the joint PDF, the marginal CDFs and $P(1 < X \le 2, Y > 1)$.

Show that $G(x, y) = 1$ if $x + y \ge 0$ and 0 otherwise is non-decreasing in each argument and has the right limits, yet is not a joint CDF. Hint: compute the rectangle probability for $(-1, 1] \times (-1, 1]$.

Marginal Distributions

Why it matters. Often we have the joint distribution but want to talk about one variable alone: from a joint model of (age, income), what is the distribution of income? Getting a single-variable distribution out of a joint one is called marginalization, and it is everywhere: summing a crosstab over its rows, integrating out nuisance parameters in Bayesian statistics (Chapter 11), computing the evidence in Bayes' theorem. The name comes from the practice of writing the row and column totals in the margins of the table.

Discrete case. Sum the joint PMF over all values of the other variable. This is the law of total probability (Section 3.4), applied to the partition ${Y = y}$:

$$p_X(x) = \sum_y p(x, y), \qquad p_Y(y) = \sum_x p(x, y).$$

Continuous case. Integrate out the other variable:

$$f_X(x) = \int_{-\infty}^{\infty} f(x, y), dy, \qquad f_Y(y) = \int_{-\infty}^{\infty} f(x, y), dx.$$

Example: dice. The row sums of the min/max table give the distribution of the minimum, the column sums that of the maximum:

  x 36                     Y = max                        row sums
                  1    2    3    4    5    6           p_X(x) x 36   (min)
             1 [  1    2    2    2    2    2  ]            11
             2 [  0    1    2    2    2    2  ]             9
  X = min    3 [  0    0    1    2    2    2  ]             7
             4 [  0    0    0    1    2    2  ]             5
             5 [  0    0    0    0    1    2  ]             3
             6 [  0    0    0    0    0    1  ]             1
  column sums     1    3    5    7    9   11                36
  p_Y(y) x 36  (max)

So $p_X(k) = (13 - 2k)/36$ and $p_Y(k) = (2k - 1)/36$: small values are likely for the minimum and large values for the maximum.

Example: continuous. For $f(x, y) = x + y$ on the unit square, $f_X(x) = \int_0^1 (x + y), dy = x + \frac12$ for $0 \le x \le 1$, and by symmetry $f_Y(y) = y + \frac12$.

Watch the limits of integration. When the support is not a rectangle, the limits of the inner integral depend on the outer variable. For $f(x, y) = 2$ on the triangle $0 < x < y < 1$:

$$f_X(x) = \int_x^1 2, dy = 2(1 - x), \qquad f_Y(y) = \int_0^y 2, dx = 2y.$$

Drawing the support first and reading off "for a fixed $x$, which $y$ are allowed?" avoids most mistakes.

Marginals do not determine the joint distribution. Many different joint distributions share the same marginals. Two fair coins can be independent, always equal, or always opposite: in all three cases each coin is fair, but the joint tables differ completely. The marginals describe each variable; how they relate is extra information. Copulas (Section 6.3) make this separation explicit.

  independent             always equal            always opposite
       0     1                 0     1                 0     1
  0 [ .25   .25 ] .5      0 [ .50   .00 ] .5      0 [ .00   .50 ] .5
  1 [ .25   .25 ] .5      1 [ .00   .50 ] .5      1 [ .50   .00 ] .5
      .5    .5                .5    .5                .5    .5        same margins, different joints

import numpy as np

rng = np.random.default_rng(63)

# Discrete marginals: row and column sums of the joint table
p = np.zeros((6, 6))
for a in range(1, 7):
    for b in range(1, 7):
        p[min(a, b) - 1, max(a, b) - 1] += 1 / 36
k = np.arange(1, 7)
print("p_min(k) x 36:", np.round(p.sum(axis=1) * 36).astype(int), " formula:", 13 - 2 * k)
print("p_max(k) x 36:", np.round(p.sum(axis=0) * 36).astype(int), " formula:", 2 * k - 1)
print(f"E[min] = {(k * p.sum(axis=1)).sum():.4f}, E[max] = {(k * p.sum(axis=0)).sum():.4f}")

# Continuous marginal of f(x, y) = 2 on 0 < x < y < 1: sample uniformly on the triangle
u = np.sort(rng.random((200_000, 2)), axis=1)         # sorting two uniforms gives a uniform point with x < y
hist, edges = np.histogram(u[:, 0], bins=5, range=(0, 1), density=True)
mid = (edges[:-1] + edges[1:]) / 2
for m_, h in zip(mid, hist):
    print(f"x = {m_:.1f}: histogram density {h:.3f}, f_X(x) = 2(1 - x) = {2 * (1 - m_):.3f}")

p_min(k) x 36: [11  9  7  5  3  1]  formula: [11  9  7  5  3  1]
p_max(k) x 36: [ 1  3  5  7  9 11]  formula: [ 1  3  5  7  9 11]
E[min] = 2.5278, E[max] = 4.4722
x = 0.1: histogram density 1.801, f_X(x) = 2(1 - x) = 1.800
x = 0.3: histogram density 1.401, f_X(x) = 2(1 - x) = 1.400
x = 0.5: histogram density 0.996, f_X(x) = 2(1 - x) = 1.000
x = 0.7: histogram density 0.600, f_X(x) = 2(1 - x) = 0.600
x = 0.9: histogram density 0.202, f_X(x) = 2(1 - x) = 0.200

Key takeaways

Marginal PMF: sum the joint PMF over the other variable; marginal PDF: integrate the joint PDF over the other variable.
With a non-rectangular support, the limits of integration depend on the variable that is kept; draw the support first.
Marginals are determined by the joint distribution, but not the other way around: the dependence structure is extra information.

🧪 Practice

A joint PMF is given by $p(0, 0) = 0.1$, $p(0, 1) = 0.3$, $p(1, 0) = 0.2$, $p(1, 1) = 0.4$. Find both marginals and $E[X]$, $E[Y]$.

For $f(x, y) = 6x$ on $0 < x < y < 1$, find $f_X$ and $f_Y$ and check that each integrates to 1.

Construct two different joint PMFs on ${0, 1, 2} \times {0, 1, 2}$ whose marginals are both uniform on ${0, 1, 2}$.

Independence of Random Variables

Why it matters. Independence is the assumption behind almost every formula in elementary statistics: the variance of a sample mean, the likelihood of an i.i.d. sample, the binomial distribution. It is also a strong assumption that fails in many datasets (repeated measurements on the same user, time series, clustered samples), and when it fails, standard errors and p-values become wrong. Knowing exactly what it means, and how to check it in a model, is essential.

Definition. Random variables $X$ and $Y$ are independent if every event about $X$ is independent of every event about $Y$ (Section 3.3):

$$P(X \in A,; Y \in B) = P(X \in A), P(Y \in B) \quad \text{for all sets } A, B.$$

Equivalent and more practical forms:

Setting	Independence holds if and only if
Any	$F(x, y) = F_X(x), F_Y(y)$ for all $x, y$
Discrete	$p(x, y) = p_X(x), p_Y(y)$ for all $x, y$
Continuous	$f(x, y) = f_X(x), f_Y(y)$ for all $x, y$ (except on a set of probability 0)

In words: the joint table is the product of its margins; every row of the table is proportional to every other row. Mutual independence of $X_1, \dots, X_n$ means that the joint distribution factors into all $n$ marginals, which is stronger than independence of every pair (Section 3.3).

The factorization criterion. A joint density is a product of independent marginals if and only if it can be written as

$$f(x, y) = g(x), h(y) \quad \text{on a support that is a product set (a rectangle).}$$

There is no need to find the normalizing constants of $g$ and $h$. Two quick checks:

Density	Independent?	Why
$4xy$ on $[0, 1]^2$	yes	factors as $(2x)(2y)$ on a square
$x + y$ on $[0, 1]^2$	no	a sum cannot be written as $g(x)h(y)$
$2$ on $0 < x < y < 1$	no	constant density, but the support is a triangle: knowing $X = 0.9$ forces $Y > 0.9$
$e^{-x-y}$ on $x, y > 0$	yes	$e^{-x} e^{-y}$ on a quadrant

The triangle case is the classic trap: a support that is not a rectangle always means dependence, because the possible values of one variable depend on the other.

Functions of independent variables are independent. If $X$ and $Y$ are independent, then so are $g(X)$ and $h(Y)$ for any functions $g, h$. For example, $X^2$ and $e^Y$ are independent.

Independence simplifies everything. For independent $X$ and $Y$:

$E[XY] = E[X],E[Y]$, and more generally $E[g(X)h(Y)] = E[g(X)],E[h(Y)]$;
$\operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y)$;
the MGF of the sum is the product of the MGFs, $M_{X+Y}(t) = M_X(t),M_Y(t)$ (Section 4.3);
the joint likelihood of independent observations is the product of the individual likelihoods.

Checking independence in a table. In data, compare each cell of the joint table with the product of its margins; the chi-square test of independence (Section 9.3) formalizes this comparison.

import numpy as np

rng = np.random.default_rng(64)

def independence_gap(p):
    """Largest |p(x, y) - p_X(x) p_Y(y)| over all cells: 0 means independent."""
    outer = np.outer(p.sum(axis=1), p.sum(axis=0))    # product of the margins
    return np.abs(p - outer).max()

# Two dice: (first, second) versus (min, max)
first_second = np.full((6, 6), 1 / 36)
min_max = np.zeros((6, 6))
for a in range(1, 7):
    for b in range(1, 7):
        min_max[min(a, b) - 1, max(a, b) - 1] += 1 / 36
print(f"(first die, second die): max gap = {independence_gap(first_second):.4f}")
print(f"(min, max)             : max gap = {independence_gap(min_max):.4f}")

# Continuous: f = 4xy (independent) versus f = 2 on the triangle 0 < x < y < 1 (dependent)
a = np.sqrt(rng.random((200_000, 2)))                 # inverse CDF of 2x gives independent X, Y with f = 4xy
t = np.sort(rng.random((200_000, 2)), axis=1)         # uniform on the triangle
for name, s in [("4xy on square", a), ("2 on triangle", t)]:
    joint = np.mean((s[:, 0] > 0.5) & (s[:, 1] < 0.6))
    prod = np.mean(s[:, 0] > 0.5) * np.mean(s[:, 1] < 0.6)
    print(f"{name}: P(X > 0.5, Y < 0.6) = {joint:.4f}, product of marginals = {prod:.4f}, "
          f"E[XY] - E[X]E[Y] = {np.mean(s[:, 0] * s[:, 1]) - s[:, 0].mean() * s[:, 1].mean():+.4f}")

(first die, second die): max gap = 0.0000
(min, max)             : max gap = 0.0378
4xy on square: P(X > 0.5, Y < 0.6) = 0.2713, product of marginals = 0.2710, E[XY] - E[X]E[Y] = -0.0002
2 on triangle: P(X > 0.5, Y < 0.6) = 0.0100, product of marginals = 0.0900, E[XY] - E[X]E[Y] = +0.0277

For the triangle, the joint probability is much smaller than the product of the marginals: a large $X$ makes a small $Y$ nearly impossible.

Key takeaways

$X$ and $Y$ are independent when the joint PMF, PDF or CDF factors into the product of the marginals, for all values.
Factorization criterion: $f(x, y) = g(x)h(y)$ on a rectangular support; a non-rectangular support always means dependence.
Functions of independent variables are independent; expectations of products factor; variances of sums add.
In data, compare the joint table with the outer product of its margins.

🧪 Practice

Are $X$ and $Y$ independent if $p(0, 0) = 0.12$, $p(0, 1) = 0.28$, $p(1, 0) = 0.18$, $p(1, 1) = 0.42$? Justify with the margins.

Decide whether each density describes independent variables: $f(x, y) = 6e^{-2x - 3y}$ on $x, y > 0$; $f(x, y) = 8xy$ on $0 < x < y < 1$; $f(x, y) = \frac{3}{2}(x^2 + y^2)$ on $[0, 1]^2$.

Roll two dice. Show that the events "the sum is 7" and "the first die shows 4" are independent, but the random variables $S = $ sum and $X_1 = $ first die are not.

Interview: A dataset contains 50 measurements from each of 200 patients. A colleague computes a standard error as $s/\sqrt{10{,}000}$. What is wrong, and in which direction is the error? Hint: measurements from the same patient are dependent; positive dependence makes the effective sample size smaller than 10,000, so the standard error is too small and p-values too optimistic.

Joint Distributions of Transformations

Why it matters. Many useful random variables are functions of several others: the ratio of two measurements, polar coordinates of a random point, the sum and difference of two scores. Finding the joint distribution of new variables built from old ones is the multivariate version of the change of variables formula (Section 4.4). It explains, among other things, how the Box-Muller algorithm turns uniform random numbers into normal ones.

The Jacobian method. Let $(X, Y)$ have joint density $f_{X,Y}$, and let $(U, V) = g(X, Y)$ be a one-to-one, differentiable transformation with inverse $x = x(u, v)$, $y = y(u, v)$. Then

$$f_{U,V}(u, v) = f_{X,Y}\big(x(u, v), y(u, v)\big), \lvert J \rvert, \qquad J = \det \begin{pmatrix} \partial x / \partial u & \partial x / \partial v \ \partial y / \partial u & \partial y / \partial v \end{pmatrix}.$$

The absolute Jacobian $\lvert J \rvert$ is the local factor by which the inverse map stretches areas: a small rectangle of area $du,dv$ in the $(u, v)$ plane corresponds to a small parallelogram of area $\lvert J \rvert, du, dv$ in the $(x, y)$ plane, and probability must be conserved.

     (u, v) plane                          (x, y) plane

     +----+                                    ____
     |    | area du dv     inverse map        /   /    area |J| du dv
     +----+               ------------->     /___/     same probability, so
                                                       f_UV du dv = f_XY |J| du dv

Recipe.

Write the inverse transformation $x(u, v)$, $y(u, v)$.
Compute the Jacobian determinant of the inverse and take its absolute value.
Substitute into $f_{X,Y}$ and multiply by $\lvert J \rvert$.
Find the support of $(U, V)$: the image of the support of $(X, Y)$.

Equivalently, $\lvert J \rvert = 1 / \lvert J_g \rvert$, where $J_g$ is the Jacobian of the forward map, evaluated at the inverse point.

Example: sum and proportion of exponentials. Let $X, Y$ be independent Exponential(1), and $U = X + Y$, $V = X/(X + Y)$. The inverse is $x = uv$, $y = u(1 - v)$, so

$$J = \det \begin{pmatrix} v & u \ 1 - v & -u \end{pmatrix} = -uv - u(1 - v) = -u, \qquad f_{U,V}(u, v) = e^{-uv} e^{-u(1 - v)} \cdot u = u e^{-u}, \quad u > 0,; 0 < v < 1.$$

The density factors into $u e^{-u}$ (a Gamma(2, 1) density) times 1 (a Uniform(0, 1) density) on a rectangle. So the total $U$ is Gamma(2, 1), the share $V$ is uniform, and, surprisingly, the total and the share are independent. The same calculation with Gamma($a$, 1) and Gamma($b$, 1) variables shows that $X/(X + Y)$ is Beta($a$, $b$), the link between the gamma and beta families of Section 5.2.

Example: Box-Muller. Let $U_1, U_2$ be independent Uniform(0, 1) and

$$R = \sqrt{-2 \ln U_1}, \quad \Theta = 2\pi U_2, \qquad X = R\cos\Theta, \quad Y = R\sin\Theta.$$

Working through the Jacobian of polar coordinates (which is $r$) shows that $X$ and $Y$ are independent standard normals: $f(x, y) = \frac{1}{2\pi} e^{-(x^2 + y^2)/2}$. The radius carries the "how far" information and the angle the "which direction", and the normal density depends only on the distance from the origin.

One function of two variables. To find the distribution of a single $U = g(X, Y)$, add an auxiliary variable (often $V = X$), find the joint density of $(U, V)$, and integrate out $V$. With $U = X + Y$, $V = X$ this gives the convolution formula of Section 6.4. Alternatively, use the CDF method: $F_U(u) = P(g(X, Y) \le u)$ is the probability of a region of the plane.

import numpy as np

rng = np.random.default_rng(65)
n = 500_000

# Sum and proportion of two independent Exponential(1) variables
x, y = rng.exponential(1, n), rng.exponential(1, n)
u, v = x + y, x / (x + y)
print(f"U = X + Y : mean {u.mean():.3f} (Gamma(2,1): 2), variance {u.var():.3f} (2)")
print(f"V = X/(X+Y): mean {v.mean():.3f} (uniform: 0.5), variance {v.var():.4f} (1/12 = {1 / 12:.4f})")
print(f"P(U < 1, V < 0.3) = {np.mean((u < 1) & (v < 0.3)):.4f}, "
      f"P(U < 1) P(V < 0.3) = {np.mean(u < 1) * np.mean(v < 0.3):.4f}")   # equal: independent

# Box-Muller: two uniforms -> two independent standard normals
u1, u2 = rng.random(n), rng.random(n)
r, theta = np.sqrt(-2 * np.log(u1)), 2 * np.pi * u2
z1, z2 = r * np.cos(theta), r * np.sin(theta)
print(f"Z1: mean {z1.mean():+.4f}, SD {z1.std():.4f}, P(Z1 <= 1.96) = {np.mean(z1 <= 1.96):.4f} (0.9750)")
print(f"corr(Z1, Z2) = {np.corrcoef(z1, z2)[0, 1]:+.4f}, corr(Z1^2, Z2^2) = {np.corrcoef(z1**2, z2**2)[0, 1]:+.4f}")

U = X + Y : mean 2.003 (Gamma(2,1): 2), variance 2.009 (2)
V = X/(X+Y): mean 0.500 (uniform: 0.5), variance 0.0833 (1/12 = 0.0833)
P(U < 1, V < 0.3) = 0.0790, P(U < 1) P(V < 0.3) = 0.0791
Z1: mean +0.0011, SD 0.9989, P(Z1 <= 1.96) = 0.9747 (0.9750)
corr(Z1, Z2) = +0.0001, corr(Z1^2, Z2^2) = -0.0029

The correlation of the squares is checked as well, because uncorrelated variables can still be dependent (Section 6.3); a zero there too is consistent with independence.

Key takeaways

For a one-to-one transformation, $f_{U,V}(u, v) = f_{X,Y}(x(u, v), y(u, v)),\lvert J \rvert$, with $J$ the Jacobian determinant of the inverse map.
$\lvert J \rvert$ is the area-scaling factor that keeps total probability equal to 1; do not forget the absolute value or the new support.
For a single function $U = g(X, Y)$, add an auxiliary variable and integrate it out, or use the CDF method.
Classic results: the sum and proportion of i.i.d. exponentials are independent (Gamma and Uniform); Box-Muller turns two uniforms into two independent normals.

🧪 Practice

$X, Y$ are independent Uniform(0, 1). Find the joint density of $U = X + Y$, $V = X - Y$, and sketch its support.

$X, Y$ are independent standard normals. Use the transformation $U = X + Y$, $V = X - Y$ to show that $U$ and $V$ are independent $N(0, 2)$.

$X, Y$ are independent Exponential(1). Find the density of $W = X/Y$ by adding the auxiliary variable $V = Y$.

Interview: Explain why Box-Muller produces normal samples, and why many libraries prefer other methods (such as the ziggurat algorithm) in practice. Hint: write the standard bivariate normal in polar coordinates and identify the distributions of $R^2$ and $\Theta$; then think about the cost of $\log$, $\sqrt{}$, $\sin$ and $\cos$ per sample.

6.2 Conditional Distributions

Conditioning on the value of one variable turns a joint distribution into a family of single-variable distributions, and the expectations and variances of these distributions are the language of prediction and regression.

Conditional PMF and PDF

Why it matters. Prediction is conditioning. "What is the distribution of a house's price, given that it has 3 bedrooms?" "How long will this user stay, given that they have already stayed 10 minutes?" The answer to each is a conditional distribution: the distribution of one variable once the value of another is known. It is the variable-level version of conditional probability (Section 3.3), and every regression and classification model is, at heart, a model of a conditional distribution.

Discrete case. For a value $x$ with $p_X(x) > 0$, the conditional PMF of $Y$ given $X = x$ is

$$p_{Y \mid X}(y \mid x) = P(Y = y \mid X = x) = \frac{p(x, y)}{p_X(x)}.$$

Graphically: take the row of the joint table for $X = x$ and divide it by its total so that it sums to 1. Conditioning is slicing and renormalizing.

Example: dice. Given that the minimum of two dice is 3, the row $X = 3$ of the min/max table is $(0, 0, 1, 2, 2, 2)/36$ with total $7/36$, so

  row X = 3 of the joint table (x 36):   y = 1   2   3   4   5   6
                                              0   0   1   2   2   2     total 7

  conditional PMF p(y | X = 3):          y = 1   2   3   4   5   6
                                              0   0  1/7 2/7 2/7 2/7    total 1

The maximum is 3 only if both dice show 3, which is half as likely as each of the other possible values.

Continuous case. For $x$ with $f_X(x) > 0$, the conditional PDF of $Y$ given $X = x$ is

$$f_{Y \mid X}(y \mid x) = \frac{f(x, y)}{f_X(x)}.$$

Geometrically, it is the vertical slice of the joint density surface at $X = x$, rescaled to have area 1. For $f(x, y) = x + y$ on the unit square, $f_X(x) = x + \frac12$, so

$$f_{Y \mid X}(y \mid x) = \frac{x + y}{x + \frac12}, \qquad 0 \le y \le 1.$$

For $x = 0$ the slice is $2y$, strongly tilted toward large $y$; for $x = 1$ it is $(1 + y)/1.5$, much flatter.

Conditioning on an event of probability 0. For a continuous $X$, the event ${X = x}$ has probability 0, so $P(Y \le y \mid X = x)$ cannot be defined by the formula $P(A \cap B)/P(B)$. The conditional density is instead understood as a limit of conditioning on the thin strip $x \le X \le x + \varepsilon$. This works for nice densities, but it depends on how the strip shrinks: conditioning on "the point lies on the equator" gives different answers depending on whether the equator is seen as a limit of latitude bands or of longitude wedges (the Borel-Kolmogorov paradox). In practice, conditional densities are always defined relative to a specific variable, and the formula above is the right one.

Multiplication rule and Bayes' rule for densities. Rearranging the definition gives the tools of Section 3.4 for random variables:

$$f(x, y) = f_X(x), f_{Y \mid X}(y \mid x), \qquad f_{X \mid Y}(x \mid y) = \frac{f_{Y \mid X}(y \mid x), f_X(x)}{\int f_{Y \mid X}(y \mid t), f_X(t), dt}.$$

The first equation is how joint models are usually built: choose a distribution for $X$, then a distribution for $Y$ given each value of $X$. The second is Bayesian inference in one line (Chapter 11).

Independence in terms of conditionals. $X$ and $Y$ are independent exactly when the conditional distribution of $Y$ given $X = x$ is the same for every $x$, and equal to the marginal of $Y$: learning $X$ changes nothing about $Y$.

import numpy as np

rng = np.random.default_rng(66)

# Discrete: condition the (min, max) table on min = 3 by slicing and renormalizing
p = np.zeros((6, 6))
for a in range(1, 7):
    for b in range(1, 7):
        p[min(a, b) - 1, max(a, b) - 1] += 1 / 36
row = p[2]                                            # the row X = 3
print("p(y | min = 3):", np.round(row / row.sum(), 4))

# Continuous: f(x, y) = x + y; approximate conditioning on X = 0.2 by a thin strip
u = rng.random((4_000_000, 2))
s = u[rng.random(len(u)) < u.sum(axis=1) / 2]         # rejection sampling (Section 6.1)
x0, eps = 0.2, 0.005
strip = s[np.abs(s[:, 0] - x0) < eps, 1]               # Y values of points with X close to 0.2
print(f"points in the strip: {len(strip)}")
hist, edges = np.histogram(strip, bins=5, range=(0, 1), density=True)
for m, h in zip((edges[:-1] + edges[1:]) / 2, hist):
    print(f"y = {m:.1f}: strip density {h:.3f}, f(y | x = 0.2) = {(x0 + m) / (x0 + 0.5):.3f}")

p(y | min = 3): [0.     0.     0.1429 0.2857 0.2857 0.2857]
points in the strip: 13918
y = 0.1: strip density 0.418, f(y | x = 0.2) = 0.429
y = 0.3: strip density 0.710, f(y | x = 0.2) = 0.714
y = 0.5: strip density 1.018, f(y | x = 0.2) = 1.000
y = 0.7: strip density 1.302, f(y | x = 0.2) = 1.286
y = 0.9: strip density 1.552, f(y | x = 0.2) = 1.571

Key takeaways

The conditional PMF or PDF is the joint distribution sliced at $X = x$ and renormalized: $f_{Y \mid X}(y \mid x) = f(x, y)/f_X(x)$.
For continuous $X$, conditioning on $X = x$ is a limit of conditioning on thin strips; the formula is well defined, but only relative to the variable being conditioned on.
Joint = marginal $\times$ conditional is the standard way to build multivariate models; Bayes' rule for densities follows from it.
Independence means every conditional distribution of $Y$ equals its marginal.

🧪 Practice

From the dice table, find the conditional PMF of the minimum given that the maximum is 4.

Let $f(x, y) = 2$ on $0 < x < y < 1$. Find $f_{Y \mid X}(y \mid x)$ and identify it as a named distribution.

$X \sim$ Uniform(0, 1), and given $X = x$, $Y \sim$ Binomial($n$, $x$). Write the joint distribution, and use Bayes' rule to find the conditional density of $X$ given $Y = k$ (up to a constant). Which family does it belong to?

Interview: A classifier outputs $P(\text{fraud} \mid \text{features})$. Explain why this is a conditional distribution, and what can go wrong if the fraud rate in production differs from the one in the training data. Hint: write $P(\text{fraud} \mid x) \propto P(x \mid \text{fraud}), P(\text{fraud})$ and ask which factor changed.

Conditional Expectation

Why it matters. If you must predict $Y$ with a single number once you know $X = x$, the best choice (in mean squared error) is the mean of the conditional distribution. This conditional expectation $E[Y \mid X = x]$, viewed as a function of $x$, is called the regression function, and estimating it is what regression and most supervised learning methods do. It is also a powerful calculation tool, as the next topics show.

Definition. The conditional expectation of $Y$ given $X = x$ is the mean of the conditional distribution:

$$E[Y \mid X = x] = \sum_y y, p_{Y \mid X}(y \mid x) \quad \text{or} \quad \int y, f_{Y \mid X}(y \mid x), dy.$$

Example: dice. Given that the minimum is 3, the maximum is 3, 4, 5, 6 with probabilities $1/7, 2/7, 2/7, 2/7$, so $E[\max \mid \min = 3] = (3 + 8 + 10 + 12)/7 = 33/7 \approx 4.71$.

Example: continuous. For $f(x, y) = x + y$:

$$E[Y \mid X = x] = \int_0^1 y, \frac{x + y}{x + \frac12}, dy = \frac{x/2 + 1/3}{x + 1/2} = \frac{3x + 2}{6x + 3},$$

which decreases from $2/3$ at $x = 0$ to $5/9$ at $x = 1$.

$E[Y \mid X]$ is a random variable. Write $h(x) = E[Y \mid X = x]$. Plugging in the random $X$ gives the random variable $E[Y \mid X] = h(X)$: a function of $X$, whose randomness comes from not yet knowing $X$. For the dice, $E[\max \mid \min]$ takes the value $33/7$ with probability $P(\min = 3) = 7/36$, and so on. This distinction between the number $E[Y \mid X = x]$ and the random variable $E[Y \mid X]$ is essential for the law of iterated expectation.

Properties.

Property	Formula
Linearity	$E[aY_1 + bY_2 \mid X] = a,E[Y_1 \mid X] + b,E[Y_2 \mid X]$
Taking out what is known	$E[g(X),Y \mid X] = g(X),E[Y \mid X]$
Independence	if $X, Y$ are independent, $E[Y \mid X] = E[Y]$
Function of $X$	$E[g(X) \mid X] = g(X)$
Tower (next topic)	$E\big[E[Y \mid X]\big] = E[Y]$

"Taking out what is known" says that once $X$ is known, any function of $X$ behaves like a constant.

The best predictor. Among all functions $g$, the conditional expectation minimizes the mean squared prediction error:

$$E\big[(Y - g(X))^2\big] \ge E\big[(Y - E[Y \mid X])^2\big] \quad \text{for every function } g.$$

The reason is the decomposition $E[(Y - g(X))^2] = E[(Y - h(X))^2] + E[(h(X) - g(X))^2]$, where the cross term vanishes because the error $Y - h(X)$ has conditional mean 0 given $X$. The error $Y - E[Y \mid X]$ is uncorrelated with every function of $X$: there is no signal in $X$ left to exploit. Linear regression finds the best linear predictor, which coincides with $E[Y \mid X]$ only when the regression function happens to be linear (as it is for the bivariate normal, Section 6.4).

Analogy. $E[Y \mid X]$ is a smart average: instead of averaging $Y$ over everyone, average it over the people who share your value of $X$. In data, that is literally df.groupby("x")["y"].mean().

import numpy as np

rng = np.random.default_rng(67)
n = 1_000_000

# X ~ Uniform(-2, 2), Y = X^2 + noise: the regression function is E[Y | X = x] = x^2
x = rng.uniform(-2, 2, n)
y = x**2 + rng.normal(0, 0.5, n)

# Estimate E[Y | X] by averaging Y within narrow bins of X (a group-by mean)
bins = np.linspace(-2, 2, 9)
idx = np.digitize(x, bins) - 1
for b in range(0, 8, 2):
    mid = (bins[b] + bins[b + 1]) / 2
    print(f"x in [{bins[b]:+.1f}, {bins[b + 1]:+.1f}): mean of Y = {y[idx == b].mean():.3f}, "
          f"E[X^2 | bin] = {(bins[b]**2 + bins[b] * bins[b + 1] + bins[b + 1]**2) / 3:.3f}")

# Mean squared error of three predictors
slope, intercept = np.polyfit(x, y, 1)                # best linear predictor (least squares)
preds = {"constant E[Y]": np.full(n, y.mean()),
         "best linear": intercept + slope * x,
         "E[Y | X] = X^2": x**2}
for name, pred in preds.items():
    print(f"{name:15s}: MSE = {np.mean((y - pred) ** 2):.4f}")
print(f"fitted slope = {slope:+.4f}: X and Y are uncorrelated, yet X predicts Y very well")

x in [-2.0, -1.5): mean of Y = 3.086, E[X^2 | bin] = 3.083
x in [-1.0, -0.5): mean of Y = 0.583, E[X^2 | bin] = 0.583
x in [+0.0, +0.5): mean of Y = 0.085, E[X^2 | bin] = 0.083
x in [+1.0, +1.5): mean of Y = 1.582, E[X^2 | bin] = 1.583
constant E[Y]  : MSE = 1.6736
best linear    : MSE = 1.6736
E[Y | X] = X^2 : MSE = 0.2502
fitted slope = -0.0021: X and Y are uncorrelated, yet X predicts Y very well

The best linear predictor is no better than a constant here, while the conditional mean leaves only the irreducible noise variance $0.5^2 = 0.25$.

Key takeaways

$E[Y \mid X = x]$ is the mean of the conditional distribution; as a function of $x$ it is the regression function.
$E[Y \mid X]$ is a random variable, a function of $X$.
Known functions of $X$ factor out: $E[g(X)Y \mid X] = g(X)E[Y \mid X]$.
$E[Y \mid X]$ is the best predictor of $Y$ from $X$ in mean squared error; the best linear predictor can be much worse when the relationship is non-linear.

🧪 Practice

Compute $E[\max \mid \min = k]$ for $k = 1, \dots, 6$ from the dice table.

For $f(x, y) = 2$ on $0 < x < y < 1$, find $E[Y \mid X = x]$ and $E[X \mid Y = y]$.

Show that if $E[Y \mid X] = a + bX$, then $b = \operatorname{Cov}(X, Y)/\operatorname{Var}(X)$. Hint: multiply by $X - E[X]$ and take expectations, using the tower property.

Interview: A gradient-boosted model trained with squared error loss outputs predictions for house prices. What quantity is it estimating, and what would change if it were trained with absolute error loss instead? Hint: squared error targets the conditional mean; absolute error targets the conditional median.

Law of Iterated Expectation

Why it matters. Many expectations are hard to compute directly but easy to compute in two stages: first condition on something that simplifies the problem, then average over it. The total revenue of a day with a random number of customers, the expected number of steps of a randomized algorithm, the mean of a mixture distribution: all of them yield to the law of iterated expectation, also called the tower property or the law of total expectation.

Statement. For any random variables $X$ and $Y$ (with $E\lvert Y \rvert < \infty$):

$$E[Y] = E\big[E[Y \mid X]\big].$$

In the discrete case this is just the law of total probability applied to means:

$$E[Y] = \sum_x E[Y \mid X = x], P(X = x).$$

The overall average is the weighted average of the group averages, weighted by the group sizes.

                       E[Y | X = x1] = 10   (weight P(X = x1) = 0.2)
                     /
  E[Y] = 0.2 * 10   ---- E[Y | X = x2] = 20   (weight 0.5)
       + 0.5 * 20    \
       + 0.3 * 40 =   E[Y | X = x3] = 40   (weight 0.3)
       = 24

Example: random sums (Wald's identity). A shop has $N$ customers per day, and customer $i$ spends $X_i$, where the $X_i$ are i.i.d. with mean $\mu$ and independent of $N$. The daily revenue is $S = X_1 + \dots + X_N$, a sum with a random number of terms. Conditioning on $N$:

$$E[S \mid N = n] = n\mu \quad\Longrightarrow\quad E[S \mid N] = N\mu \quad\Longrightarrow\quad E[S] = E[N\mu] = E[N],\mu.$$

Example: hatching eggs. A hen lays $N \sim$ Poisson($\lambda$) eggs and each hatches independently with probability $p$. Given $N$, the number hatched $H$ is Binomial($N$, $p$), so $E[H \mid N] = Np$ and $E[H] = \lambda p$. (In fact $H \sim$ Poisson($\lambda p$), the Poisson thinning property of Section 5.1.)

Example: first-step analysis. The expected number of fair coin tosses until the first head, $m = E[T]$, satisfies $m = \frac12 \cdot 1 + \frac12 (1 + m)$ by conditioning on the first toss, so $m = 2$. All first-step analyses (Section 3.4, Section 5.5) are the law of iterated expectation in disguise.

A more general version. For nested information, the smaller information set wins: $E\big[E[Y \mid X, Z] \mid X\big] = E[Y \mid X]$. Averaging a detailed forecast over the details you do not have gives the coarse forecast. This is the property that makes martingales and sequential prediction work.

Pitfall: Simpson's paradox. Because the overall mean is a weighted average of group means, the overall comparison between two treatments can reverse the comparison within every group when the group weights differ between treatments (Section 2.5). The law of iterated expectation is the arithmetic behind the paradox.

import numpy as np

rng = np.random.default_rng(68)
days = 200_000

# Random sum: N ~ Poisson(20) customers, each spends Exponential with mean 30
N = rng.poisson(20, days)
spend = rng.exponential(30, N.sum())                  # all purchases of all days, in one array
day_of = np.repeat(np.arange(days), N)                # which day each purchase belongs to
S = np.bincount(day_of, weights=spend, minlength=days)
print(f"E[S]: simulated {S.mean():.2f}, Wald E[N] * mu = {20 * 30}")

# Hatching eggs: H | N ~ Binomial(N, p), N ~ Poisson(lambda)
lam, p = 12, 0.3
H = rng.binomial(rng.poisson(lam, days), p)
print(f"E[H]: simulated {H.mean():.3f}, lambda * p = {lam * p:.3f}; Var(H) = {H.var():.3f} (Poisson: equal to the mean)")

# Tower property on the dice: E[max] = sum over k of E[max | min = k] P(min = k)
p_min = (13 - 2 * np.arange(1, 7)) / 36
cond = [(k + 2 * sum(range(k + 1, 7))) / (13 - 2 * k) for k in range(1, 7)]   # E[max | min = k]
print("E[max | min = k]:", np.round(cond, 3))
print(f"weighted average = {np.dot(cond, p_min):.4f}, E[max] = 161/36 = {161 / 36:.4f}")

E[S]: simulated 600.76, Wald E[N] * mu = 600
E[H]: simulated 3.608, lambda * p = 3.600; Var(H) = 3.591 (Poisson: equal to the mean)
E[max | min = k]: [3.727 4.222 4.714 5.2   5.667 6.   ]
weighted average = 4.4722, E[max] = 161/36 = 4.4722

Key takeaways

$E[Y] = E[E[Y \mid X]]$: average the conditional means, weighted by the distribution of $X$.
Choose $X$ so that $E[Y \mid X]$ is easy: the number of terms of a random sum, the first step of a process, the component of a mixture.
Wald's identity: for a random sum with $N$ independent of the i.i.d. terms, $E[S] = E[N],E[X]$.
The weighting by group size is what makes Simpson's paradox possible.

🧪 Practice

30% of customers are on the premium plan and spend on average 80 dollars; the others spend on average 25 dollars. What is the average spend of a random customer?

Roll a die, then toss a fair coin as many times as the die shows. Find the expected number of heads.

A miner is trapped behind three doors: one leads out in 3 hours, the others lead back to the start after 5 and 7 hours. If the miner picks a door uniformly at random each time (forgetting past choices), what is the expected time to get out?

Interview: The average order value rose in every country last quarter, yet the global average order value fell. Is this possible, and how would you explain it to a stakeholder? Hint: write the global average as a weighted average of country averages and think about a shift in the weights.

Law of Total Variance

Why it matters. The law of iterated expectation splits a mean into pieces; the law of total variance does the same for the spread. It says that the total variability of $Y$ is the sum of the variability within groups and the variability between group means. This decomposition is the idea behind analysis of variance (Section 9.4), the $R^2$ of regression (Chapter 10), random effects models, and the variance of random sums used in insurance and queueing.

Statement. Using the conditional variance $\operatorname{Var}(Y \mid X)$, the variance of the conditional distribution of $Y$ given $X$ (studied in detail in the next topic),

$$\operatorname{Var}(Y) = \underbrace{E\big[\operatorname{Var}(Y \mid X)\big]}{\text{within groups: unexplained}} + \underbrace{\operatorname{Var}\big(E[Y \mid X]\big)}{\text{between groups: explained by } X}.$$

It is sometimes called Eve's law (EVVE: expectation of variance plus variance of expectation), in parallel with Adam's law for the tower property.

Proof sketch. Write $h(X) = E[Y \mid X]$ and $Y - \mu = (Y - h(X)) + (h(X) - \mu)$. Squaring and taking expectations, the cross term vanishes because $Y - h(X)$ has conditional mean 0 given $X$; the first square averages to $E[\operatorname{Var}(Y \mid X)]$ and the second to $\operatorname{Var}(h(X))$.

  Y values in three groups

  group A:  |--x--x-x--x-|                    within-group spread  -> E[Var(Y | X)]
  group B:              |--x-x--x--x-|
  group C:                         |-x--x--x-x--|
            ^ mean A       ^ mean B       ^ mean C
            |<-------- spread of the group means -------->|   -> Var(E[Y | X])

  total spread of all points = within + between

Example: a mixture of two groups. Adult heights: half the population has mean 163 cm and the other half mean 177 cm, both with SD 7 cm. Then $E[\operatorname{Var}(Y \mid G)] = 49$ and $\operatorname{Var}(E[Y \mid G]) = 7^2 = 49$ (the group means are $170 \pm 7$), so $\operatorname{Var}(Y) = 98$ and $\operatorname{SD}(Y) \approx 9.9$ cm. The fraction of the variance explained by the group is $49/98 = 50%$.

Example: random sums. For $S = X_1 + \dots + X_N$ with i.i.d. $X_i$ (mean $\mu$, variance $\sigma^2$) independent of $N$: given $N$, $E[S \mid N] = N\mu$ and $\operatorname{Var}(S \mid N) = N\sigma^2$. Hence

$$\operatorname{Var}(S) = E[N],\sigma^2 + \operatorname{Var}(N),\mu^2.$$

The first term is the uncertainty in the amounts; the second, the uncertainty in how many there are. For Poisson $N$ (compound Poisson), $\operatorname{Var}(N) = E[N] = \lambda$ and $\operatorname{Var}(S) = \lambda E[X^2]$.

Explained variance. The ratio $\operatorname{Var}(E[Y \mid X]) / \operatorname{Var}(Y)$ is the fraction of the variance of $Y$ that $X$ can explain with the best possible predictor. It is the population version of $R^2$, sometimes called the correlation ratio $\eta^2$. Because $\operatorname{Var}(E[Y \mid X]) \le \operatorname{Var}(Y)$, conditioning never increases variance on average.

import numpy as np

rng = np.random.default_rng(69)
n = 1_000_000

# Mixture of two groups: within + between
g = rng.integers(0, 2, n)                             # group label, each with probability 1/2
h = rng.normal(np.where(g == 0, 163, 177), 7)
within = np.mean([h[g == k].var() for k in (0, 1)])   # groups have equal size, so a plain mean is the weighted mean
between = np.var(np.where(g == 0, h[g == 0].mean(), h[g == 1].mean()))
print(f"Var(Y) = {h.var():.2f}; within = {within:.2f}, between = {between:.2f}, sum = {within + between:.2f}")
print(f"explained fraction = {between / h.var():.3f}")

# Random sum: N ~ Poisson(20) customers, each spends Exponential(mean 30): Var(X) = 900, E[X^2] = 1800
days = 200_000
N = rng.poisson(20, days)
S = np.bincount(np.repeat(np.arange(days), N), weights=rng.exponential(30, N.sum()), minlength=days)
print(f"Var(S): simulated {S.var():.0f}, E[N] Var(X) + Var(N) mu^2 = {20 * 900 + 20 * 900}")

# Same mean number of customers but overdispersed: N ~ Negative binomial with variance 120
N2 = rng.negative_binomial(4, 4 / 24, days)          # mean 20, variance 20 + 20^2/4 = 120
S2 = np.bincount(np.repeat(np.arange(days), N2), weights=rng.exponential(30, N2.sum()), minlength=days)
print(f"overdispersed N: Var(N) = {N2.var():.1f}, SD(S) = {S2.std():.1f} vs {S.std():.1f} with Poisson N")

Var(Y) = 98.05; within = 48.98, between = 49.07, sum = 98.05
explained fraction = 0.500
Var(S): simulated 36333, E[N] Var(X) + Var(N) mu^2 = 36000
overdispersed N: Var(N) = 120.3, SD(S) = 355.7 vs 190.6 with Poisson N

Customer counts that vary more from day to day (the second case) make the daily revenue noticeably less predictable, even though the expected revenue is unchanged.

Key takeaways

$\operatorname{Var}(Y) = E[\operatorname{Var}(Y \mid X)] + \operatorname{Var}(E[Y \mid X])$: within-group plus between-group variance.
The between part is the variance explained by $X$; its share of the total is the population $R^2$ (correlation ratio).
Random sums: $\operatorname{Var}(S) = E[N]\sigma^2 + \operatorname{Var}(N)\mu^2$; uncertainty in the count adds to uncertainty in the amounts.
Conditioning reduces variance on average: $E[\operatorname{Var}(Y \mid X)] \le \operatorname{Var}(Y)$.

🧪 Practice

Exam scores in three classes of equal size have means 60, 70, 80 and SDs 10, 12, 8. Compute the overall variance and the fraction explained by the class.

$X \sim$ Uniform(0, 1) and, given $X = x$, $Y \sim$ Binomial(10, $x$). Find $E[Y]$ and $\operatorname{Var}(Y)$, and compare $\operatorname{Var}(Y)$ with the binomial variance $10 \cdot 0.5 \cdot 0.5$.

An insurer receives $N \sim$ Poisson(100) claims per year with i.i.d. sizes of mean 5,000 and SD 8,000 dollars. Find the mean and SD of the total claims.

Interview: A model has an $R^2$ of 0.3 on a large test set. Using the law of total variance, explain what that means about the remaining 70%, and whether a better model could push it to 1. Hint: separate $E[\operatorname{Var}(Y \mid X)]$, which no model using only $X$ can remove, from the part of $\operatorname{Var}(E[Y \mid X])$ the model failed to capture.

Conditional Variance

Why it matters. A prediction without a measure of its uncertainty is incomplete. The conditional variance tells how much $Y$ still varies once $X$ is known, and that amount often depends on $X$: the spending of high-income customers is more variable than that of low-income ones, and errors in forecasting far in the future are larger than in the near term. Modeling this heteroscedasticity is essential for honest prediction intervals, for weighted least squares, and for volatility models in finance.

Definition. The conditional variance of $Y$ given $X = x$ is the variance of the conditional distribution:

$$\operatorname{Var}(Y \mid X = x) = E\big[(Y - E[Y \mid X = x])^2 \mid X = x\big] = E[Y^2 \mid X = x] - \big(E[Y \mid X = x]\big)^2.$$

As with the conditional mean, plugging in $X$ gives a random variable $\operatorname{Var}(Y \mid X)$, a function of $X$.

Properties.

Property	Formula
Shift by a function of $X$	$\operatorname{Var}(Y + g(X) \mid X) = \operatorname{Var}(Y \mid X)$
Scale by a function of $X$	$\operatorname{Var}(g(X),Y \mid X) = g(X)^2 \operatorname{Var}(Y \mid X)$
Independence	if $X, Y$ are independent, $\operatorname{Var}(Y \mid X) = \operatorname{Var}(Y)$
Function of $X$	$\operatorname{Var}(g(X) \mid X) = 0$
Average	$E[\operatorname{Var}(Y \mid X)] \le \operatorname{Var}(Y)$ (law of total variance)

Note that the conditional variance for a specific $x$ can be larger than $\operatorname{Var}(Y)$; it is only on average that conditioning reduces variance.

Example: dice. Given $\min = 3$, the maximum takes values 3, 4, 5, 6 with probabilities $1/7, 2/7, 2/7, 2/7$. $E[\max^2 \mid \min = 3] = (9 + 32 + 50 + 72)/7 = 163/7$, so $\operatorname{Var}(\max \mid \min = 3) = 163/7 - (33/7)^2 = 52/49 \approx 1.06$.

Homoscedastic and heteroscedastic models.

Model	Conditional mean	Conditional variance
Bivariate normal (Section 6.4)	linear in $x$	constant: $\sigma_Y^2(1 - \rho^2)$
Classical linear regression	$\beta_0 + \beta_1 x$	constant $\sigma^2$ (homoscedastic)
Poisson count with mean $\lambda(x)$	$\lambda(x)$	$\lambda(x)$: grows with the mean
Multiplicative noise $Y = m(x)(1 + \varepsilon)$	$m(x)$	$m(x)^2 \sigma_\varepsilon^2$: SD proportional to the mean

  homoscedastic                         heteroscedastic (fan shape)
  y |        .  . .  .  .               y |                    .  .
    |    . .  .  . .  .  .                |               .  .   . .
    |  .  . . .  .  . .                   |          . . . .  .  .
    | . . .  . .  .                       |     . .. . .  .
    |. .  . .                             | ......  .
    +--------------------- x              +--------------------- x
     same spread at every x                 spread grows with x

Prediction intervals. If $Y$ given $X = x$ is approximately normal, a 95% prediction interval is

$$E[Y \mid X = x] \pm 1.96 \sqrt{\operatorname{Var}(Y \mid X = x)}.$$

Using a single, average variance for every $x$ in a heteroscedastic situation gives intervals that are too wide where the noise is small and too narrow where it is large, even though their average coverage may look fine.

import numpy as np

rng = np.random.default_rng(70)
n = 1_000_000

# Heteroscedastic model: Y = 2 + 0.5 X + noise with SD 0.2 X
x = rng.uniform(1, 10, n)
y = 2 + 0.5 * x + rng.normal(0, 0.2 * x)

# Estimate Var(Y | X) in bins of X
edges = np.arange(1, 11, 3)
for lo, hi in zip(edges[:-1], edges[1:]):
    m = (x >= lo) & (x < hi)
    resid = y[m] - (2 + 0.5 * x[m])                   # subtract the known conditional mean
    print(f"x in [{lo}, {hi}): Var(Y | X) ~ {resid.var():.3f}, model E[(0.2X)^2 | bin] = {0.04 * (lo*lo + lo*hi + hi*hi) / 3:.3f}")

# Coverage of 95% intervals: one global SD versus the conditional SD
resid = y - (2 + 0.5 * x)
global_sd = resid.std()
inside_global = np.abs(resid) <= 1.96 * global_sd
inside_cond = np.abs(resid) <= 1.96 * 0.2 * x
for lo, hi in [(1, 4), (4, 7), (7, 10)]:
    m = (x >= lo) & (x < hi)
    print(f"x in [{lo}, {hi}): coverage with global SD {inside_global[m].mean():.3f}, with conditional SD {inside_cond[m].mean():.3f}")
print(f"overall coverage: global {inside_global.mean():.3f}, conditional {inside_cond.mean():.3f}")

x in [1, 4): Var(Y | X) ~ 0.279, model E[(0.2X)^2 | bin] = 0.280
x in [4, 7): Var(Y | X) ~ 1.244, model E[(0.2X)^2 | bin] = 1.240
x in [7, 10): Var(Y | X) ~ 2.935, model E[(0.2X)^2 | bin] = 2.920
x in [1, 4): coverage with global SD 1.000, with conditional SD 0.950
x in [4, 7): coverage with global SD 0.964, with conditional SD 0.949
x in [7, 10): coverage with global SD 0.839, with conditional SD 0.949
overall coverage: global 0.934, conditional 0.950

The global interval covers 93% overall, not far from the target, but it does so by covering everything for small $x$ and only 84% for large $x$, often exactly where the stakes are highest. The conditional interval covers 95% everywhere.

Key takeaways

$\operatorname{Var}(Y \mid X = x) = E[Y^2 \mid X = x] - (E[Y \mid X = x])^2$ is the spread of $Y$ left once $X$ is known.
Functions of $X$ behave like constants: they shift without changing the conditional variance and scale it by their square.
Conditional variance that changes with $x$ (heteroscedasticity) is common; it calls for prediction intervals whose width depends on $x$.
Only on average does conditioning reduce variance; for a particular $x$ the conditional variance can exceed $\operatorname{Var}(Y)$.

🧪 Practice

For $f(x, y) = x + y$ on the unit square, compute $\operatorname{Var}(Y \mid X = 0)$ and $\operatorname{Var}(Y \mid X = 1)$.

Given $X = x$, $Y \sim$ Poisson($3x$), with $X \sim$ Uniform(0, 2). Find $\operatorname{Var}(Y \mid X)$, $E[\operatorname{Var}(Y \mid X)]$ and $\operatorname{Var}(Y)$.

Find a joint distribution and a value $x$ for which $\operatorname{Var}(Y \mid X = x) > \operatorname{Var}(Y)$.

Interview: Residuals of a regression of sales on store size fan out as size grows. What are two ways to handle this, and what does each change? Hint: think about modeling $\log(\text{sales})$, and about weighted least squares with weights proportional to $1/\operatorname{Var}(Y \mid X)$, plus heteroscedasticity-robust standard errors.

6.3 Dependence and Covariance

This section measures how much and in which way two random variables move together, from the single number of covariance to the full dependence structure captured by a copula.

Covariance of Random Variables

Why it matters. The full joint distribution is often too much information; we want a single number that says whether two variables tend to move in the same direction. Covariance is that number. It appears in the variance of every sum (portfolio risk, the error of an average), in the slope of linear regression, and in principal component analysis. Section 2.5 introduced the sample covariance of a dataset; here we define its population counterpart for random variables.

Definition. The covariance of $X$ and $Y$, with means $\mu_X$ and $\mu_Y$, is

$$\operatorname{Cov}(X, Y) = E\big[(X - \mu_X)(Y - \mu_Y)\big] = E[XY] - E[X],E[Y].$$

The second form, the "mean of the product minus the product of the means", is usually easier to compute.

Reading the sign. The product $(X - \mu_X)(Y - \mu_Y)$ is positive when both variables are on the same side of their means and negative when they are on opposite sides. Covariance is the average of this product:

              Y
              ^
   (-)  X below, Y above   |   (+)  both above
                           |
   ------------------------+------------------------> X
                       (mu_X, mu_Y)
   (+)  both below         |   (-)  X above, Y below

   mass mostly in (+) quadrants -> Cov > 0;  mostly in (-) -> Cov < 0;  balanced -> Cov = 0

Properties. Covariance is bilinear, like a dot product:

Property	Formula
Symmetry	$\operatorname{Cov}(X, Y) = \operatorname{Cov}(Y, X)$
Variance as a special case	$\operatorname{Cov}(X, X) = \operatorname{Var}(X)$
Constants	$\operatorname{Cov}(X, c) = 0$
Scaling and shifting	$\operatorname{Cov}(aX + b, cY + d) = ac,\operatorname{Cov}(X, Y)$
Additivity	$\operatorname{Cov}(X + Y, Z) = \operatorname{Cov}(X, Z) + \operatorname{Cov}(Y, Z)$
Sums	$\operatorname{Cov}\big(\sum_i a_i X_i, \sum_j b_j Y_j\big) = \sum_i \sum_j a_i b_j \operatorname{Cov}(X_i, Y_j)$
Independence	$X, Y$ independent $\Rightarrow \operatorname{Cov}(X, Y) = 0$ (the converse is false, see below)

The most used consequence is the variance of a sum:

$$\operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y) + 2\operatorname{Cov}(X, Y).$$

Units. Covariance has the units of $X$ times the units of $Y$ (for example, cm times kg), so its size is hard to interpret: measuring height in mm instead of cm multiplies it by 10. Its sign is meaningful; its magnitude is not, until it is standardized into a correlation (next topic).

Example: dice. For the minimum $X$ and maximum $Y$ of two dice, $XY$ is the product of the two dice, so $E[XY] = E[D_1 D_2] = 3.5^2 = 12.25$. With $E[X] = 91/36$ and $E[Y] = 161/36$:

$$\operatorname{Cov}(X, Y) = 12.25 - \frac{91 \cdot 161}{36^2} = \frac{1225}{1296} \approx 0.945.$$

Positive, as expected: a large minimum forces a large maximum.

Example: continuous. For $f(x, y) = x + y$ on the unit square, $E[XY] = \iint xy(x + y) = 1/3$ and $E[X] = E[Y] = 7/12$, so $\operatorname{Cov}(X, Y) = 1/3 - 49/144 = -1/144$: a very slight negative dependence, consistent with $E[Y \mid X = x]$ decreasing in $x$ (Section 6.2).

import numpy as np

rng = np.random.default_rng(71)

# Exact covariance of (min, max) of two dice from the joint table
p = np.zeros((6, 6))
for a in range(1, 7):
    for b in range(1, 7):
        p[min(a, b) - 1, max(a, b) - 1] += 1 / 36
k = np.arange(1, 7)
EX, EY = k @ p.sum(axis=1), k @ p.sum(axis=0)
EXY = k @ p @ k                                       # sum over x, y of x * y * p(x, y)
print(f"Cov(min, max) = {EXY - EX * EY:.4f} (1225/1296 = {1225 / 1296:.4f})")

# Bilinearity check by simulation: Cov(2X + 1, -3Y) = -6 Cov(X, Y)
d = rng.integers(1, 7, (1_000_000, 2))
x, y = d.min(axis=1), d.max(axis=1)
cov = lambda a, b: np.mean((a - a.mean()) * (b - b.mean()))
print(f"Cov(X, Y) = {cov(x, y):.4f}, Cov(2X + 1, -3Y) = {cov(2 * x + 1, -3 * y):.4f} (-6 x {cov(x, y):.4f} = {-6 * cov(x, y):.4f})")

# Variance of a sum: Var(X + Y) = Var(X) + Var(Y) + 2 Cov(X, Y)
print(f"Var(X + Y) = {np.var(x + y):.4f}; Var X + Var Y + 2 Cov = {x.var() + y.var() + 2 * cov(x, y):.4f}")
print(f"Var(X + Y) = Var(D1 + D2) = 2 x 35/12 = {2 * 35 / 12:.4f}")   # min + max = sum of the dice

Cov(min, max) = 0.9452 (1225/1296 = 0.9452)
Cov(X, Y) = 0.9445, Cov(2X + 1, -3Y) = -5.6672 (-6 x 0.9445 = -5.6672)
Var(X + Y) = 5.8322; Var X + Var Y + 2 Cov = 5.8322
Var(X + Y) = Var(D1 + D2) = 2 x 35/12 = 5.8333

The last line is a useful sanity check: the minimum plus the maximum is just the sum of the two dice, whose variance is known exactly.

Key takeaways

$\operatorname{Cov}(X, Y) = E[(X - \mu_X)(Y - \mu_Y)] = E[XY] - E[X]E[Y]$; its sign shows the direction of the (linear) co-movement.
Covariance is bilinear: constants factor out, shifts disappear, sums expand term by term.
$\operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y) + 2\operatorname{Cov}(X, Y)$.
Its magnitude depends on the units; use correlation to judge strength.

🧪 Practice

$\operatorname{Var}(X) = 4$, $\operatorname{Var}(Y) = 9$ and $\operatorname{Cov}(X, Y) = -2$. Compute $\operatorname{Var}(X + Y)$, $\operatorname{Var}(X - Y)$ and $\operatorname{Cov}(3X - Y, X + 2Y)$.

Toss a fair coin 3 times. Let $X$ be the number of heads in the first two tosses and $Y$ the number of heads in the last two. Compute $\operatorname{Cov}(X, Y)$.

For $f(x, y) = 2$ on $0 < x < y < 1$, compute $\operatorname{Cov}(X, Y)$.

Interview: How would you compute the covariance of two columns over a stream of billions of rows, in one pass and without numerical trouble? Hint: the formula $E[XY] - E[X]E[Y]$ suffers from catastrophic cancellation when the means are large; look at Welford-style online updates for the co-moment.

Correlation and Its Bounds

Why it matters. Is a covariance of 0.95 strong? It depends on the units. Correlation removes the units by standardizing both variables, giving a number between $-1$ and $1$ that can be compared across problems: the correlation between features in a model, between two assets in a portfolio, between a biomarker and an outcome. It is also one of the most misinterpreted numbers in statistics.

Definition. The (Pearson) correlation coefficient is

$$\rho(X, Y) = \operatorname{Corr}(X, Y) = \frac{\operatorname{Cov}(X, Y)}{\sigma_X, \sigma_Y} = \operatorname{Cov}!\left(\frac{X - \mu_X}{\sigma_X}, \frac{Y - \mu_Y}{\sigma_Y}\right),$$

the covariance of the standardized variables. It is the population counterpart of the sample correlation $r$ of Section 2.5.

The bounds. For any $X, Y$ with positive finite variances,

$$-1 \le \rho(X, Y) \le 1.$$

Proof. Let $Z_X$ and $Z_Y$ be the standardized variables, each with variance 1. Variances are non-negative, so

$$0 \le \operatorname{Var}(Z_X \pm Z_Y) = 1 + 1 \pm 2\rho \quad\Longrightarrow\quad \mp\rho \le 1.$$

The bound is reached exactly when $\operatorname{Var}(Z_X - Z_Y) = 0$ (or $Z_X + Z_Y$), that is, when $Y = a + bX$ with probability 1: $\rho = 1$ if $b > 0$ and $\rho = -1$ if $b < 0$. The same bound also follows from the Cauchy-Schwarz inequality (Section 6.5).

Properties.

Property	Statement
Unit-free	$\rho(aX + b, cY + d) = \operatorname{sign}(ac),\rho(X, Y)$
Extreme values	$\lvert \rho \rvert = 1$ if and only if $Y$ is an exact linear function of $X$
Independence	independent $\Rightarrow \rho = 0$, but not conversely
Explained variance	for the best linear predictor of $Y$ from $X$, the fraction of $\operatorname{Var}(Y)$ explained is $\rho^2$

What correlation does not measure.

Non-linear dependence: $Y = X^2$ with $X$ symmetric about 0 has $\rho = 0$ although $Y$ is a function of $X$.
Slope: $\rho = 0.9$ can come with a slope of 0.01 or 100; correlation measures how closely points follow a line, not how steep it is.
Causation: correlation can come from a common cause, from selection, or from chance (Section 12.4).
Robustness: a single extreme point can create or destroy a large correlation. Rank correlations (Spearman, Kendall) are more robust and capture any monotone relationship.

Range restriction. Restricting the range of $X$ (for example, studying only admitted students) lowers the correlation even if the underlying relationship is unchanged, because the signal $\operatorname{Var}(E[Y \mid X])$ shrinks while the noise does not.

import numpy as np

rng = np.random.default_rng(72)
n = 200_000
corr = lambda a, b: np.corrcoef(a, b)[0, 1]

x = rng.normal(0, 1, n)
cases = {
    "Y = 2X + 1 (exact line)":      2 * x + 1,
    "Y = -0.5X (exact line)":       -0.5 * x,
    "Y = X + noise (SD 1)":         x + rng.normal(0, 1, n),
    "Y = 0.01X + tiny noise":       0.01 * x + rng.normal(0, 0.003, n),
    "Y = X^2 (exact, non-linear)":  x**2,
    "Y independent of X":           rng.normal(0, 1, n),
}
for name, y in cases.items():
    print(f"{name:28s} rho = {corr(x, y):+.3f}")

# Units do not matter; signs of the scale factors do
y = x + rng.normal(0, 1, n)
print(f"rho(X, Y) = {corr(x, y):+.3f}, rho(100X + 7, -3Y) = {corr(100 * x + 7, -3 * y):+.3f}")

# Range restriction: keep only the top 20% of X
top = x > np.quantile(x, 0.8)
print(f"full range rho = {corr(x, y):.3f}, top 20% of X only: rho = {corr(x[top], y[top]):.3f}")

Y = 2X + 1 (exact line)      rho = +1.000
Y = -0.5X (exact line)       rho = -1.000
Y = X + noise (SD 1)         rho = +0.708
Y = 0.01X + tiny noise       rho = +0.958
Y = X^2 (exact, non-linear)  rho = -0.007
Y independent of X           rho = -0.001
rho(X, Y) = +0.709, rho(100X + 7, -3Y) = -0.709
full range rho = 0.709, top 20% of X only: rho = 0.426

Notice the fourth case: a correlation of almost 0.96 with a slope of only 0.01. Strength of association and size of effect are different questions.

Key takeaways

$\rho = \operatorname{Cov}(X, Y)/(\sigma_X \sigma_Y)$ is unit-free and always lies in $[-1, 1]$.
$\lvert \rho \rvert = 1$ exactly when $Y$ is a linear function of $X$; $\rho^2$ is the variance explained by the best linear predictor.
Correlation captures linear association only; it says nothing about slope, non-linear dependence or causation.
Outliers and range restriction can change correlation dramatically.

🧪 Practice

Using $\operatorname{Cov}(\min, \max) = 1225/1296$ and $\operatorname{Var}(\min) = \operatorname{Var}(\max) = 2555/1296$, compute the correlation of the minimum and maximum of two dice.

$X$ and $Y$ have variances 1 and 4 and correlation 0.5. Find $\operatorname{Var}(X - Y)$ and $\rho(X, X - Y)$.

Show that for any three variables with $\rho(X, Y) = \rho(Y, Z) = 0.9$, $\rho(X, Z)$ must be at least 0.62. Hint: the $3 \times 3$ correlation matrix must be positive semidefinite (next topic).

Interview: Two features in your model have a correlation of 0.98. What problems could this cause, and what would you do? Hint: think about the variance of the coefficient estimates in linear regression (multicollinearity), interpretability, and options such as dropping one, combining them, or regularization (Section 10.4).

Covariance Matrices

Why it matters. With more than two variables there is one covariance for every pair, and bookkeeping becomes impossible without matrix notation. The covariance matrix collects them all, and linear algebra then answers the important questions in one line: the risk of any portfolio, the distribution of any linear combination, the directions of greatest variation (PCA), and how to simulate correlated variables. It is the central object of multivariate statistics (Section 12.2).

Definition. For a random vector $\mathbf{X} = (X_1, \dots, X_d)^\top$ with mean vector $\boldsymbol{\mu} = E[\mathbf{X}]$, the covariance matrix is the $d \times d$ matrix

$$\Sigma = \operatorname{Cov}(\mathbf{X}) = E\big[(\mathbf{X} - \boldsymbol{\mu})(\mathbf{X} - \boldsymbol{\mu})^\top\big], \qquad \Sigma_{ij} = \operatorname{Cov}(X_i, X_j).$$

The diagonal holds the variances and the off-diagonal entries the covariances:

            X1          X2          X3
  X1  [ Var(X1)    Cov(X1,X2)  Cov(X1,X3) ]
  X2  [ Cov(X2,X1) Var(X2)     Cov(X2,X3) ]      symmetric: mirror image across the diagonal
  X3  [ Cov(X3,X1) Cov(X3,X2)  Var(X3)    ]

The correlation matrix $R$ is the covariance matrix of the standardized variables: $R = D^{-1} \Sigma D^{-1}$, with $D$ the diagonal matrix of standard deviations. It has ones on the diagonal.

Linear transformations. If $\mathbf{Y} = A\mathbf{X} + \mathbf{b}$ for a fixed matrix $A$ and vector $\mathbf{b}$, then

$$E[\mathbf{Y}] = A\boldsymbol{\mu} + \mathbf{b}, \qquad \operatorname{Cov}(\mathbf{Y}) = A, \Sigma, A^\top.$$

In particular, for a single linear combination $\mathbf{a}^\top \mathbf{X} = \sum_i a_i X_i$,

$$\operatorname{Var}(\mathbf{a}^\top \mathbf{X}) = \mathbf{a}^\top \Sigma, \mathbf{a} = \sum_i \sum_j a_i a_j \Sigma_{ij}.$$

This is the portfolio variance formula: with weights $\mathbf{w}$ and return covariance $\Sigma$, the portfolio risk is $\sqrt{\mathbf{w}^\top \Sigma \mathbf{w}}$.

Positive semidefinite. Because the variance of any linear combination is non-negative, $\mathbf{a}^\top \Sigma \mathbf{a} \ge 0$ for every $\mathbf{a}$: every covariance matrix is symmetric positive semidefinite (PSD), and all its eigenvalues are $\ge 0$. A zero eigenvalue means some linear combination of the variables is constant, so one variable is an exact linear function of the others. Conversely, a symmetric matrix with a negative eigenvalue cannot be a covariance matrix, which is why correlations cannot be chosen freely (practice 3 of the previous topic).

Three decompositions and what they are for.

Decomposition	Form	Use
Eigendecomposition	$\Sigma = Q \Lambda Q^\top$	principal components: directions $Q$ and their variances $\Lambda$
Cholesky	$\Sigma = L L^\top$, $L$ lower triangular	simulation: if $\mathbf{Z}$ has covariance $I$, then $L\mathbf{Z}$ has covariance $\Sigma$
Inverse $\Sigma^{-1}$ (precision matrix)		Mahalanobis distance, normal density, partial correlations

Mahalanobis distance. The distance $\sqrt{(\mathbf{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})}$ measures how unusual a point is, taking scale and correlation into account. A person 190 cm tall and 50 kg is not extreme on either variable alone, but very unusual for the combination; the Mahalanobis distance catches this, the Euclidean distance does not. It is a standard tool for multivariate outlier detection.

import numpy as np

rng = np.random.default_rng(73)

# Annual returns of three assets: SDs and correlations
sd = np.array([0.20, 0.15, 0.05])                     # stocks, real estate, bonds
R = np.array([[1.0, 0.6, -0.2],
              [0.6, 1.0,  0.1],
              [-0.2, 0.1, 1.0]])
Sigma = np.diag(sd) @ R @ np.diag(sd)                 # covariance = D R D
print("eigenvalues of Sigma:", np.round(np.linalg.eigvalsh(Sigma), 5), "(all >= 0: valid)")

w = np.array([0.5, 0.2, 0.3])                          # portfolio weights
print(f"portfolio SD = sqrt(w' Sigma w) = {np.sqrt(w @ Sigma @ w):.4f}; "
      f"if all correlations were 1: {w @ sd:.4f}")

# Simulate correlated returns with the Cholesky factor
L = np.linalg.cholesky(Sigma)
Z = rng.standard_normal((1_000_000, 3))               # independent, covariance I
X = Z @ L.T                                           # each row is L z, covariance L L' = Sigma
print("sample correlation matrix:\n", np.round(np.corrcoef(X, rowvar=False), 3))
print(f"simulated portfolio SD = {(X @ w).std():.4f}")

# An "impossible" correlation matrix: not PSD
bad = np.array([[1, 0.9, -0.9], [0.9, 1, 0.9], [-0.9, 0.9, 1]])
print("eigenvalues of a bad 'correlation matrix':", np.round(np.linalg.eigvalsh(bad), 3))

# Mahalanobis vs Euclidean distance for (height cm, weight kg)
mu, S = np.array([175, 75]), np.array([[49, 42], [42, 100]])   # SD 7 and 10, correlation 0.6
Sinv = np.linalg.inv(S)
for pt in [(190, 50), (190, 95)]:
    d = np.array(pt) - mu
    print(f"{pt}: Euclidean {np.linalg.norm(d):5.1f}, Mahalanobis {np.sqrt(d @ Sinv @ d):.2f}")

eigenvalues of Sigma: [0.00215 0.01155 0.0513 ] (all >= 0: valid)
portfolio SD = sqrt(w' Sigma w) = 0.1192; if all correlations were 1: 0.1450
sample correlation matrix:
 [[ 1.     0.6   -0.199]
 [ 0.6    1.     0.101]
 [-0.199  0.101  1.   ]]
simulated portfolio SD = 0.1192
eigenvalues of a bad 'correlation matrix': [-0.8  1.9  1.9]
(190, 50): Euclidean  29.2, Mahalanobis 5.19
(190, 95): Euclidean  25.0, Mahalanobis 2.32

The tall, light person is closer to the mean in Euclidean terms than the tall, heavy one, but far more unusual once the correlation between height and weight is taken into account.

Key takeaways

The covariance matrix $\Sigma$ holds variances on the diagonal and covariances off it; it is symmetric and positive semidefinite.
$\operatorname{Cov}(A\mathbf{X} + \mathbf{b}) = A \Sigma A^\top$ and $\operatorname{Var}(\mathbf{a}^\top \mathbf{X}) = \mathbf{a}^\top \Sigma \mathbf{a}$: the portfolio variance formula.
Eigenvectors give principal components; the Cholesky factor generates correlated samples; the inverse defines the Mahalanobis distance.
Not every symmetric matrix with ones on the diagonal is a valid correlation matrix: check the eigenvalues.

🧪 Practice

$X_1, X_2$ have variances 4 and 1 and covariance 1. Write $\Sigma$ and compute $\operatorname{Var}(X_1 - 2X_2)$ with the formula $\mathbf{a}^\top \Sigma \mathbf{a}$.

For $\mathbf{Y} = (X_1 + X_2, X_1 - X_2)$ with $\Sigma$ from exercise 1, find $\operatorname{Cov}(\mathbf{Y})$ as $A \Sigma A^\top$.

Find all values of $c$ for which the matrix with ones on the diagonal and $c$ in every off-diagonal entry is a valid $3 \times 3$ correlation matrix.

Interview: You estimate a $500 \times 500$ covariance matrix from 250 days of stock returns and the optimizer produces extreme portfolio weights. What is going on, and how would you fix it? Hint: with fewer observations than variables, the sample covariance matrix is singular and its smallest eigenvalues are badly underestimated; look at shrinkage (Ledoit-Wolf) or factor models.

Uncorrelated vs. Independent

Why it matters. "The variables are uncorrelated, so they are independent" is one of the most common mistakes in data analysis. Zero correlation only rules out a linear trend; the variables can still be strongly, even perfectly, dependent. Confusing the two leads to dropped features that actually carry signal, to wrong variance calculations, and to wrong conclusions about risk.

The one-way implication. If $X$ and $Y$ are independent and have finite variances, then $E[XY] = E[X]E[Y]$, so $\operatorname{Cov}(X, Y) = 0$:

$$\text{independent} ;\Longrightarrow; \text{uncorrelated}, \qquad \text{uncorrelated} ;\not\Longrightarrow; \text{independent}.$$

Counterexample 1: a parabola. Let $X \sim$ Uniform$(-1, 1)$ and $Y = X^2$. Then $\operatorname{Cov}(X, Y) = E[X^3] - E[X]E[X^2] = 0$ because both $E[X^3]$ and $E[X]$ vanish by symmetry. Yet $Y$ is completely determined by $X$. The positive and negative trends on the two halves cancel in the average.

Counterexample 2: random sign. Let $X \sim N(0, 1)$ and let $S = \pm 1$ with probability 1/2 each, independent of $X$. Then $Y = SX$ is also $N(0, 1)$, and $\operatorname{Cov}(X, Y) = E[S]E[X^2] = 0$. But $\lvert Y \rvert = \lvert X \rvert$ always, so they are highly dependent. This example also shows that two normal variables need not be jointly normal: $X + Y$ equals 0 with probability 1/2, which no normal variable does.

When uncorrelated does imply independent.

Situation	Uncorrelated $\Rightarrow$ independent?
$(X, Y)$ jointly (bivariate) normal	yes (Section 6.4)
$X$ and $Y$ both Bernoulli (binary)	yes: the 2 x 2 table is fixed by the margins and one covariance
each variable separately normal, but not jointly normal	no (random sign example)
general distributions	no

Stronger notions that sit in between. Mean independence, $E[Y \mid X] = E[Y]$, is stronger than zero correlation and weaker than independence: it rules out any trend in the mean of $Y$, not only linear ones, but allows the spread to depend on $X$ (as in heteroscedastic errors, Section 6.2). The full hierarchy is

$$\text{independent} ;\Rightarrow; E[Y \mid X] = E[Y] ;\Rightarrow; \operatorname{Cov}(X, Y) = 0.$$

Detecting non-linear dependence. Compare the correlation of transformed variables (such as $X^2$ and $Y$, or the ranks), compare conditional distributions across bins of $X$, or use measures designed to detect any dependence: mutual information, distance correlation, or the maximal information coefficient.

import numpy as np

rng = np.random.default_rng(74)
n = 1_000_000
corr = lambda a, b: np.corrcoef(a, b)[0, 1]

# Counterexample 1: Y = X^2 with X symmetric
x = rng.uniform(-1, 1, n)
y = x**2
print(f"parabola:    corr(X, Y) = {corr(x, y):+.4f}, corr(|X|, Y) = {corr(np.abs(x), y):+.4f}")
print(f"             P(Y > 0.5) = {np.mean(y > 0.5):.3f}, P(Y > 0.5 | |X| < 0.5) = {np.mean(y[np.abs(x) < 0.5] > 0.5):.3f}")

# Counterexample 2: random sign, both marginals N(0, 1)
x = rng.standard_normal(n)
y = rng.choice([-1, 1], n) * x
print(f"random sign: corr(X, Y) = {corr(x, y):+.4f}, corr(X^2, Y^2) = {corr(x**2, y**2):+.4f}")
print(f"             SD of Y = {y.std():.3f}, P(X + Y = 0) = {np.mean(x + y == 0):.3f} (a normal sum would give 0)")

# Mean independent but not independent: E[Y | X] = 0, spread grows with |X|
y = rng.standard_normal(n) * (0.2 + np.abs(x))
inner, outer = np.abs(x) < 0.5, np.abs(x) > 1.5
print(f"heteroscedastic: corr = {corr(x, y):+.4f}; mean of Y: {y[inner].mean():+.3f} vs {y[outer].mean():+.3f}; "
      f"SD of Y: {y[inner].std():.3f} vs {y[outer].std():.3f}")

parabola:    corr(X, Y) = -0.0006, corr(|X|, Y) = +0.9682
             P(Y > 0.5) = 0.292, P(Y > 0.5 | |X| < 0.5) = 0.000
random sign: corr(X, Y) = +0.0014, corr(X^2, Y^2) = +1.0000
             SD of Y = 1.000, P(X + Y = 0) = 0.500 (a normal sum would give 0)
heteroscedastic: corr = -0.0009; mean of Y: +0.001 vs -0.003; SD of Y: 0.469 vs 2.175

In all three cases the correlation is essentially zero, yet knowing $X$ changes the distribution of $Y$: the conditional probabilities, the squared values or the spread all differ.

Key takeaways

Independence implies zero correlation; zero correlation does not imply independence.
Correlation only detects linear trends; symmetric non-linear relationships (like $Y = X^2$) can have $\rho = 0$.
For jointly normal variables (and for pairs of binary variables) uncorrelated does mean independent; marginally normal is not enough.
Hierarchy: independent $\Rightarrow$ mean independent $\Rightarrow$ uncorrelated.

🧪 Practice

Let $X$ take values $-1, 0, 1$ with probability $1/3$ each, and $Y = \mathbf{1}{X = 0}$. Show that $X$ and $Y$ are uncorrelated but not independent.

Prove that two Bernoulli variables with $\operatorname{Cov}(X, Y) = 0$ are independent. Hint: show that $P(X = 1, Y = 1) = P(X = 1)P(Y = 1)$, then deduce the other three cells.

Let $\Theta \sim$ Uniform$(0, 2\pi)$, $X = \cos\Theta$, $Y = \sin\Theta$. Are $X$ and $Y$ correlated? Are they independent?

Interview: A feature has near-zero correlation with the target, so a teammate proposes dropping it. What would you check first? Hint: plot the target against binned values of the feature, look at rank or non-linear dependence measures, and think about interactions with other features.

Copulas

Why it matters. A joint distribution mixes two kinds of information: how each variable behaves on its own (the marginals) and how they move together (the dependence). Copulas separate the two. This lets you model, for example, insurance losses with heavy-tailed marginals and a dependence structure chosen separately, or simulate realistic correlated inputs for a risk model. Copulas also carry a cautionary tale: the Gaussian copula, used widely to price mortgage-backed securities before 2008, assumes that extreme losses almost never happen together, and they did.

The idea: transform to uniforms. By the probability integral transform (Section 4.4), if $X$ has a continuous CDF $F_X$, then $U = F_X(X)$ is Uniform(0, 1). Applying each marginal CDF to its own variable turns $(X, Y)$ into a pair $(U, V)$ with uniform marginals but the same dependence: the transformation is increasing, so it does not change which values go together. The joint distribution of $(U, V)$ is the copula.

  (X, Y) with any marginals             (U, V) = (F_X(X), F_Y(Y))          (X', Y') with new marginals
  e.g. income and spend        --->     uniform marginals on [0,1]^2  --->  e.g. exponential and lognormal
                                        dependence only: the COPULA          same dependence, new marginals

Definition and Sklar's theorem. A (bivariate) copula is a joint CDF $C(u, v)$ on $[0, 1]^2$ with Uniform(0, 1) marginals. Sklar's theorem says that every joint CDF can be written as

$$F(x, y) = C\big(F_X(x), F_Y(y)\big)$$

for some copula $C$, which is unique when the marginals are continuous. Conversely, plugging any marginals into any copula gives a valid joint distribution. Marginals and dependence can be chosen separately.

Important copulas.

Copula	Formula or construction	Dependence
Independence	$C(u, v) = uv$	none
Comonotone (upper bound)	$C(u, v) = \min(u, v)$	perfect positive: $V = U$
Countermonotone (lower bound)	$C(u, v) = \max(u + v - 1, 0)$	perfect negative: $V = 1 - U$
Gaussian	$(\Phi(Z_1), \Phi(Z_2))$ with $(Z_1, Z_2)$ bivariate normal, correlation $\rho$	symmetric, no tail dependence
Student $t$	same with a bivariate $t$ distribution	symmetric, extremes occur together
Clayton	$(u^{-\theta} + v^{-\theta} - 1)^{-1/\theta}$	strong dependence in the lower tail
Gumbel	Archimedean family	strong dependence in the upper tail

Every copula lies between the countermonotone and comonotone bounds (the Fréchet-Hoeffding bounds).

Tail dependence. The upper tail dependence coefficient $\lambda_U = \lim_{q \to 1} P(V > q \mid U > q)$ measures whether extreme values occur together. For the Gaussian copula with $\rho < 1$, $\lambda_U = 0$: however strong the correlation, joint extremes become relatively rare far enough in the tail. For the $t$ copula, $\lambda_U > 0$. In finance and insurance, where joint extreme losses are what matter, this difference is decisive.

Rank correlations live in the copula. Pearson correlation depends on the marginals, but Spearman's rank correlation and Kendall's tau depend only on the copula, so they do not change under increasing transformations of either variable.

Sampling from a copula model (Gaussian copula).

Draw $(Z_1, Z_2)$ bivariate normal with correlation $\rho$.
Set $U = \Phi(Z_1)$, $V = \Phi(Z_2)$ (uniform marginals, Gaussian dependence).
Set $X = F_X^{-1}(U)$, $Y = F_Y^{-1}(V)$ with any desired marginals (inverse transform sampling, Section 4.4).

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(75)
n, rho = 400_000, 0.7
Phi = np.vectorize(NormalDist().cdf)                  # standard normal CDF

def ranks(a):
    """Pseudo-observations: ranks scaled to (0, 1), the empirical version of F(X)."""
    return (np.argsort(np.argsort(a)) + 1) / (len(a) + 1)

def spearman(a, b):
    return np.corrcoef(ranks(a), ranks(b))[0, 1]

# Gaussian copula with exponential and lognormal marginals
z = rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]], n)
u, v = Phi(z[:, 0]), Phi(z[:, 1])                     # uniform marginals, Gaussian dependence
x = -np.log(1 - u) * 10                               # Exponential with mean 10 (inverse CDF)
y = np.exp(1 + 0.8 * z[:, 1])                         # Lognormal: exp of a normal, same as F^-1(V)
print(f"Pearson:  normal pair {np.corrcoef(z.T)[0, 1]:.3f}, transformed pair {np.corrcoef(x, y)[0, 1]:.3f}")
print(f"Spearman: normal pair {spearman(z[:, 0], z[:, 1]):.3f}, transformed pair {spearman(x, y):.3f}")

# Tail dependence: Gaussian copula versus t copula (3 degrees of freedom), same correlation matrix
w = np.sqrt(3 / rng.chisquare(3, n))                  # common random scale makes extremes cluster
t = z * w[:, None]                                    # bivariate t sample
gu, gv = ranks(z[:, 0]), ranks(z[:, 1])
tu, tv = ranks(t[:, 0]), ranks(t[:, 1])
print("q       P(V > q | U > q): Gaussian   t(3)")
for q in [0.9, 0.99, 0.999]:
    g = np.mean(gv[gu > q] > q)
    s = np.mean(tv[tu > q] > q)
    print(f"{q:<7} {'':18s}{g:.3f}     {s:.3f}")

Pearson:  normal pair 0.701, transformed pair 0.641
Spearman: normal pair 0.684, transformed pair 0.684
q       P(V > q | U > q): Gaussian   t(3)
0.9                       0.469     0.533
0.99                      0.271     0.461
0.999                     0.180     0.482

Changing the marginals changed the Pearson correlation but left the Spearman correlation intact, because the copula did not change. Further into the tail, the conditional probability of a joint extreme keeps falling under the Gaussian copula, while under the $t$ copula it stays high: same correlation, very different risk of simultaneous disasters.

Key takeaways

A copula is a joint distribution with uniform marginals; it captures dependence separately from the marginals.
Sklar's theorem: $F(x, y) = C(F_X(x), F_Y(y))$, and any marginals can be combined with any copula.
Rank correlations (Spearman, Kendall) depend only on the copula; Pearson correlation also depends on the marginals.
The Gaussian copula has no tail dependence; the $t$, Clayton and Gumbel copulas can model extremes that happen together.

🧪 Practice

Verify that $C(u, v) = uv$ has uniform marginals, and that plugging in marginal CDFs gives $F(x, y) = F_X(x)F_Y(y)$, the independence case.

If $Y = e^X$, what is the copula of $(X, Y)$? What about $Y = -X^3$?

Simulate a Clayton copula with $\theta = 2$ (draw $W \sim$ Gamma(1/2, 1) and independent $E_1, E_2 \sim$ Exponential(1), then $U_i = (1 + E_i/W)^{-1/\theta}$) and estimate $P(V < q \mid U < q)$ for $q = 0.01$.

Interview: Why was the Gaussian copula criticized after the 2008 financial crisis? Hint: think about what a zero tail dependence coefficient implies for the probability that many mortgages default at the same time, and about correlations estimated from calm periods.

6.4 Sums and Multivariate Models

Sums of random variables are everywhere in statistics, and a few multivariate models (the multinomial, the multivariate normal, the joint behavior of sorted values) cover a large share of practical problems.

Expectation and Variance of Sums

Why it matters. A sample mean is a sum divided by $n$; a portfolio return is a weighted sum; the total time of a multi-step job is a sum of step times; the number of successes is a sum of indicators. Being able to compute the mean and variance of a sum, without finding its distribution, is the single most useful skill in applied probability. It also explains the most important fact of statistics: averaging reduces noise, by a factor of $\sqrt n$ when the terms are independent, and by much less when they are not.

Mean of a sum. Expectation is linear, always, with no independence required:

$$E\Big[\sum_{i=1}^n a_i X_i\Big] = \sum_{i=1}^n a_i, E[X_i].$$

Variance of a sum. Expanding the square of a sum (bilinearity of covariance, Section 6.3):

$$\operatorname{Var}\Big(\sum_{i=1}^n a_i X_i\Big) = \sum_{i=1}^n a_i^2 \operatorname{Var}(X_i) + 2 \sum_{i < j} a_i a_j \operatorname{Cov}(X_i, X_j) = \mathbf{a}^\top \Sigma, \mathbf{a}.$$

There are $n$ variance terms and $n(n - 1)/2$ covariance terms, so for large $n$ the covariances usually dominate unless they are zero.

Situation	$\operatorname{Var}(X_1 + \dots + X_n)$	$\operatorname{Var}(\bar X)$
Independent (or just uncorrelated), common variance $\sigma^2$	$n\sigma^2$	$\sigma^2/n$
Common variance $\sigma^2$, all pairwise correlations $\rho$	$n\sigma^2 + n(n - 1)\rho\sigma^2$	$\dfrac{\sigma^2}{n}\big(1 + (n - 1)\rho\big)$
Perfectly correlated ($\rho = 1$)	$n^2\sigma^2$	$\sigma^2$: averaging does not help at all

Two lessons.

The $\sqrt n$ law. For independent terms, the SD of the sum grows like $\sqrt n$ and the SD of the mean shrinks like $\sigma/\sqrt n$. To halve the standard error, you need four times as much data.
Correlation puts a floor under the noise. With equicorrelated terms, $\operatorname{Var}(\bar X) \to \rho\sigma^2$ as $n \to \infty$. Averaging 1,000 opinions from people who all read the same newspaper is far less informative than averaging 1,000 independent opinions, and diversifying over 500 stocks that all move with the market does not remove market risk. The design effect of cluster samples (Section 7.1) is this formula in action.

The indicator method. To find the mean and variance of a count, write it as a sum of indicators $X = \sum_i \mathbf{1}_{A_i}$. Then $E[X] = \sum_i P(A_i)$, and

$$\operatorname{Var}(X) = \sum_i P(A_i)\big(1 - P(A_i)\big) + 2\sum_{i < j} \big(P(A_i \cap A_j) - P(A_i)P(A_j)\big).$$

Example: the matching problem. $n$ people throw their hats in a pile and each takes one at random. Let $X$ be the number who get their own hat, with $A_i$ = "person $i$ gets their own hat". Then $P(A_i) = 1/n$ and $P(A_i \cap A_j) = \frac{1}{n(n - 1)}$, so

$$E[X] = n \cdot \frac1n = 1, \qquad \operatorname{Var}(X) = n \cdot \frac1n\Big(1 - \frac1n\Big) + n(n - 1)\Big(\frac{1}{n(n - 1)} - \frac{1}{n^2}\Big) = 1,$$

for every $n \ge 2$. Computing the distribution of $X$ directly requires inclusion-exclusion; the mean and variance need none of it.

import numpy as np

rng = np.random.default_rng(76)

# Averaging equicorrelated measurements: Var(mean) = sigma^2/n (1 + (n - 1) rho)
sigma, reps = 1.0, 20_000
for rho in [0.0, 0.1, 0.5]:
    out = []
    for n in [1, 10, 100, 1000]:
        common = rng.standard_normal((reps, 1))       # shared component creates correlation rho
        own = rng.standard_normal((reps, n))
        x = sigma * (np.sqrt(rho) * common + np.sqrt(1 - rho) * own)
        sim = x.mean(axis=1).std()
        theory = sigma * np.sqrt((1 + (n - 1) * rho) / n)
        out.append(f"n={n:<4d} {sim:.3f} ({theory:.3f})")
    print(f"rho = {rho}: SD of mean  " + "   ".join(out))

# Matching problem: number of fixed points of a random permutation
for n in [3, 10, 1000]:
    perms = rng.permuted(np.tile(np.arange(n), (20_000, 1)), axis=1)   # each row a random permutation
    fixed = (perms == np.arange(n)).sum(axis=1)
    print(f"n = {n:4d}: E[X] = {fixed.mean():.3f}, Var(X) = {fixed.var():.3f}")

rho = 0.0: SD of mean  n=1    0.997 (1.000)   n=10   0.316 (0.316)   n=100  0.100 (0.100)   n=1000 0.032 (0.032)
rho = 0.1: SD of mean  n=1    1.000 (1.000)   n=10   0.433 (0.436)   n=100  0.327 (0.330)   n=1000 0.317 (0.318)
rho = 0.5: SD of mean  n=1    0.999 (1.000)   n=10   0.742 (0.742)   n=100  0.711 (0.711)   n=1000 0.709 (0.707)
n =    3: E[X] = 1.001, Var(X) = 0.994
n =   10: E[X] = 1.007, Var(X) = 1.008
n = 1000: E[X] = 0.986, Var(X) = 0.997

With $\rho = 0.1$, averaging 1,000 measurements is barely better than averaging 100: the SD of the mean is stuck near $\sqrt{0.1} \approx 0.32$.

Key takeaways

$E[\sum a_i X_i] = \sum a_i E[X_i]$ always; $\operatorname{Var}(\sum a_i X_i) = \sum a_i^2 \operatorname{Var}(X_i) + 2\sum_{i<j} a_i a_j \operatorname{Cov}(X_i, X_j)$.
For independent terms, $\operatorname{SD}(\bar X) = \sigma/\sqrt n$; correlated terms put a floor of $\rho\sigma^2$ under $\operatorname{Var}(\bar X)$.
Write counts as sums of indicators to get their mean and variance without their distribution.

🧪 Practice

A job has three stages with mean times 2, 5, 3 minutes and SDs 1, 2, 1 minutes. Find the mean and SD of the total time if the stages are independent, and if stages 1 and 2 have correlation 0.5.

Roll a die 60 times. Use indicators to find the mean and variance of the number of distinct faces that never appear.

In a group of 30 people, let $X$ be the number of pairs that share a birthday. Find $E[X]$ and $\operatorname{Var}(X)$ (365 equally likely days). Are the pair indicators independent? Are they uncorrelated?

Interview: You average the predictions of 10 models, each with error variance $\sigma^2$, and the error reduction is far smaller than the factor of 10 you hoped for. Why? Hint: the errors of models trained on the same data are positively correlated; apply the equicorrelated formula, and think about how bagging and random feature subsets try to lower $\rho$.

Convolution of Distributions

Why it matters. The mean and variance of a sum are easy, but sometimes we need its whole distribution: the probability that total demand exceeds capacity, the distribution of the sum of two dice, the total delay of a request crossing several services. For independent variables, the distribution of the sum is the convolution of the individual distributions. Convolution also explains why some families are closed under addition and why sums tend toward the normal distribution.

Discrete convolution. If $X$ and $Y$ are independent and integer-valued, the event $X + Y = s$ splits by the value of $X$:

$$p_{X+Y}(s) = \sum_k P(X = k), P(Y = s - k) = \sum_k p_X(k), p_Y(s - k).$$

Continuous convolution. For independent continuous $X$ and $Y$ (this is the transformation $U = X + Y$, $V = X$ of Section 6.1 with $V$ integrated out):

$$f_{X+Y}(s) = \int_{-\infty}^{\infty} f_X(x), f_Y(s - x), dx.$$

Picture: flip and slide. Flip $f_Y$, shift it by $s$, multiply by $f_X$ and add up the overlap. As $s$ moves, the overlap traces the density of the sum.

  Sum of two dice: two flat distributions convolve into a triangle

   s   36 p(s)
   2     1     #
   3     2     ##
   4     3     ###
   5     4     ####
   6     5     #####
   7     6     ######
   8     5     #####
   9     4     ####
  10     3     ###
  11     2     ##
  12     1     #

  Sum of two Uniform(0, 1): flat * flat = triangle on [0, 2], peak at 1
  Sum of three:             piecewise quadratic, already bell-shaped

Families closed under addition. For independent variables:

$X$	$Y$	$X + Y$
Binomial($n$, $p$)	Binomial($m$, $p$)	Binomial($n + m$, $p$)
Poisson($\lambda$)	Poisson($\mu$)	Poisson($\lambda + \mu$)
Normal($\mu_1$, $\sigma_1^2$)	Normal($\mu_2$, $\sigma_2^2$)	Normal($\mu_1 + \mu_2$, $\sigma_1^2 + \sigma_2^2$)
Gamma($\alpha$, $\beta$)	Gamma($\gamma$, $\beta$)	Gamma($\alpha + \gamma$, $\beta$)
Exponential($\lambda$), $n$ copies		Gamma($n$, $\lambda$) (Erlang)
$\chi^2_m$	$\chi^2_k$	$\chi^2_{m+k}$
Geometric($p$), $r$ copies		Negative binomial($r$, $p$)

Note the conditions: binomials need the same $p$, gammas the same rate. A uniform plus a uniform is not uniform.

Convolution becomes multiplication. For independent $X$ and $Y$, the MGF (and characteristic function) of the sum is the product: $M_{X+Y}(t) = M_X(t),M_Y(t)$ (Section 4.3). Every row of the table above can be proved in one line this way, for example $e^{\lambda(e^t - 1)} e^{\mu(e^t - 1)} = e^{(\lambda + \mu)(e^t - 1)}$. Probability generating functions do the same for integer variables, and the fast Fourier transform makes numerical convolution fast.

Repeated convolution and the CLT. Convolving a distribution with itself many times smooths it and, after centering and scaling, it approaches the normal shape. This is the central limit theorem (Section 6.5) seen through convolution.

import numpy as np

# Sum of dice by repeated discrete convolution
die = np.full(6, 1 / 6)                               # P(1), ..., P(6)
dist = die.copy()
for n_dice in range(2, 5):
    dist = np.convolve(dist, die)                     # distribution of the sum of n_dice dice
    if n_dice == 2:
        print("two dice, 36 x p(s) for s = 2..12:", np.round(dist * 36).astype(int))
support = np.arange(4, 25)
print(f"four dice: P(sum >= 20) = {dist[support >= 20].sum():.4f}, mean = {support @ dist:.2f}")

# Continuous convolution on a grid: sum of uniforms
h = 0.001
u = np.ones(int(1 / h))                              # Uniform(0, 1) density sampled on a grid
f2 = np.convolve(u, u) * h                           # density of U1 + U2 (Riemann sum of the integral)
f3 = np.convolve(f2, u) * h                          # density of U1 + U2 + U3
for s in [0.5, 1.0, 1.5]:
    print(f"f_(U1+U2)({s}) = {f2[int(s / h)]:.3f}   (triangle: {min(s, 2 - s):.3f})")
print(f"f_(U1+U2+U3)(1.5) = {f3[int(1.5 / h)]:.3f}   (exact 0.75)")

# Closure: Poisson(2) + Poisson(3) = Poisson(5)
from math import exp, factorial
pois = lambda lam, k: np.array([exp(-lam) * lam**i / factorial(i) for i in range(k)])
conv = np.convolve(pois(2, 30), pois(3, 30))[:10]
print("max |Poisson(2) * Poisson(3) - Poisson(5)| on 0..9:", f"{np.abs(conv - pois(5, 10)).max():.1e}")

two dice, 36 x p(s) for s = 2..12: [1 2 3 4 5 6 5 4 3 2 1]
four dice: P(sum >= 20) = 0.0540, mean = 14.00
f_(U1+U2)(0.5) = 0.501   (triangle: 0.500)
f_(U1+U2)(1.0) = 0.999   (triangle: 1.000)
f_(U1+U2)(1.5) = 0.499   (triangle: 0.500)
f_(U1+U2+U3)(1.5) = 0.750   (exact 0.75)
max |Poisson(2) * Poisson(3) - Poisson(5)| on 0..9: 2.8e-17

Key takeaways

For independent variables, the distribution of the sum is the convolution: $p_{X+Y}(s) = \sum_k p_X(k)p_Y(s - k)$, or the integral analogue for densities.
Binomial (same $p$), Poisson, normal, gamma (same rate) and chi-square families are closed under addition of independent members.
MGFs and characteristic functions turn convolution into multiplication, which is the fastest way to identify the distribution of a sum.
Repeated convolution produces the bell shape of the central limit theorem.

🧪 Practice

Compute the distribution of the sum of a fair die and a fair coin scored 0/1. Is it uniform?

Derive the density of the sum of two independent Exponential($\lambda$) variables with the convolution integral, and check that it is Gamma(2, $\lambda$).

Use MGFs to show that the sum of independent $N(\mu_1, \sigma_1^2)$ and $N(\mu_2, \sigma_2^2)$ variables is normal, and find its parameters.

Interview: A request passes through three microservices whose latency distributions you know from histograms. How would you compute the distribution of the end-to-end latency, and its 99th percentile? Hint: if the latencies are independent, convolve the discretized histograms (or simulate by resampling); then ask whether independence is realistic when all services share the same load spikes.

Multinomial Distribution

Why it matters. The binomial counts successes when each trial has two outcomes. Many trials have more: a die has six faces, a survey answer has five categories, a word in a document is one of thousands of vocabulary entries, a customer chooses one of several products. The multinomial distribution counts how many trials fall into each category, and it is the model behind chi-square goodness-of-fit tests (Section 9.3), bag-of-words text models, and categorical outcomes in general.

Setup. Perform $n$ independent trials, each landing in one of $k$ categories with probabilities $p_1, \dots, p_k$ (summing to 1). Let $N_i$ be the number of trials in category $i$, so $N_1 + \dots + N_k = n$. Then $(N_1, \dots, N_k) \sim$ Multinomial($n$; $p_1, \dots, p_k$), with

$$P(N_1 = n_1, \dots, N_k = n_k) = \frac{n!}{n_1!, n_2! \cdots n_k!}, p_1^{n_1} p_2^{n_2} \cdots p_k^{n_k}, \qquad n_1 + \dots + n_k = n.$$

The multinomial coefficient counts the arrangements of the trials with the given counts (Section 1.3), and each arrangement has probability $p_1^{n_1} \cdots p_k^{n_k}$. With $k = 2$ this is the binomial distribution.

Example. Roll a die 12 times. The probability that every face appears exactly twice is $\frac{12!}{(2!)^6} (1/6)^{12} \approx 0.0034$: even the "most balanced" outcome is rare.

Moments and dependence.

Quantity	Value	Why
Marginal of $N_i$	Binomial($n$, $p_i$)	each trial is "category $i$" or "not"
$E[N_i]$	$np_i$
$\operatorname{Var}(N_i)$	$np_i(1 - p_i)$
$\operatorname{Cov}(N_i, N_j)$, $i \ne j$	$-np_ip_j$	the counts compete for a fixed total $n$
$\operatorname{Corr}(N_i, N_j)$	$-\sqrt{\dfrac{p_ip_j}{(1 - p_i)(1 - p_j)}}$	does not depend on $n$

The negative covariance follows from the indicator method: $N_i = \sum_t \mathbf{1}{\text{trial } t \text{ in } i}$, and within one trial the indicators of categories $i$ and $j$ cannot both be 1, giving a covariance of $0 - p_ip_j$ per trial.

Useful structural facts.

Merging categories: combining categories adds their probabilities and counts; the result is still multinomial.
Conditioning: given $N_1 = n_1$, the remaining counts are Multinomial($n - n_1$; $p_2, \dots, p_k$ renormalized to sum to 1).
Poissonization: if the total number of trials is $N \sim$ Poisson($\lambda$) instead of fixed, the counts become independent Poisson($\lambda p_i$). Conversely, independent Poisson counts conditioned on their total are multinomial. This is often the easiest way to analyze multinomial problems.

The chi-square connection. For large $n$, $\sum_i (N_i - np_i)^2 / (np_i)$ is approximately $\chi^2_{k-1}$: the Pearson goodness-of-fit statistic. One degree of freedom is lost because the counts must add to $n$.

import numpy as np
from math import factorial, prod

rng = np.random.default_rng(77)

# Probability that 12 die rolls show each face exactly twice
p_exact = factorial(12) / 2**6 * (1 / 6) ** 12
counts = rng.multinomial(12, [1 / 6] * 6, size=1_000_000)
print(f"P(each face twice): exact {p_exact:.5f}, simulated {np.mean((counts == 2).all(axis=1)):.5f}")

# Moments: n = 100 customers choosing among 3 products
n, p = 100, np.array([0.5, 0.3, 0.2])
N = rng.multinomial(n, p, size=500_000)
print("means:", np.round(N.mean(axis=0), 2), " theory:", n * p)
C = np.cov(N, rowvar=False)
print(f"Cov(N1, N2) = {C[0, 1]:.2f} (theory {-n * p[0] * p[1]:.2f}), "
      f"Corr(N1, N2) = {np.corrcoef(N[:, 0], N[:, 1])[0, 1]:.3f} "
      f"(theory {-np.sqrt(p[0] * p[1] / ((1 - p[0]) * (1 - p[1]))):.3f})")

# Poissonization: with a Poisson total, the category counts become independent
total = rng.poisson(n, 500_000)
Np = np.array([rng.multinomial(t, p) for t in total[:100_000]])
print(f"Poisson total: Corr(N1, N2) = {np.corrcoef(Np[:, 0], Np[:, 1])[0, 1]:+.3f}, "
      f"Var(N1) = {Np[:, 0].var():.1f} (Poisson mean {n * p[0]:.0f})")

# Pearson statistic is approximately chi-square with k - 1 = 2 degrees of freedom
X2 = (((N - n * p) ** 2) / (n * p)).sum(axis=1)
print(f"Pearson X^2: mean {X2.mean():.3f} (chi2_2 mean 2), P(X^2 > 5.991) = {np.mean(X2 > 5.991):.4f} (0.05)")

P(each face twice): exact 0.00344, simulated 0.00338
means: [49.99 30.01 20.  ]  theory: [50. 30. 20.]
Cov(N1, N2) = -15.02 (theory -15.00), Corr(N1, N2) = -0.655 (theory -0.655)
Poisson total: Corr(N1, N2) = +0.007, Var(N1) = 49.9 (Poisson mean 50)
Pearson X^2: mean 1.999 (chi2_2 mean 2), P(X^2 > 5.991) = 0.0490 (0.05)

Key takeaways

Multinomial($n$; $p_1, \dots, p_k$) counts $n$ independent trials over $k$ categories; with $k = 2$ it is the binomial.
Each count is binomial, $E[N_i] = np_i$, and different counts are negatively correlated: $\operatorname{Cov}(N_i, N_j) = -np_ip_j$.
Merging and conditioning keep the multinomial form; a Poisson total makes the counts independent Poisson.
The Pearson statistic $\sum (N_i - np_i)^2/(np_i)$ is approximately $\chi^2_{k-1}$.

🧪 Practice

A survey question has answers A, B, C with probabilities 0.5, 0.3, 0.2. In 10 responses, what is the probability of exactly 5 A, 3 B and 2 C? What is the probability of at least one C?

In the setting of exercise 1 with $n = 200$, find the mean and SD of $N_A - N_B$. Hint: use $\operatorname{Var}(N_A) + \operatorname{Var}(N_B) - 2\operatorname{Cov}(N_A, N_B)$.

Show that if $X_1 \sim$ Poisson($\lambda_1$) and $X_2 \sim$ Poisson($\lambda_2$) are independent, then $X_1$ given $X_1 + X_2 = m$ is Binomial($m$, $\lambda_1/(\lambda_1 + \lambda_2)$).

Interview: An A/B/C test splits traffic 1/3 each, and you observe 3,412, 3,300 and 3,288 users. How would you check that the split is working? Hint: this is a sample ratio mismatch check; compare the counts with a multinomial expectation using a chi-square goodness-of-fit test with 2 degrees of freedom.

Bivariate and Multivariate Normal

Why it matters. The multivariate normal is to several variables what the normal is to one: the default model for continuous measurements that are influenced by many small factors, the limit of the multivariate central limit theorem, and the distribution for which everything can be computed in closed form. Linear regression, Kalman filters, Gaussian processes, linear discriminant analysis and portfolio theory all rest on it. Its conditional distributions also explain regression to the mean.

Bivariate normal. $(X, Y)$ is bivariate normal with means $\mu_X, \mu_Y$, standard deviations $\sigma_X, \sigma_Y$ and correlation $\rho$ ($\lvert \rho \rvert < 1$) if its density is

$$f(x, y) = \frac{1}{2\pi \sigma_X \sigma_Y \sqrt{1 - \rho^2}} \exp!\left(-\frac{1}{2(1 - \rho^2)}\left[z_x^2 - 2\rho z_x z_y + z_y^2\right]\right), \quad z_x = \frac{x - \mu_X}{\sigma_X},; z_y = \frac{y - \mu_Y}{\sigma_Y}.$$

Multivariate normal. A random vector $\mathbf{X} \in \mathbb{R}^d$ is $N_d(\boldsymbol{\mu}, \Sigma)$ if

$$f(\mathbf{x}) = \frac{1}{(2\pi)^{d/2} \lvert \Sigma \rvert^{1/2}} \exp!\left(-\frac12 (\mathbf{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})\right).$$

The exponent is minus one half the squared Mahalanobis distance (Section 6.3). An equivalent definition that also covers singular $\Sigma$: $\mathbf{X}$ is multivariate normal if every linear combination $\mathbf{a}^\top \mathbf{X}$ is univariate normal.

Shape: elliptical contours. The density is constant on ellipses (ellipsoids in $d$ dimensions) centered at $\boldsymbol{\mu}$, with axes along the eigenvectors of $\Sigma$ and half-lengths proportional to the square roots of the eigenvalues. Positive $\rho$ tilts the ellipses up and to the right.

  rho = 0                 rho = 0.8                     rho = -0.8
     .-''-.                        .-'                  '-.
   /        \                   .-'  /                  \  '-.
  |    +     |               .-'  + /                    \  +  '-.
   \        /               /   .-'                       '-.   \
     '-..-'                '-.-'                             '-.-'
  circles (equal SDs)      tilted along y = x           tilted along y = -x

The key properties.

Property	Statement
Marginals	each $X_i$ is normal; any subvector is multivariate normal
Linear maps	$A\mathbf{X} + \mathbf{b} \sim N(A\boldsymbol{\mu} + \mathbf{b},, A\Sigma A^\top)$
Uncorrelated = independent	components with zero covariance are independent
Conditionals	normal, with a mean linear in the conditioning values and a constant variance
Mahalanobis distance	$(\mathbf{X} - \boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{X} - \boldsymbol{\mu}) \sim \chi^2_d$
Simulation	$\boldsymbol{\mu} + L\mathbf{Z}$ with $\Sigma = LL^\top$ and $\mathbf{Z}$ i.i.d. $N(0, 1)$

Remember the warning of Section 6.3: normal marginals do not make a vector jointly normal.

Conditional distribution (bivariate). Given $X = x$,

$$Y \mid X = x ;\sim; N!\left(\mu_Y + \rho,\frac{\sigma_Y}{\sigma_X}(x - \mu_X),;; \sigma_Y^2 (1 - \rho^2)\right).$$

The conditional mean is a straight line with slope $\rho\sigma_Y/\sigma_X$: the regression line. The conditional variance does not depend on $x$ and is reduced by the factor $1 - \rho^2$. In general dimension, partitioning $\mathbf{X} = (\mathbf{X}_1, \mathbf{X}_2)$:

$$\mathbf{X}1 \mid \mathbf{X}2 = \mathbf{x}2 ;\sim; N\big(\boldsymbol{\mu}1 + \Sigma{12}\Sigma{22}^{-1}(\mathbf{x}2 - \boldsymbol{\mu}2),;; \Sigma{11} - \Sigma{12}\Sigma{22}^{-1}\Sigma{21}\big).$$

Regression to the mean. In standardized units, the conditional mean is $E[Z_Y \mid Z_X = z] = \rho z$. With $\lvert \rho \rvert < 1$, the prediction is always closer to the mean than the observed value. Children of very tall parents are tall, but on average less extreme than their parents; the best-performing salespeople this year will on average do less well next year; patients enrolled in a study because their blood pressure was very high will on average show lower readings at the next visit, even with no treatment. No causal force is involved: an extreme observation is partly signal and partly luck, and the luck does not repeat.

import numpy as np

rng = np.random.default_rng(78)
n = 2_000_000

# Heights of fathers and sons: mean 175, SD 7, correlation 0.5
mu, sd, rho = 175.0, 7.0, 0.5
Sigma = np.array([[sd**2, rho * sd**2], [rho * sd**2, sd**2]])
L = np.linalg.cholesky(Sigma)
X = mu + rng.standard_normal((n, 2)) @ L.T            # columns: father, son
father, son = X[:, 0], X[:, 1]

# Conditional distribution of the son's height given a 190 cm father
near = np.abs(father - 190) < 0.5
print(f"fathers near 190 cm: {near.sum()}")
print(f"E[son | father = 190]: simulated {son[near].mean():.2f}, formula {mu + rho * (190 - mu):.2f}")
print(f"SD[son | father = 190]: simulated {son[near].std():.2f}, formula {sd * np.sqrt(1 - rho**2):.2f}")

# Regression to the mean works in both directions
near_son = np.abs(son - 190) < 0.5
print(f"E[father | son = 190] = {father[near_son].mean():.2f}: tall sons have less extreme fathers too")

# Mahalanobis distance squared ~ chi-square with 2 df: P(D^2 <= -2 ln(1 - p)) = p
d = X - mu
D2 = np.einsum("ij,jk,ik->i", d, np.linalg.inv(Sigma), d)
for p in [0.5, 0.9, 0.95]:
    r2 = -2 * np.log(1 - p)                           # chi2_2 quantile in closed form
    print(f"ellipse with D^2 <= {r2:.3f} contains {np.mean(D2 <= r2):.4f} of the points (target {p})")

fathers near 190 cm: 11454
E[son | father = 190]: simulated 182.53, formula 182.50
SD[son | father = 190]: simulated 6.06, formula 6.06
E[father | son = 190] = 182.61: tall sons have less extreme fathers too
ellipse with D^2 <= 1.386 contains 0.5000 of the points (target 0.5)
ellipse with D^2 <= 4.605 contains 0.8997 of the points (target 0.9)
ellipse with D^2 <= 5.991 contains 0.9497 of the points (target 0.95)

Key takeaways

$N_d(\boldsymbol{\mu}, \Sigma)$ is determined by its mean vector and covariance matrix; its density contours are ellipses defined by the Mahalanobis distance.
Marginals, linear combinations and conditionals of a multivariate normal are normal; for jointly normal variables, uncorrelated means independent.
$E[Y \mid X = x] = \mu_Y + \rho\frac{\sigma_Y}{\sigma_X}(x - \mu_X)$ is linear and $\operatorname{Var}(Y \mid X) = \sigma_Y^2(1 - \rho^2)$ is constant.
Regression to the mean follows from $\lvert \rho \rvert < 1$: extreme values are followed, on average, by less extreme ones, without any causal mechanism.

🧪 Practice

Exam scores in math and physics are bivariate normal with means 70 and 65, SDs 10 and 12, and correlation 0.6. Find the distribution of the physics score of a student who scored 85 in math, and the probability that this student scores above 80 in physics.

In exercise 1, find the distribution of the total score and of the difference (math minus physics).

$(X_1, X_2, X_3)$ is multivariate normal with $\Sigma_{12} = 0$. Are $X_1$ and $X_2$ independent? Are they conditionally independent given $X_3$? Hint: use the conditional covariance formula $\Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$ with $X_3$ as the conditioning block.

Interview: A new training program is given to the 50 worst-performing employees, and their scores improve significantly the following quarter. The manager wants to roll it out company-wide. What would you say? Hint: the bottom 50 were selected partly for bad luck; regression to the mean predicts improvement with no program at all; a randomized control group is needed (Section 7.2).

Order Statistics

Why it matters. Many quantities of interest are not sums but extremes or ranks: the largest flood in 100 years, the slowest of 50 parallel tasks (which determines when the job finishes), the lifetime of a system that fails when its first component fails, the median of a sample. These are order statistics. Their distributions explain why tail latency grows with fan-out, how sample medians and percentiles behave, and how to estimate the size of a population from serial numbers.

Definition. Sort a sample $X_1, \dots, X_n$ in increasing order:

$$X_{(1)} \le X_{(2)} \le \dots \le X_{(n)}.$$

$X_{(k)}$ is the $k$-th order statistic: $X_{(1)}$ is the minimum, $X_{(n)}$ the maximum, and the sample median is the middle one (or the average of the two middle ones). The order statistics are dependent even when the sample is i.i.d., since $X_{(1)} \le X_{(2)}$ always.

Maximum and minimum. For an i.i.d. sample with CDF $F$:

$$P(X_{(n)} \le x) = P(\text{all } X_i \le x) = F(x)^n, \qquad P(X_{(1)} > x) = P(\text{all } X_i > x) = \big(1 - F(x)\big)^n.$$

So the maximum has CDF $F^n$ and the minimum has CDF $1 - (1 - F)^n$. Two immediate consequences:

Minimum of exponentials: if $X_i \sim$ Exponential($\lambda_i$) are independent, $X_{(1)} \sim$ Exponential($\lambda_1 + \dots + \lambda_n$). A series system of components fails at the rate of the sum of the failure rates.
Tail latency at scale: if each of $n$ parallel calls is slower than its 99th percentile with probability 0.01, the probability that the slowest of $n = 100$ calls is that slow is $1 - 0.99^{100} \approx 0.63$.

The $k$-th order statistic. For a continuous distribution with density $f$, $X_{(k)}$ falls near $x$ when $k - 1$ observations fall below $x$, one near $x$, and $n - k$ above:

$$f_{X_{(k)}}(x) = \frac{n!}{(k - 1)!,(n - k)!}, F(x)^{k - 1}, f(x), \big(1 - F(x)\big)^{n - k}.$$

     k - 1 points below          one point at x           n - k points above
  ----x---x--x------x---------------- * --------x----x-----x---x-----x-----> 
      probability F(x)^(k-1)       f(x) dx           (1 - F(x))^(n-k)
      times the number of ways to choose which observation plays which role

Uniform order statistics are beta. For a Uniform(0, 1) sample, the formula gives $X_{(k)} \sim$ Beta($k$, $n - k + 1$), so

$$E[X_{(k)}] = \frac{k}{n + 1}.$$

The $n$ points split $[0, 1]$ into $n + 1$ gaps of equal expected length $1/(n + 1)$. Combined with the probability integral transform, this says that $F(X_{(k)})$ is Beta($k$, $n - k + 1$) for any continuous $F$, which is the basis of plotting positions such as $k/(n + 1)$ in Q-Q plots (Section 5.3) and of distribution-free confidence intervals for quantiles.

Example: the German tank problem. Serial numbers $1, \dots, N$ are observed on $n$ captured tanks, drawn without replacement. The sample maximum $M$ underestimates $N$; the gap argument suggests adding the average gap, giving the unbiased estimator

$$\hat N = M\Big(1 + \frac1n\Big) - 1.$$

In World War II, estimates of this kind were far more accurate than intelligence reports.

import numpy as np

rng = np.random.default_rng(79)

# Uniform order statistics: X_(k) ~ Beta(k, n - k + 1), mean k / (n + 1)
n = 9
s = np.sort(rng.random((500_000, n)), axis=1)
for k in [1, 3, 5, 9]:
    a, b = k, n - k + 1
    print(f"X_({k}): mean {s[:, k - 1].mean():.4f} (k/(n+1) = {k / (n + 1):.4f}), "
          f"var {s[:, k - 1].var():.5f} (beta: {a * b / ((a + b)**2 * (a + b + 1)):.5f})")

# Minimum of exponentials with rates 1, 2, 3 is Exponential(6)
t = rng.exponential(1 / np.array([1, 2, 3]), (500_000, 3)).min(axis=1)
print(f"min of Exp(1), Exp(2), Exp(3): mean {t.mean():.4f} (1/6 = {1 / 6:.4f})")

# Tail latency: slowest of n parallel calls, each above its own p99 with probability 0.01
for fan_out in [1, 10, 100]:
    print(f"fan-out {fan_out:3d}: P(slowest call exceeds the single-call p99) = {1 - 0.99**fan_out:.3f}")

# German tank problem: N = 1000 tanks, n = 15 serial numbers observed
N, k, reps = 1000, 15, 100_000
serials = np.array([rng.choice(N, k, replace=False) + 1 for _ in range(reps)])   # one sample per row
M = serials.max(axis=1)
estimators = {"sample max": M,
              "2 x mean - 1": 2 * serials.mean(axis=1) - 1,
              "M(1 + 1/n) - 1": M * (1 + 1 / k) - 1}
for name, est in estimators.items():
    print(f"{name:15s}: mean {est.mean():7.1f}, SD {est.std():5.1f}")

X_(1): mean 0.1002 (k/(n+1) = 0.1000), var 0.00822 (beta: 0.00818)
X_(3): mean 0.3000 (k/(n+1) = 0.3000), var 0.01911 (beta: 0.01909)
X_(5): mean 0.5001 (k/(n+1) = 0.5000), var 0.02277 (beta: 0.02273)
X_(9): mean 0.9001 (k/(n+1) = 0.9000), var 0.00817 (beta: 0.00818)
min of Exp(1), Exp(2), Exp(3): mean 0.1668 (1/6 = 0.1667)
fan-out   1: P(slowest call exceeds the single-call p99) = 0.010
fan-out  10: P(slowest call exceeds the single-call p99) = 0.096
fan-out 100: P(slowest call exceeds the single-call p99) = 0.634
sample max     : mean   938.5, SD  58.3
2 x mean - 1   : mean  1000.1, SD 147.9
M(1 + 1/n) - 1 : mean  1000.0, SD  62.2

Both $2\bar X - 1$ and the corrected maximum are unbiased for $N = 1000$, but the one based on the maximum has a much smaller spread: the maximum carries almost all of the information about $N$ (it is a sufficient statistic, Section 8.1).

Key takeaways

The maximum of an i.i.d. sample has CDF $F^n$; the minimum has survival function $(1 - F)^n$.
The $k$-th order statistic has density $\frac{n!}{(k-1)!(n-k)!} F^{k-1} f (1 - F)^{n-k}$; for uniforms, $X_{(k)} \sim$ Beta($k$, $n - k + 1$) with mean $k/(n + 1)$.
The minimum of independent exponentials is exponential with the summed rate; extremes of many parallel tasks make rare delays common.
Order statistics give estimators such as the sample median and the German tank estimator $M(1 + 1/n) - 1$.

🧪 Practice

Five independent light bulbs have Exponential lifetimes with mean 1,000 hours. Find the distribution and mean of the time until the first failure, and the CDF of the time until the last failure.

For a sample of 3 Uniform(0, 1) values, find the density and mean of the median $X_{(2)}$.

Show that the range $X_{(n)} - X_{(1)}$ of $n$ uniforms has mean $(n - 1)/(n + 1)$.

Interview: A service fans out each user request to 200 backend shards and waits for all of them. Each shard's latency has a median of 5 ms and a 99.9th percentile of 50 ms. What latency should a user expect at the median, and what can you do about it? Hint: the median user request waits for the maximum of 200 latencies; compute $P(\max \le x) = F(x)^{200}$ near the 99.9th percentile, then think about hedged requests and tail-tolerant designs.

6.5 Inequalities and Limit Theorems

Inequalities give guaranteed bounds on probabilities and expectations with very little information about a distribution, and limit theorems describe what happens to sums and averages as the number of terms grows; together they justify almost every large-sample method in statistics.

Markov's Inequality

Why it matters. Often we know only the mean of a quantity, such as the average response time or the average loss, and want a guarantee about how often it can be large. Markov's inequality gives such a guarantee for any non-negative random variable, with no other assumption. The bound is usually loose, but it is the root from which Chebyshev's inequality, Chernoff bounds and many proofs in probability and computer science (analysis of randomized algorithms, hashing, sketching) grow.

Statement. If $X \ge 0$ and $a > 0$, then

$$P(X \ge a) \le \frac{E[X]}{a}.$$

Equivalently, $P(X \ge k,E[X]) \le 1/k$: a non-negative variable exceeds $k$ times its mean with probability at most $1/k$.

Proof in one line. For $X \ge 0$, $a,\mathbf{1}{X \ge a} \le X$ always (if $X \ge a$ the left side is $a \le X$; otherwise it is $0 \le X$). Take expectations: $a,P(X \ge a) \le E[X]$.

     X
     ^                     /  X
     |                   /
     |                 /
   a +- - - - - - -  +========  a * 1{X >= a}  (step function, always below the line X)
     |             / :
     |           /   :
     |         /     :
   0 +=======/=======+-------------> outcome, sorted by X
          the step function never exceeds X, so its mean a P(X >= a) never exceeds E[X]

Analogy. If the average income in a town is 40,000, then at most a quarter of the people can earn 160,000 or more: otherwise those people alone would push the average above 40,000 (incomes cannot be negative to compensate).

Tightness. The bound cannot be improved without more information: if $X$ equals $a$ with probability $E[X]/a$ and 0 otherwise, equality holds. For typical distributions it is very loose, because it uses only the mean.

Non-negativity is essential. For a variable that can be negative, $E[X]$ can be small or zero while large values remain likely; the bound then fails. Apply it to $\lvert X \rvert$, $X^2$ or $e^{tX}$ instead.

Generalizations. Applying Markov to an increasing function of $\lvert X \rvert$ gives a whole family of bounds:

Apply Markov to	Result	Name
$(X - \mu)^2$	$P(\lvert X - \mu \rvert \ge a) \le \sigma^2/a^2$	Chebyshev (next topic)
$\lvert X \rvert^p$	$P(\lvert X \rvert \ge a) \le E\lvert X \rvert^p / a^p$	moment bounds
$e^{tX}$, $t > 0$	$P(X \ge a) \le e^{-ta} M_X(t)$, then minimize over $t$	Chernoff bound: exponentially small tails for sums

import numpy as np

rng = np.random.default_rng(80)
n = 1_000_000
samples = {
    "Exponential(mean 1)": rng.exponential(1, n),
    "Uniform(0, 2)":       rng.uniform(0, 2, n),
    "Pareto, alpha = 1.5": rng.pareto(1.5, n) * (1 / 3) + 1 / 3,   # shifted and scaled to mean 1
    "0 or 4 (extremal)":   np.where(rng.random(n) < 0.25, 4.0, 0.0),
}
print(f"{'distribution':22s} {'E[X]':>6s} {'P(X>=2)':>8s} {'P(X>=4)':>8s}   Markov bounds: 1/2 and 1/4")
for name, x in samples.items():
    print(f"{name:22s} {x.mean():6.3f} {np.mean(x >= 2):8.4f} {np.mean(x >= 4):8.4f}")

# Chernoff vs Markov for a sum of 100 fair coin flips: P(S >= 75)
from math import comb
exact = sum(comb(100, k) for k in range(75, 101)) / 2**100
t = np.linspace(0.01, 3, 3000)
chernoff = np.min(np.exp(-75 * t) * ((1 + np.exp(t)) / 2) ** 100)   # M_S(t) = ((1 + e^t)/2)^100
print(f"P(S >= 75): exact {exact:.1e}, Markov {50 / 75:.3f}, Chernoff {chernoff:.1e}")

distribution             E[X]  P(X>=2)  P(X>=4)   Markov bounds: 1/2 and 1/4
Exponential(mean 1)     1.002   0.1358   0.0186
Uniform(0, 2)           1.000   0.0000   0.0000
Pareto, alpha = 1.5     0.987   0.0680   0.0237
0 or 4 (extremal)       1.000   0.2500   0.2500
P(S >= 75): exact 2.8e-07, Markov 0.667, Chernoff 2.1e-06

The extremal two-point distribution meets the bound $P(X \ge 4) \le 1/4$ exactly; the others are far below it. For the coin flips, Markov's bound of 0.67 is useless, while the Chernoff bound, which uses the whole MGF, is within an order of magnitude of the truth.

Key takeaways

For $X \ge 0$: $P(X \ge a) \le E[X]/a$; at most a fraction $1/k$ of the mass can lie beyond $k$ times the mean.
The proof is the pointwise inequality $a,\mathbf{1}{X \ge a} \le X$.
The bound is tight in the worst case but loose for typical distributions; it requires non-negativity.
Applying it to $(X - \mu)^2$ or $e^{tX}$ gives Chebyshev and Chernoff bounds.

🧪 Practice

The average waiting time at a call center is 3 minutes. What does Markov's inequality say about the fraction of callers who wait at least 15 minutes? Compute the actual fraction if waiting times are Exponential with mean 3.

A randomized algorithm has expected running time 10 seconds. Show that if you stop it after 20 seconds and restart, the probability of not having finished after $k$ attempts is at most $2^{-k}$.

Show by example why Markov's inequality fails for a random variable that can take negative values.

Interview: A hash table uses 1,000 buckets for 1,000 keys with a good hash function. Bound the probability that a particular bucket contains 10 or more keys, first with Markov's inequality and then with a better tool. Hint: the bucket size is Binomial(1000, 0.001) with mean 1; compare $1/10$ with a Chernoff bound or the Poisson approximation.

Chebyshev's Inequality

Why it matters. Markov's inequality uses only the mean. If we also know the standard deviation, we can bound the probability of being far from the mean in either direction, for any distribution whatsoever. Chebyshev's inequality makes precise the intuition that "most values are within a few standard deviations of the mean", it gives distribution-free sample size guarantees, and it provides the simplest proof of the law of large numbers.

Statement. For any random variable with mean $\mu$ and finite variance $\sigma^2$, and any $k > 0$,

$$P\big(\lvert X - \mu \rvert \ge k\sigma\big) \le \frac{1}{k^2}, \qquad \text{equivalently} \qquad P\big(\lvert X - \mu \rvert \ge a\big) \le \frac{\sigma^2}{a^2}.$$

So at least $1 - 1/k^2$ of the probability lies within $k$ standard deviations of the mean: at least 75% within 2 SDs, at least 89% within 3 SDs.

Proof. Apply Markov's inequality to the non-negative variable $(X - \mu)^2$: $P((X - \mu)^2 \ge a^2) \le E[(X - \mu)^2]/a^2 = \sigma^2/a^2$.

Chebyshev versus reality.

$k$	Chebyshev bound	Normal	Uniform	Exponential (two-sided)
1.5	0.444	0.134	0.134	0.082
2	0.250	0.046	0.000	0.050
3	0.111	0.0027	0.000	0.018

The bound holds for every distribution, so it must cover the worst one; for any specific, well-behaved distribution it is conservative. The empirical rule (68-95-99.7) is far sharper, but only for approximately normal data. The bound is tight for a three-point distribution that puts mass $1/(2k^2)$ at each of $\mu \pm k\sigma$ and the rest at $\mu$.

One-sided version (Cantelli). For a bound in one direction only,

$$P(X - \mu \ge k\sigma) \le \frac{1}{1 + k^2},$$

which improves on the two-sided bound: for $k = 2$ it gives 0.2 rather than 0.25.

Sample sizes without distributional assumptions. The sample mean of $n$ i.i.d. observations has variance $\sigma^2/n$, so Chebyshev gives

$$P\big(\lvert \bar X_n - \mu \rvert \ge \varepsilon\big) \le \frac{\sigma^2}{n\varepsilon^2}.$$

To guarantee that this probability is at most $\delta$, it suffices that $n \ge \sigma^2/(\varepsilon^2\delta)$. This is the weak law of large numbers in quantitative form. The normal approximation (via the CLT) gives much smaller sample sizes, but only approximately.

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(81)
n = 2_000_000
dists = {
    "normal":      rng.standard_normal(n),
    "uniform":     rng.uniform(-np.sqrt(3), np.sqrt(3), n),       # mean 0, SD 1
    "exponential": rng.exponential(1, n) - 1,                      # mean 0, SD 1
    "t (3 df)":    rng.standard_t(3, n) / np.sqrt(3),              # mean 0, SD 1, heavy tails
}
print(f"{'k':>4s} {'bound':>7s} " + " ".join(f"{name:>12s}" for name in dists))
for k in [1.5, 2, 3, 5]:
    row = " ".join(f"{np.mean(np.abs(x) >= k):12.4f}" for x in dists.values())
    print(f"{k:4.1f} {1 / k**2:7.4f} {row}")

# Sample size for P(|mean - mu| >= 0.1 sigma) <= 0.05: Chebyshev vs normal approximation
eps, delta = 0.1, 0.05
n_cheb = 1 / (eps**2 * delta)
n_clt = (NormalDist().inv_cdf(1 - delta / 2) / eps) ** 2
print(f"required n: Chebyshev {n_cheb:.0f}, normal approximation {n_clt:.0f}")

   k   bound       normal      uniform  exponential     t (3 df)
 1.5  0.4444       0.1337       0.1341       0.0821       0.0807
 2.0  0.2500       0.0453       0.0000       0.0499       0.0406
 3.0  0.1111       0.0026       0.0000       0.0184       0.0139
 5.0  0.0400       0.0000       0.0000       0.0025       0.0032
required n: Chebyshev 2000, normal approximation 384

Even the heavy-tailed $t$ distribution stays well inside the bound, as it must; the price of Chebyshev's generality is a sample size more than five times larger than the normal approximation requires.

Key takeaways

$P(\lvert X - \mu \rvert \ge k\sigma) \le 1/k^2$ for every distribution with finite variance: at least 75% within 2 SDs, 89% within 3 SDs.
It is Markov's inequality applied to $(X - \mu)^2$; it is tight only for an extreme three-point distribution.
Cantelli's one-sided bound is $1/(1 + k^2)$.
Applied to the sample mean, it proves the weak law of large numbers and gives conservative, assumption-free sample sizes $n \ge \sigma^2/(\varepsilon^2\delta)$.

🧪 Practice

A machine fills bottles with mean 500 ml and SD 4 ml. Without assuming normality, bound the fraction of bottles outside 490 to 510 ml. What is the fraction if the fill is normal?

How many fair coin tosses guarantee, by Chebyshev, that the proportion of heads is within 0.01 of 1/2 with probability at least 0.99?

Construct a distribution with mean 0 and variance 1 for which $P(\lvert X \rvert \ge 2) = 1/4$ exactly.

Interview: A monitoring system flags any metric more than 3 standard deviations from its mean. A colleague claims this gives a 0.3% false alarm rate. Under what assumptions is that right, and what can you guarantee without them? Hint: 0.27% is the normal tail; Chebyshev guarantees only 11.1%; heavy-tailed or skewed metrics fall in between, and the mean and SD themselves must be estimated.

Jensen's Inequality

Why it matters. Is the average of a function the same as the function of the average? Almost never, and Jensen's inequality tells you in which direction the error goes. It explains why the average speed on a round trip is not the average of the two speeds, why volatility lowers long-run investment growth, why the sample standard deviation underestimates $\sigma$, why the log-likelihood is used in EM, and why the Kullback-Leibler divergence is never negative.

Convex functions. A function $\varphi$ is convex if the chord between any two points of its graph lies above the graph: $\varphi(\lambda a + (1 - \lambda)b) \le \lambda\varphi(a) + (1 - \lambda)\varphi(b)$. Examples: $x^2$, $e^x$, $\lvert x \rvert$, $1/x$ for $x > 0$, $-\log x$. A function is concave if $-\varphi$ is convex: $\log x$, $\sqrt x$.

Statement. If $\varphi$ is convex and $E[X]$ and $E[\varphi(X)]$ exist, then

$$\varphi\big(E[X]\big) \le E\big[\varphi(X)\big].$$

For concave $\varphi$ the inequality reverses. Equality holds when $X$ is constant or $\varphi$ is linear on the range of $X$.

  phi(x)
    ^
    |  *                               *      phi is convex (bowl-shaped); X = a or b, each with probability 1/2
    |   \~ ~ ~ ~ ~ ~ ~ ~o~ ~ ~ ~ ~ ~ ~/       ~ ~ ~ : chord from (a, phi(a)) to (b, phi(b))
    |    \              |            /
    |     \             |           /         o = chord at E[X]  = (phi(a) + phi(b)) / 2 = E[phi(X)]
    |      '.           |         .'          + = curve at E[X]  = phi(E[X])
    |        '-.        |      .-'
    |           '--.____+__.--'               o lies above +, so E[phi(X)] >= phi(E[X])
    +----+--------------+--------------+----> x
         a            E[X]             b

Proof sketch. A convex function lies above each of its tangent lines. Take the tangent at $\mu = E[X]$: $\varphi(x) \ge \varphi(\mu) + c(x - \mu)$ for some slope $c$. Plug in $X$ and take expectations; the linear term has mean zero, leaving $E[\varphi(X)] \ge \varphi(\mu)$.

Consequences.

Choose $\varphi$	Inequality	Meaning
$x^2$	$E[X]^2 \le E[X^2]$	variance is non-negative
$1/x$ on $x > 0$	$1/E[X] \le E[1/X]$	averaging rates and times are different things
$\log x$ (concave)	$E[\log X] \le \log E[X]$	geometric mean $\le$ arithmetic mean; volatility drag
$\sqrt x$ (concave), applied to $S^2$	$E[S] \le \sqrt{E[S^2]} = \sigma$	the sample SD is biased downward
$-\log x$, applied to $q(X)/p(X)$	$D_{KL}(p ,\Vert, q) \ge 0$	Gibbs' inequality, foundation of information theory

Example: average speed. You drive 60 km at 30 km/h and return at 90 km/h. The average speed is not $60$ km/h but total distance over total time: $120 / (2 + 2/3) = 45$ km/h. Time is distance times $1/\text{speed}$, and $1/x$ is convex, so averaging speeds overstates the true average speed.

Example: volatility drag. An investment gains 50% or loses 50% with equal probability each year. The expected one-year growth factor is $E[R] = 1$, but the long-run growth is governed by $E[\log R] = \frac12(\log 1.5 + \log 0.5) = \log\sqrt{0.75} < 0$: the typical investor loses money, about 13% per year.

import numpy as np

rng = np.random.default_rng(82)

# Average speed: 1/x is convex, so the mean of speeds overstates the true average speed
speeds = np.array([30, 90])
print(f"mean of speeds = {speeds.mean():.1f} km/h, true average = {1 / np.mean(1 / speeds):.1f} km/h (harmonic mean)")

# Volatility drag: E[R] = 1 but E[log R] < 0
R = rng.choice([1.5, 0.5], size=(100_000, 30))        # 30 years of +50% / -50%
wealth = R.prod(axis=1)
print(f"E[R] = {R.mean():.3f}, exp(E[log R]) = {np.exp(np.log(R).mean()):.3f}")
print(f"after 30 years: mean wealth {wealth.mean():.2f}, median wealth {np.median(wealth):.4f}, "
      f"P(wealth < 1) = {np.mean(wealth < 1):.3f}")

# Sample SD is biased low: E[S] < sigma (Jensen with the concave square root)
for n in [2, 5, 30]:
    x = rng.standard_normal((200_000, n))
    s = x.std(axis=1, ddof=1)                         # ddof = 1 makes S^2 unbiased for sigma^2 = 1
    print(f"n = {n:2d}: E[S^2] = {np.mean(s**2):.4f}, E[S] = {s.mean():.4f} (sigma = 1)")

mean of speeds = 60.0 km/h, true average = 45.0 km/h (harmonic mean)
E[R] = 1.000, exp(E[log R]) = 0.866
after 30 years: mean wealth 1.11, median wealth 0.0134, P(wealth < 1) = 0.899
n =  2: E[S^2] = 0.9994, E[S] = 0.7982 (sigma = 1)
n =  5: E[S^2] = 1.0015, E[S] = 0.9407 (sigma = 1)
n = 30: E[S^2] = 1.0007, E[S] = 0.9918 (sigma = 1)

The 30-year simulation shows the gap between the mean and the typical outcome: the mean wealth is pulled up by a handful of extremely lucky paths, while the median investor is left with almost nothing.

Key takeaways

For convex $\varphi$: $\varphi(E[X]) \le E[\varphi(X)]$; for concave $\varphi$ the inequality reverses.
The function of the mean is not the mean of the function; the direction of the error follows from the curvature.
Consequences include $\operatorname{Var}(X) \ge 0$, AM-GM, the downward bias of the sample SD, volatility drag and $D_{KL} \ge 0$.

🧪 Practice

A server handles requests whose sizes have mean 10 MB. Is the mean transfer time at 2 MB/s equal to 5 seconds? Is the mean throughput (size divided by time) of requests with random times equal to size divided by mean time? Use Jensen to decide the direction.

Show that $E[e^X] \ge e^{E[X]}$, and compute both sides for $X \sim N(\mu, \sigma^2)$ using the lognormal mean $e^{\mu + \sigma^2/2}$.

Prove that $\frac{1}{n}\sum \log x_i \le \log\big(\frac1n \sum x_i\big)$ for positive $x_i$, and deduce the AM-GM inequality.

Interview: A model predicts $\log(\text{price})$ with an unbiased error. You exponentiate the predictions to get prices. Are the price predictions biased, and in which direction? Hint: $\exp$ is convex, so $E[\exp(\cdot)] \ge \exp(E[\cdot])$; under normal errors the correction factor is $e^{\sigma^2/2}$ (smearing estimators handle the general case).

Cauchy-Schwarz for Random Variables

Why it matters. The Cauchy-Schwarz inequality is the reason a correlation can never exceed 1 in absolute value. More generally, it bounds the expectation of a product by second moments, which is useful whenever two random quantities interact: bounding a covariance you cannot compute, proving that estimators have finite variance, or deriving the Cramér-Rao lower bound on the variance of estimators (Section 8.2).

Statement. For random variables with finite second moments,

$$\big(E[XY]\big)^2 \le E[X^2], E[Y^2].$$

Equality holds if and only if one variable is a constant multiple of the other with probability 1 ($Y = cX$, or $X = 0$).

Proof. For every real $t$, the quantity $E[(Y - tX)^2] = E[Y^2] - 2tE[XY] + t^2E[X^2]$ is non-negative. A quadratic in $t$ that never goes below zero has a non-positive discriminant: $4E[XY]^2 - 4E[X^2]E[Y^2] \le 0$. Equality means the quadratic touches zero, so $Y - tX = 0$ with probability 1 for some $t$.

Geometric meaning. Random variables with finite second moments form a vector space with the inner product $\langle X, Y \rangle = E[XY]$ and the length $\lVert X \rVert = \sqrt{E[X^2]}$. Cauchy-Schwarz is the familiar $\lvert \langle X, Y \rangle \rvert \le \lVert X \rVert \lVert Y \rVert$, and for centered variables:

Geometry	Probability (centered variables)
inner product $\langle X, Y \rangle$	covariance $\operatorname{Cov}(X, Y)$
squared length $\lVert X \rVert^2$	variance $\operatorname{Var}(X)$
$\cos$ of the angle between them	correlation $\rho(X, Y)$
orthogonal	uncorrelated
projection onto a subspace	best linear prediction (regression)

The covariance inequality. Applying Cauchy-Schwarz to $X - \mu_X$ and $Y - \mu_Y$:

$$\lvert \operatorname{Cov}(X, Y) \rvert \le \sigma_X \sigma_Y \quad\Longleftrightarrow\quad -1 \le \rho \le 1,$$

with $\lvert \rho \rvert = 1$ only for exact linear relationships, as found in Section 6.3.

Other uses.

Bounding unknown covariances. If you know only the variances of two risks, $\operatorname{Var}(X + Y) \le (\sigma_X + \sigma_Y)^2$: the worst case is perfect correlation.
Bounding a mean with an indicator. $E[X\mathbf{1}_A]^2 \le E[X^2],P(A)$: an event of small probability cannot carry much of the mean of a variable with a moderate second moment.
Paley-Zygmund inequality. For $X \ge 0$, $P(X > 0) \ge E[X]^2/E[X^2]$: a lower bound on the probability that a count is positive, used in the second moment method.

import numpy as np

rng = np.random.default_rng(83)
n = 1_000_000

x = rng.exponential(1, n)
cases = {"Y = X^2 (dependent)": x**2,
         "Y = 3X (proportional)": 3 * x,
         "Y independent": rng.exponential(1, n)}
for name, y in cases.items():
    lhs, rhs = np.mean(x * y) ** 2, np.mean(x**2) * np.mean(y**2)
    print(f"{name:22s}: E[XY]^2 = {lhs:8.3f} <= E[X^2]E[Y^2] = {rhs:8.3f}, ratio {lhs / rhs:.3f}")

# Bounding the share of the mean carried by a rare event: E[X 1_A]^2 <= E[X^2] P(A)
A = x > 5                                             # rare event, P(A) = e^-5
print(f"P(A) = {A.mean():.5f}; E[X 1_A] = {np.mean(x * A):.4f} <= sqrt(E[X^2] P(A)) = {np.sqrt(np.mean(x**2) * A.mean()):.4f}")

# Paley-Zygmund: number of heads among 10 coins with p = 0.05
k = rng.binomial(10, 0.05, n)
print(f"P(K > 0) = {np.mean(k > 0):.4f} >= E[K]^2 / E[K^2] = {k.mean()**2 / np.mean(k**2):.4f}")

Y = X^2 (dependent)   : E[XY]^2 =   36.195 <= E[X^2]E[Y^2] =   48.208, ratio 0.751
Y = 3X (proportional) : E[XY]^2 =   36.142 <= E[X^2]E[Y^2] =   36.142, ratio 1.000
Y independent         : E[XY]^2 =    1.001 <= E[X^2]E[Y^2] =    4.018, ratio 0.249
P(A) = 0.00677; E[X 1_A] = 0.0406 <= sqrt(E[X^2] P(A)) = 0.1165
P(K > 0) = 0.4013 >= E[K]^2 / E[K^2] = 0.3449

Only the proportional case reaches ratio 1, as the equality condition predicts.

Key takeaways

$(E[XY])^2 \le E[X^2]E[Y^2]$, with equality only when one variable is a multiple of the other.
The proof: a non-negative quadratic $E[(Y - tX)^2]$ has a non-positive discriminant.
Applied to centered variables it gives $\lvert \operatorname{Cov}(X, Y) \rvert \le \sigma_X\sigma_Y$, that is $\lvert \rho \rvert \le 1$.
Geometrically: covariance is an inner product, variance a squared length, correlation a cosine.

🧪 Practice

Use Cauchy-Schwarz to show that $E\lvert X \rvert \le \sqrt{E[X^2]}$, and compare the two sides for $X \sim N(0, 1)$, where $E\lvert X \rvert = \sqrt{2/\pi}$.

Two assets have return SDs of 20% and 30%. Without knowing their correlation, give the range of possible SDs of an equally weighted portfolio.

Use Paley-Zygmund to show that if $X \sim$ Binomial($n$, $p$) then $P(X > 0) \ge \frac{np}{1 + (n - 1)p}$.

Interview: Explain the connection between linear regression and orthogonal projection in the space of random variables. Hint: the best linear predictor $a + bX$ of $Y$ makes the error orthogonal (uncorrelated) to 1 and to $X$; write those two conditions and solve for $a$ and $b$.

Modes of Convergence

Why it matters. Statements such as "the sample mean converges to the population mean" or "the distribution of the test statistic converges to a normal" are the backbone of large-sample statistics. But a sequence of random variables can converge in several different senses, and the senses are not equivalent. Knowing which mode a theorem guarantees tells you what you can and cannot conclude: whether a single long simulation will settle down, whether an estimator's mean converges, or only its distribution.

The four main modes. Let $X_1, X_2, \dots$ and $X$ be random variables.

Mode	Notation	Definition	Intuition
Almost surely	$X_n \xrightarrow{a.s.} X$	$P(\lim_n X_n = X) = 1$	every path (except a negligible set) eventually settles at $X$
In probability	$X_n \xrightarrow{P} X$	$P(\lvert X_n - X \rvert > \varepsilon) \to 0$ for every $\varepsilon > 0$	for large $n$, a big deviation is unlikely at that $n$
In mean square ($L^2$)	$X_n \xrightarrow{L^2} X$	$E[(X_n - X)^2] \to 0$	the average squared error vanishes
In distribution	$X_n \xrightarrow{d} X$	$F_{X_n}(x) \to F_X(x)$ at every continuity point of $F_X$	the histograms look alike; the variables need not be close

How they relate.

   almost surely  ----+
                      +---->  in probability  ---->  in distribution
   in mean square ----+
   (L^p, p >= 1)

   No other arrows hold in general. Special case: convergence in distribution to a CONSTANT
   implies convergence in probability to that constant.

Counterexamples that separate the modes.

In probability but not almost surely. Let $X_n$ be independent with $P(X_n = 1) = 1/n$ and $X_n = 0$ otherwise. $P(X_n \ne 0) = 1/n \to 0$, so $X_n \to 0$ in probability. But $\sum 1/n = \infty$, so by the second Borel-Cantelli lemma a 1 appears infinitely often on almost every path: no path settles at 0. With $P(X_n = 1) = 1/n^2$ instead, the sum is finite, only finitely many 1s occur, and convergence is almost sure.
In probability but not in mean. Let $X_n = n$ with probability $1/n$ and 0 otherwise. $X_n \to 0$ in probability, but $E[X_n] = 1$ for all $n$: a rare but huge value keeps the mean from converging.
In distribution but not in probability. Let $X \sim N(0, 1)$ and $X_n = -X$ for all $n$. Each $X_n$ has exactly the distribution of $X$, so convergence in distribution holds trivially, yet $\lvert X_n - X \rvert = 2\lvert X \rvert$ never shrinks.

Tools for combining limits.

Continuous mapping theorem: if $X_n \to X$ in any of the three senses (a.s., in probability, in distribution) and $g$ is continuous, then $g(X_n) \to g(X)$ in the same sense.
Slutsky's theorem: if $X_n \xrightarrow{d} X$ and $Y_n \xrightarrow{P} c$ (a constant), then $X_n + Y_n \xrightarrow{d} X + c$ and $X_nY_n \xrightarrow{d} cX$. This is why the $t$ statistic, with the estimated SD in place of the true one, is still asymptotically normal.

import numpy as np

rng = np.random.default_rng(84)
paths, N = 1_000, 20_000
n = np.arange(1, N + 1)

# Independent X_n ~ Bernoulli(1/n): converges in probability, not almost surely
ones = rng.random((paths, N)) < 1 / n                  # each row is one path X_1, ..., X_N
late = ones[:, 1000:].any(axis=1)                     # does a 1 still appear after n = 1000?
print(f"P(X_n = 1) at n = {N:,}: {1 / N:.0e}; fraction of paths with a 1 after n = 1000: {late.mean():.3f}")
print(f"mean number of 1s between n = 1001 and {N:,}: {ones[:, 1000:].sum(axis=1).mean():.2f} "
      f"(sum of 1/n = {np.sum(1 / n[1000:]):.2f})")

# Independent X_n ~ Bernoulli(1/n^2): converges almost surely
ones2 = rng.random((paths, N)) < 1 / n**2
print(f"with P(X_n = 1) = 1/n^2: fraction of paths with a 1 after n = 1000: {ones2[:, 1000:].any(axis=1).mean():.4f}")

# X_n = n with probability 1/n: converges in probability to 0, but E[X_n] = 1
for m in [10, 1000, 100_000]:
    x = np.where(rng.random(10_000_000) < 1 / m, m, 0)
    print(f"n = {m:6d}: P(X_n > 0.5) = {np.mean(x > 0.5):.5f}, E[X_n] = {x.mean():.3f}")

P(X_n = 1) at n = 20,000: 5e-05; fraction of paths with a 1 after n = 1000: 0.962
mean number of 1s between n = 1001 and 20,000: 3.02 (sum of 1/n = 3.00)
with P(X_n = 1) = 1/n^2: fraction of paths with a 1 after n = 1000: 0.0020
n =     10: P(X_n > 0.5) = 0.10004, E[X_n] = 1.000
n =   1000: P(X_n > 0.5) = 0.00099, E[X_n] = 0.995
n = 100000: P(X_n > 0.5) = 0.00001, E[X_n] = 0.970

In the first experiment, the probability of a 1 at any given late time is tiny, yet almost every path still contains 1s after time 1000: the expected count, $\sum_{n=1001}^{20000} 1/n \approx 3.0$, keeps growing like $\log n$ as the horizon extends. That is the difference between "unlikely at each time" and "eventually never".

Key takeaways

Almost sure convergence is about whole paths; convergence in probability is about each fixed large $n$; mean-square convergence is about average squared error; convergence in distribution is only about CDFs.
Almost sure and $L^2$ convergence each imply convergence in probability, which implies convergence in distribution; the converses fail.
Rare large values can break convergence of means while convergence in probability still holds.
Slutsky's theorem and the continuous mapping theorem let you combine limits, for example replacing an unknown $\sigma$ by a consistent estimate.

🧪 Practice

Let $\bar X_n$ be the mean of $n$ i.i.d. variables with variance $\sigma^2$. Show that $\bar X_n \to \mu$ in mean square, and deduce convergence in probability.

Let $X_n \sim$ Uniform$(0, 1/n)$. In which of the four senses does $X_n$ converge, and to what?

Let $X_n \sim N(0, 1 + 1/n)$. Does $X_n$ converge in distribution? Does it converge in probability to anything?

Interview: A colleague says "my estimator is consistent, so with a large enough sample its bias is negligible". Is that implied? Hint: consistency is convergence in probability; the "$X_n = n$ with probability $1/n$" example shows that the mean need not converge, so bias and MSE require a separate argument.

Weak and Strong Laws of Large Numbers

Why it matters. Why does the proportion of heads approach 1/2? Why do casinos and insurance companies, which face pure randomness on every bet or policy, earn predictable profits? Why does a Monte Carlo simulation give the right answer if run long enough? The laws of large numbers are the mathematical answer: the average of many independent draws converges to the expected value. They justify interpreting probability as long-run frequency (Section 3.1) and estimating population means by sample means.

Setup. Let $X_1, X_2, \dots$ be i.i.d. with mean $\mu = E[X_1]$, and $\bar X_n = (X_1 + \dots + X_n)/n$.

Weak law (WLLN). For every $\varepsilon > 0$,

$$P\big(\lvert \bar X_n - \mu \rvert > \varepsilon\big) \to 0 \quad \text{as } n \to \infty, \qquad \text{that is, } \bar X_n \xrightarrow{P} \mu.$$

With finite variance, the proof is Chebyshev's inequality: $P(\lvert \bar X_n - \mu \rvert > \varepsilon) \le \sigma^2/(n\varepsilon^2) \to 0$. (The law holds even without a finite variance, as long as $E\lvert X_1 \rvert < \infty$.)

Strong law (SLLN). If $E\lvert X_1 \rvert < \infty$, then

$$P\Big(\lim_{n \to \infty} \bar X_n = \mu\Big) = 1, \qquad \text{that is, } \bar X_n \xrightarrow{a.s.} \mu.$$

Weak versus strong.

	Weak law	Strong law
Mode	in probability	almost surely
Says	at any fixed large $n$, the average is probably close to $\mu$	the sequence of averages of (almost) every single path converges to $\mu$
Typical proof	Chebyshev's inequality	Borel-Cantelli with fourth moments, or more refined arguments
Practical reading	a large sample gives a good estimate	one long simulation run settles to the right answer

What the laws do not say.

No compensation (gambler's fallacy). After 10 heads in a row, the next toss is still 50/50. The proportion converges because early deviations are diluted by many new tosses, not corrected. In fact the count of heads minus $n/2$ typically grows like $\sqrt n$ (Section 5.5).
No rate. The laws do not say how large $n$ must be; that comes from Chebyshev (conservative) or from the CLT (approximate): the error is of order $\sigma/\sqrt n$.
No finite mean, no law. For the Cauchy distribution, $E\lvert X \rvert = \infty$, and the average of $n$ standard Cauchy variables has the same Cauchy distribution as a single one: averaging does not help at all.

Applications.

Monte Carlo integration: $E[g(X)] \approx \frac1n \sum g(X_i)$, with error of order $1/\sqrt n$ regardless of dimension.
Frequencies estimate probabilities: apply the law to indicators, $\frac1n \sum \mathbf{1}_A(X_i) \to P(A)$; this is also why the empirical CDF converges to the true CDF (and, uniformly, by the Glivenko-Cantelli theorem).
Insurance and casinos: the average result per policy or bet converges to its expectation, so a large book of independent risks has a predictable total per policy. Dependence (a hurricane, a pandemic) breaks this.

import numpy as np

rng = np.random.default_rng(85)
n = 1_000_000

# Running mean of die rolls for 5 independent paths
rolls = rng.integers(1, 7, (5, n))
running = rolls.cumsum(axis=1) / np.arange(1, n + 1)
print("running means of 5 paths (mu = 3.5)")
for m in [10, 100, 10_000, 1_000_000]:
    print(f"  n = {m:9,d}: " + "  ".join(f"{r:.4f}" for r in running[:, m - 1]))

# The proportion converges, the count difference does not: heads - n/2 grows like sqrt(n)
flips = rng.integers(0, 2, (2_000, 10_000))
heads = flips.cumsum(axis=1)
for m in [100, 10_000]:
    print(f"  n = {m:6d}: SD of (proportion - 1/2) = {(heads[:, m - 1] / m - 0.5).std():.4f}, "
          f"SD of (heads - n/2) = {(heads[:, m - 1] - m / 2).std():.1f}")

# Cauchy: no mean, no law of large numbers
c = rng.standard_cauchy((5, n))
running_c = c.cumsum(axis=1) / np.arange(1, n + 1)
print("running means of 5 Cauchy paths")
for m in [100, 10_000, 1_000_000]:
    print(f"  n = {m:9,d}: " + "  ".join(f"{r:+8.3f}" for r in running_c[:, m - 1]))

# Monte Carlo integration: E[exp(-U^2)] for U ~ Uniform(0, 1), exact 0.746824
for m in [100, 10_000, 1_000_000]:
    est = np.exp(-rng.random(m) ** 2).mean()
    print(f"  Monte Carlo with n = {m:9,d}: {est:.6f} (error {abs(est - 0.746824):.1e})")

running means of 5 paths (mu = 3.5)
  n =        10: 3.6000  3.4000  2.7000  3.7000  3.2000
  n =       100: 3.5200  3.6900  3.5500  3.4100  3.7200
  n =    10,000: 3.5213  3.5128  3.4874  3.5121  3.4742
  n = 1,000,000: 3.5006  3.4991  3.5002  3.4988  3.5018
  n =    100: SD of (proportion - 1/2) = 0.0498, SD of (heads - n/2) = 5.0
  n =  10000: SD of (proportion - 1/2) = 0.0051, SD of (heads - n/2) = 50.6
running means of 5 Cauchy paths
  n =       100:   +0.927    -0.340    -0.429    +0.311    -0.965
  n =    10,000:   -0.289    -0.597    -1.026    -0.130    -0.779
  n = 1,000,000:   -0.320    +0.416    +2.133    -0.179    -1.356
  Monte Carlo with n =       100: 0.743852 (error 3.0e-03)
  Monte Carlo with n =    10,000: 0.749700 (error 2.9e-03)
  Monte Carlo with n = 1,000,000: 0.746845 (error 2.1e-05)

The die averages tighten around 3.5 at the rate $1/\sqrt n$, while the Cauchy averages jump around at every scale: a single enormous observation can move the average arbitrarily far, no matter how many observations came before.

Key takeaways

WLLN: $\bar X_n \to \mu$ in probability; SLLN: $\bar X_n \to \mu$ almost surely; both need only $E\lvert X \rvert < \infty$ for i.i.d. samples.
The proportion converges because early deviations are diluted, not compensated; the gambler's fallacy is a misreading of the law.
The laws give no rate; Chebyshev and the CLT supply an error of order $\sigma/\sqrt n$.
Without a finite mean (Cauchy), sample averages do not converge.

🧪 Practice

Use the WLLN to explain why the empirical frequency of an event in i.i.d. trials converges to its probability.

Write a Monte Carlo estimate of $\pi$ by sampling points in the unit square, and estimate how many points are needed for an error below 0.001 with high probability.

Let $X_i$ be i.i.d. with $P(X_i = 2^k) = 2^{-k}$, $k = 1, 2, \dots$ (the St. Petersburg game). Does the law of large numbers apply to $\bar X_n$? What happens to $\bar X_n$ in a simulation?

Interview: An insurer argues that writing many more policies in the same coastal city will make its results more predictable, by the law of large numbers. What is wrong with the argument? Hint: the law needs independence (or at least weak dependence); with a common shock such as a hurricane, the equicorrelated variance formula of Section 6.4 shows that $\operatorname{Var}(\bar X)$ does not go to 0.

Central Limit Theorem

Why it matters. The law of large numbers says that the sample mean approaches $\mu$; the central limit theorem (CLT) describes the fluctuations around $\mu$, and says that they are approximately normal, whatever the distribution of the individual observations. This is why the normal distribution is everywhere, why confidence intervals of the form $\bar x \pm 1.96, s/\sqrt n$ work for non-normal data, and why most test statistics have known approximate distributions. It is arguably the most important theorem in statistics.

Statement (Lindeberg-Lévy). Let $X_1, X_2, \dots$ be i.i.d. with mean $\mu$ and finite variance $\sigma^2 > 0$. Then

$$Z_n = \frac{\bar X_n - \mu}{\sigma/\sqrt n} = \frac{X_1 + \dots + X_n - n\mu}{\sigma\sqrt n} \xrightarrow{d} N(0, 1).$$

In practice: for large $n$, $\bar X_n \approx N(\mu, \sigma^2/n)$ and $X_1 + \dots + X_n \approx N(n\mu, n\sigma^2)$.

  distribution of the mean of n exponential observations (strongly right-skewed parent)

  n = 1        n = 5            n = 30                n = 200
  #            .##.             .###.                  .#.
  ##          .####.           .#####.                .###.
  ###.       .######..        .#######.              .#####.
  #####....  ##########...   ###########.           .#######.
  -------    -----------     ------------           ---------
  skewed     still skewed    nearly symmetric       normal, narrow (SD sigma/sqrt n)

Why it happens (sketch). The MGF of the standardized sum is $\big(M_Y(t/\sqrt n)\big)^n$ with $Y = (X - \mu)/\sigma$. A second-order Taylor expansion gives $M_Y(s) = 1 + s^2/2 + o(s^2)$, because $E[Y] = 0$ and $E[Y^2] = 1$, so

$$\Big(1 + \frac{t^2}{2n} + o(1/n)\Big)^n \to e^{t^2/2},$$

the MGF of $N(0, 1)$. Only the mean and variance survive; every other feature of the parent distribution is washed out. (The rigorous proof uses characteristic functions, which always exist.)

How large must $n$ be? There is no universal answer; "$n \ge 30$" is a rule of thumb, not a theorem. The Berry-Esseen theorem bounds the error:

$$\sup_x \lvert P(Z_n \le x) - \Phi(x) \rvert \le \frac{C, E\lvert X - \mu \rvert^3}{\sigma^3 \sqrt n}, \qquad C < 0.48.$$

The approximation error shrinks like $1/\sqrt n$ and grows with the skewness of the parent. Symmetric, light-tailed data (uniform) are nearly normal after a handful of terms; skewed data (exponential, rare-event Bernoulli, income) may need hundreds or thousands, and the tails converge slowest of all.

Variants and extensions.

Version	What it adds
Binomial (de Moivre-Laplace)	Binomial($n$, $p$) $\approx N(np, np(1 - p))$ (Section 5.3)
Lindeberg / Lyapunov	independent but not identically distributed terms, as long as no single term dominates
Multivariate CLT	$\sqrt n(\bar{\mathbf{X}}_n - \boldsymbol{\mu}) \xrightarrow{d} N_d(\mathbf{0}, \Sigma)$
Dependent data	weakly dependent sequences (time series, Markov chains) with an adjusted variance
Infinite variance	sums converge to stable distributions, not the normal (Pareto tails with $\alpha < 2$)

Common misconceptions. The CLT is about the distribution of the mean (or sum), not of the data: a large sample of skewed incomes stays skewed. It also does not make individual observations normal, and it does not apply to the maximum or to other statistics that are not averages.

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(86)
z = NormalDist()
reps = 100_000

def tail_check(draw, mu, sigma, label):
    """Compare P(Z_n <= -1.645) and P(Z_n >= 1.645) with the normal value 0.05."""
    print(label)
    for n in [1, 5, 30, 200]:
        means = draw((reps, n)).mean(axis=1)
        zn = (means - mu) / (sigma / np.sqrt(n))
        skew = np.mean(zn**3)                         # skewness of the standardized mean
        print(f"  n = {n:3d}: P(Z <= -1.645) = {np.mean(zn <= -1.645):.4f}, "
              f"P(Z >= 1.645) = {np.mean(zn >= 1.645):.4f}, skewness {skew:+.3f}")

tail_check(lambda s: rng.uniform(0, 1, s), 0.5, np.sqrt(1 / 12), "Uniform parent (symmetric)")
tail_check(lambda s: rng.exponential(1, s), 1.0, 1.0, "Exponential parent (skewness 2)")

# 95% interval for a rare-event rate: coverage depends on n * p, not on n alone
p = 0.01
for n in [100, 1000, 10_000]:
    phat = rng.binomial(n, p, reps) / n
    half = z.inv_cdf(0.975) * np.sqrt(phat * (1 - phat) / n)
    cover = np.mean(np.abs(phat - p) <= half)
    print(f"p = 0.01, n = {n:6d} (np = {n * p:5.0f}): coverage of the 95% normal interval = {cover:.3f}")

Uniform parent (symmetric)
  n =   1: P(Z <= -1.645) = 0.0251, P(Z >= 1.645) = 0.0250, skewness -0.001
  n =   5: P(Z <= -1.645) = 0.0503, P(Z >= 1.645) = 0.0504, skewness -0.005
  n =  30: P(Z <= -1.645) = 0.0502, P(Z >= 1.645) = 0.0497, skewness -0.010
  n = 200: P(Z <= -1.645) = 0.0501, P(Z >= 1.645) = 0.0501, skewness +0.004
Exponential parent (skewness 2)
  n =   1: P(Z <= -1.645) = 0.0000, P(Z >= 1.645) = 0.0712, skewness +2.029
  n =   5: P(Z <= -1.645) = 0.0112, P(Z >= 1.645) = 0.0662, skewness +0.896
  n =  30: P(Z <= -1.645) = 0.0373, P(Z >= 1.645) = 0.0593, skewness +0.374
  n = 200: P(Z <= -1.645) = 0.0458, P(Z >= 1.645) = 0.0535, skewness +0.123
p = 0.01, n =    100 (np =     1): coverage of the 95% normal interval = 0.633
p = 0.01, n =   1000 (np =    10): coverage of the 95% normal interval = 0.927
p = 0.01, n =  10000 (np =   100): coverage of the 95% normal interval = 0.948

For the uniform parent the normal tail probabilities are essentially right at $n = 5$. For the exponential parent, the two tails are still visibly unequal at $n = 30$, and the skewness of the mean decays only like $2/\sqrt n$. For a rare event, the normal interval undercovers until the expected count $np$ is large: what matters is the number of events, not the number of trials.

Key takeaways

For i.i.d. data with finite variance, $(\bar X_n - \mu)/(\sigma/\sqrt n) \to N(0, 1)$ in distribution, whatever the parent distribution.
The error of the approximation shrinks like $1/\sqrt n$ and grows with skewness; $n \ge 30$ is a rule of thumb that fails for skewed data and rare events.
The CLT describes means and sums, not individual observations or extremes.
Extensions cover non-identical terms, vectors and weakly dependent data; infinite variance breaks it.

🧪 Practice

The weight of a parcel has mean 2 kg and SD 1.5 kg. A truck carries 400 parcels. Approximate the probability that the total exceeds 830 kg.

A fair die is rolled 1,000 times. Use the CLT to find an interval that contains the sum with probability about 95%.

Simulate the CLT for Pareto data with $\alpha = 1.5$ (finite mean, infinite variance) and describe what goes wrong.

Interview: Revenue per user in an A/B test is extremely skewed (most users spend nothing, a few spend thousands). The test has 50,000 users per arm. Can you trust a t-test on the mean? Hint: the CLT applies to the mean, but the effective sample size is driven by the few large spenders; check by bootstrapping (Section 8.4), consider trimming or capping, or transform the problem (conversion rate times spend per converter).

Delta Method

Why it matters. The CLT gives the approximate distribution of a sample mean. But we often report a function of means: a log odds ratio, a coefficient of variation, a ratio of two averages (revenue per user), a relative lift in an A/B test, a transformed rate. The delta method extends the CLT to smooth functions of asymptotically normal estimators, giving their standard errors and confidence intervals with one derivative.

The idea: linearize. If $\hat\theta_n$ is close to $\theta$, a first-order Taylor expansion gives

$$g(\hat\theta_n) \approx g(\theta) + g'(\theta)(\hat\theta_n - \theta).$$

A smooth function is approximately linear over the small range where $\hat\theta_n$ varies, and a linear function of an approximately normal variable is approximately normal, with the SD multiplied by $\lvert g'(\theta) \rvert$.

  g(theta)
     ^                     .-'          g is nearly a straight line over the narrow range
     |                 .-'              where theta-hat varies, with slope g'(theta):
     |            .--'
     |        .-' |                     SD of g(theta-hat)  ~  |g'(theta)| x SD of theta-hat
     |     .-'    |
     |  .-'       |
     +---------[--|--]---------> theta
                 theta-hat varies here

Statement (univariate). If $\sqrt n(\hat\theta_n - \theta) \xrightarrow{d} N(0, \sigma^2)$ and $g$ is differentiable at $\theta$ with $g'(\theta) \ne 0$, then

$$\sqrt n\big(g(\hat\theta_n) - g(\theta)\big) \xrightarrow{d} N\big(0,; g'(\theta)^2 \sigma^2\big), \qquad \operatorname{SE}\big(g(\hat\theta)\big) \approx \lvert g'(\hat\theta) \rvert, \operatorname{SE}(\hat\theta).$$

Multivariate version. If $\sqrt n(\hat{\boldsymbol\theta}_n - \boldsymbol\theta) \xrightarrow{d} N(\mathbf{0}, \Sigma)$ and $g: \mathbb{R}^d \to \mathbb{R}$ has gradient $\nabla g$, then $\sqrt n(g(\hat{\boldsymbol\theta}_n) - g(\boldsymbol\theta)) \xrightarrow{d} N(0, \nabla g^\top \Sigma, \nabla g)$. For a ratio $\bar Y/\bar X$, this gives the standard error of metrics like revenue per session, where both numerator and denominator are random.

Worked examples.

Estimator	$g$	Approximate variance
$\log \hat p$ (log of a proportion)	$\log p$	$\dfrac{1 - p}{np}$
log odds $\log\frac{\hat p}{1 - \hat p}$	$g'(p) = \frac{1}{p(1 - p)}$	$\dfrac{1}{np(1 - p)}$
$\sqrt{\bar X}$, Poisson data	$g'(\lambda) = \frac{1}{2\sqrt\lambda}$	$\dfrac{1}{4n}$: does not depend on $\lambda$
$1/\bar X$	$g'(\mu) = -1/\mu^2$	$\dfrac{\sigma^2}{n\mu^4}$
ratio $\bar Y/\bar X$	gradient $(1/\mu_X, -\mu_Y/\mu_X^2)$	$\dfrac{1}{n\mu_X^2}\Big(\sigma_Y^2 - 2R\sigma_{XY} + R^2\sigma_X^2\Big)$, $R = \mu_Y/\mu_X$

Variance-stabilizing transformations. The Poisson row shows a useful trick: when the variance of an estimator depends on its mean, choose $g$ with $g'(\theta) \propto 1/\text{SD}(\theta)$ so that $g(\hat\theta)$ has approximately constant variance. Square roots stabilize Poisson counts, $\arcsin\sqrt{p}$ stabilizes proportions, and logs stabilize data whose SD is proportional to the mean.

Limitations. The approximation is first order: it can be poor for small samples, for strongly curved $g$, or near points where $g'(\theta) = 0$ (then a second-order delta method gives a scaled chi-square limit). It also requires $\hat\theta$ to be approximately normal in the first place. The bootstrap (Section 8.4) is the main alternative when these conditions are doubtful.

import numpy as np

rng = np.random.default_rng(87)
reps = 200_000

# Log odds of a proportion: delta-method SE versus simulation
n, p = 400, 0.2
phat = rng.binomial(n, p, reps) / n
logodds = np.log(phat / (1 - phat))
print(f"log odds: simulated SD {logodds.std():.4f}, delta method {np.sqrt(1 / (n * p * (1 - p))):.4f}")

# Square root of a Poisson mean: variance about 1/(4n) whatever lambda is
n = 50
for lam in [1, 10, 100]:
    xbar = rng.poisson(lam, (reps, n)).mean(axis=1)
    print(f"lambda = {lam:3d}: Var(xbar) = {xbar.var():.4f}, Var(sqrt(xbar)) = {np.sqrt(xbar).var():.5f} (1/(4n) = {1 / (4 * n):.5f})")

# Ratio metric: revenue per session, where users have a random number of sessions
users = 2_000
sessions = rng.poisson(3, (reps // 20, users)) + 1     # X: sessions per user (at least 1)
revenue = rng.gamma(2, 5, sessions.shape) * sessions   # Y: revenue grows with sessions
ratio = revenue.mean(axis=1) / sessions.mean(axis=1)
x, y = sessions[0], revenue[0]                         # one observed experiment
R = y.mean() / x.mean()
var_delta = (y.var() - 2 * R * np.cov(x, y)[0, 1] + R**2 * x.var()) / (users * x.mean() ** 2)
naive = (y / x).std() / np.sqrt(users)                # wrong: treats per-user ratios as the metric
print(f"revenue per session: simulated SE {ratio.std():.4f}, delta method SE {np.sqrt(var_delta):.4f}, "
      f"naive SE of per-user ratios {naive:.4f}")

log odds: simulated SD 0.1259, delta method 0.1250
lambda =   1: Var(xbar) = 0.0199, Var(sqrt(xbar)) = 0.00501 (1/(4n) = 0.00500)
lambda =  10: Var(xbar) = 0.2002, Var(sqrt(xbar)) = 0.00501 (1/(4n) = 0.00500)
lambda = 100: Var(xbar) = 2.0019, Var(sqrt(xbar)) = 0.00500 (1/(4n) = 0.00500)
revenue per session: simulated SE 0.1729, delta method SE 0.1641, naive SE of per-user ratios 0.1515

The delta-method standard errors are close to the simulated ones, even in the ratio case, where the delta SE is computed from a single simulated experiment. In the ratio example, the naive SE answers a different question (the mean of per-user ratios), and it underestimates the uncertainty of the revenue-per-session metric; analysis platforms for A/B tests use the delta method for exactly this kind of metric.

Key takeaways

If $\hat\theta$ is approximately $N(\theta, \sigma^2/n)$, then $g(\hat\theta)$ is approximately $N(g(\theta), g'(\theta)^2\sigma^2/n)$: multiply the SE by $\lvert g' \rvert$.
The multivariate version uses the gradient: $\operatorname{Var} \approx \nabla g^\top \Sigma \nabla g / n$; it gives SEs for ratios, lifts and other functions of several means.
Choosing $g$ with $g' \propto 1/\text{SD}$ gives variance-stabilizing transformations such as $\sqrt{\cdot}$ for Poisson counts.
It is a first-order approximation: check small samples, strong curvature or $g'(\theta) = 0$, or use the bootstrap.

🧪 Practice

A sample of $n = 500$ users has a conversion rate of $\hat p = 0.04$. Use the delta method to find an approximate 95% confidence interval for $\log p$, then transform it back to a confidence interval for $p$.

$\bar X$ is the mean of $n$ i.i.d. Exponential observations with mean $\mu$. Find the approximate variance of the rate estimator $1/\bar X$.

Show that $g(p) = \arcsin\sqrt p$ stabilizes the variance of a sample proportion, and find the constant asymptotic variance.

Interview: The relative lift of an A/B test is reported as $(\bar Y_B - \bar Y_A)/\bar Y_A$. How would you compute its standard error, and when could the result be unreliable? Hint: treat it as $g(\bar Y_A, \bar Y_B) = \bar Y_B/\bar Y_A - 1$ with independent arms and apply the multivariate delta method; trouble arises when $\bar Y_A$ is small or noisy, or the metric is heavy-tailed.

← Probability And Statistics