Lecture 4: Functions of Random Variables, Linearity, Variance

Jamie Haddock

A linear function of a random variable

Let C be a random variable — today’s high temperature in Celsius, say, with some PMF p_C. In Fahrenheit, F = g(C) = \frac{9}{5}C + 32. F is a new random variable on the same sample space as C, defined by F(\omega) = g(C(\omega)).

The PMF of a transformed RV

If g is one-to-one (invertible) on the range of X, the PMF of Y = g(X) is easy: p_Y(y) = P(Y=y) = P(g(X)=y) = P\big(X = g^{-1}(y)\big) = p_X\big(g^{-1}(y)\big).

For Fahrenheit, g^{-1}(f) = \tfrac59(f-32), so p_F(f) = p_C\!\left(\tfrac59(f-32)\right) — the probabilities don’t change, only the labels on the axis do.

LOTUS: the Law Of The Unconscious Statistician

What if we only want E[Y] for Y = g(X), and g isn’t one-to-one — so finding p_Y takes real work (recall grouping outcomes by value, as we did for the dice-game payout in Lecture 3)?

Theorem: LOTUS

For a discrete RV X and any function g, E[g(X)] = \sum_x g(x)\, p_X(x).

The name is a joke: it feels too good to be true, but it is!

Distributions

Definition: Distribution

The distribution of a random variable X describes how the total probability (which always sums to 1) is spread out across the values X can take. Formally, it’s specified by giving P(X \in B) for every relevant subset B of values that X can take (or more generally, every subset B \subset \mathbb{R}).

For a discrete RV, the distribution is fully captured by the PMF p_X(x) = P(X=x) — knowing p_X at every possible value tells you P(X\in B) for any B, just by summing. (For RVs with uncountably many values, a PMF alone won’t do — we’ll need new tools, coming in a later lecture.)

LOTUS example 1: the dice-game payout, revisited

Recall X = die roll (uniform on \{1,\dots,6\}), and Y=g(X) pays $2 on even, -\$1 on odd. In Lecture 3 we found p_Y first, then computed E[Y]. LOTUS skips straight to the answer: E[Y] = \sum_{x=1}^{6} g(x)\cdot \frac16 = - \frac16 + \frac26 - \frac16 + \frac26 - \frac16 + \frac26 = \frac12.

LOTUS example 2: randomized Kaczmarz, revisited

Recall randomized Kaczmarz (RK), solving Ax = b: at each step we sample a row index I \in \{1,\dots,m\} at random and update x_{k+1} = x_k + \frac{b_I - a_I^\top x_k}{\lVert a_I \rVert^2}\, a_I .

A popular distribution over the rows is the row-norm PMF p_I(i) = \|a_i\|^2/\|A\|_F^2. We’ll see why now!

Fix the current iterate x_k and let v = x_k - x^\star. LOTUS lets us compute the expected value of any function of I straight from p_I — for instance E\!\left[\frac{(a_I^\top v)^2}{\|a_I\|^2}\right] = \sum_{i=1}^m \frac{(a_i^\top v)^2}{\|a_i\|^2}\, p_I(i).

Plugging in p_I(i) = \|a_i\|^2/\|A\|_F^2, the \|a_i\|^2 cancels: \sum_{i=1}^m \frac{(a_i^\top v)^2}{\|a_i\|^2}\cdot\frac{\|a_i\|^2}{\|A\|_F^2} = \frac{1}{\|A\|_F^2}\sum_{i=1}^m (a_i^\top v)^2 = \frac{\|Av\|^2}{\|A\|_F^2}.

If I were sampled uniformly instead, the \|a_i\|^2 would not cancel, and the final form wouldn’t be so nice and simple!

Linearity of expectation

Theorem: Linearity of expectation

For random variables X,Y on the same probability space — possibly dependent — and constants a,b,c: E[aX+bY+c] = aE[X]+bE[Y]+c.

No independence required!

This holds no matter how X and Y are related. It’s what makes linearity of expectation so powerful. (Variance will not be so forgiving — stay tuned.)

The mean of a Binomial, the slick way

Recall from Lecture 3 we derived E[X]=np for X\sim\text{Binomial}(n,p) via a combinatorial identity. Linearity gives a one-line proof instead.

Definition: Indicator random variable

For an event A, the indicator \mathbf 1_A equals 1 if A occurs, 0 otherwise. Note E[\mathbf 1_A] = 1\cdot P(A) + 0\cdot P(A^c) = P(A).

Definition: i.i.d.

Random variables X_1,\dots,X_n are independent and identically distributed (i.i.d.) if they are mutually independent and all share the same distribution.

Write binomial random variable X = \sum_{i=1}^n X_i, where X_i = \mathbf 1_{\{\text{trial } i \text{ succeeds}\}} — i.i.d. Bernoulli(p) indicators. By linearity (independence isn’t even needed here!): E[X] = \sum_{i=1}^n E[X_i] = \sum_{i=1}^n p = np.

Linearity also applies to functions of X and Y

Combining LOTUS with linearity, for any functions g,h: E[g(X) + h(Y)] = E[g(X)] + E[h(Y)].

We’ll use exactly this form below.

The coupon collector problem

There are n distinct coupon types; each box of cereal contains one coupon, uniformly random among the n types, independent of previous boxes. Let T = number of boxes you need to buy in order to collect all n types of coupons. What’s E[T]?

Write T = T_1+T_2+\cdots+T_n, where T_i = number of additional boxes needed to find a new type, right after collecting the (i-1)-th distinct type. With i-1 types already in hand, each new box is a new type with probability p_i = \frac{n-i+1}{n}, so T_i\sim\text{Geometric}(p_i) and E[T_i] = \frac{n}{n-i+1}.

By linearity, E[T] = \sum_{i=1}^n E[T_i] = \sum_{i=1}^n \frac{n}{n-i+1} = n\sum_{k=1}^n \frac1k = nH_n \approx n\ln n.

Code

import numpy as np

def collect_coupons(n, rng):
    seen = set()
    draws = 0
    while len(seen) < n:
        seen.add(rng.integers(n))
        draws += 1
    return draws

rng = np.random.default_rng(151)
n = 20
trials = 3000
results = [collect_coupons(n, rng) for _ in range(trials)]

H_n = sum(1 / k for k in range(1, n + 1))
print(f"n = {n}")
print(f"simulated average T over {trials} trials: {np.mean(results):.2f}")
print(f"formula n*H_n = {n*H_n:.2f}")

n = 20
simulated average T over 3000 trials: 72.50
formula n*H_n = 71.95

Linearity + LOTUS: average one-step progress of RK

Recall randomized Kaczmarz’s update, and note that the Pythagorean theorem (the update is an orthogonal projection) provides, \|x_{k+1}-x^\star\|^2 = \|v\|^2 - \frac{(a_I^\top v)^2}{\|a_I\|^2}. where v = x_k-x^\star.

\|v\|^2 is a constant given x_k — only I is random. By linearity, E\big[\|x_{k+1}-x^\star\|^2\big] = \|v\|^2 - E\!\left[\frac{(a_I^\top v)^2}{\|a_I\|^2}\right] = \|v\|^2 - \frac{\|Av\|^2}{\|A\|_F^2}, reusing the LOTUS computation from a few slides ago!

Since \|Av\| \ge \sigma_{\min}(A)\|v\|, E\big[\|x_{k+1}-x^\star\|^2\big] \;\le\; \left(1 - \frac{\sigma_{\min}(A)^2}{\|A\|_F^2}\right)\|x_k-x^\star\|^2. Every step shrinks the expected squared error by a fixed factor — this a very nice (and very well-cited) result from Thomas Strohmer and Roman Vershynin in 2009!

Watching the bound play out over many steps

The one-step bound says the expected squared error shrinks by a fixed factor every step. Let’s actually run RK and watch \|x_k-x^\star\|^2 fall, averaged over many independent random runs.

Code

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(151)
m, n = 20, 8
A = rng.standard_normal((m, n))
x_star = rng.standard_normal(n)
b = A @ x_star
row_norm_sq = (A ** 2).sum(axis=1)
p_row_norm = row_norm_sq / row_norm_sq.sum()
sigma_min = np.linalg.svd(A, compute_uv=False)[-1]
frob_sq = row_norm_sq.sum()

num_steps = 60
num_runs = 300
errors_sq = np.zeros((num_runs, num_steps + 1))

for r in range(num_runs):
    x = np.zeros(n)
    errors_sq[r, 0] = np.linalg.norm(x - x_star) ** 2
    idx_sequence = rng.choice(m, size=num_steps, p=p_row_norm)
    for t, i in enumerate(idx_sequence, start=1):
        a_i, b_i = A[i], b[i]
        x = x + (b_i - a_i @ x) / row_norm_sq[i] * a_i
        errors_sq[r, t] = np.linalg.norm(x - x_star) ** 2

mean_error_sq = errors_sq.mean(axis=0)
bound_curve = errors_sq[:, 0].mean() * (1 - sigma_min**2 / frob_sq) ** np.arange(num_steps + 1)

fig, ax = plt.subplots(figsize=(6, 3.5))
ax.semilogy(mean_error_sq, label="average of $\\|x_k-x^\\star\\|^2$ over runs")
ax.semilogy(bound_curve, "--", label="one-step bound, applied repeatedly")
ax.set_xlabel("step $k$")
ax.set_ylabel("squared error (log scale)")
ax.legend()
plt.tight_layout()
plt.show()

Moments

Definition: Moments

The k-th moment of a random variable X is E[X^k]. The mean E[X] is the 1st moment. The 2nd moment, E[X^2], is about to become useful.

Why do we need more than the mean?

Suppose you only know E[X] where X = \|x_k-x^\star\|^2 is small (e.g., from the bound we just derived). What does that tell you about the error on any one particular run of RK iterations?

Not much, on its own! The mean is an average over all randomness; a single run could still land far from it. To say anything about how spread out outcomes are around the mean, we need a new tool: variance.

Aside: the probabilistic method

Expectation alone can still be a powerful existence tool: if E[X]\le c, then some outcome \omega must satisfy X(\omega)\le c (otherwise every outcome would exceed c, forcing E[X]>c). This “probabilistic method” is often used to prove existence results in combinatorics — but it says nothing about how many outcomes are near c, or how far the rest can stray. That’s what variance is for.

Variance and standard deviation

Definition: Variance

\mathrm{Var}(X) = E\big[(X-E[X])^2\big]. The standard deviation is \sigma_X = \sqrt{\mathrm{Var}(X)} — same units as X itself (variance is in squared units), which is often why it’s the more interpretable number.

Variance of a linear function

For constants a,b: \mathrm{Var}(aX+b) = a^2\,\mathrm{Var}(X).

Proof: E[aX+b] = aE[X]+b, so \mathrm{Var}(aX+b) = E\big[(aX+b-aE[X]-b)^2\big] = E\big[a^2(X-E[X])^2\big] = a^2\,\mathrm{Var}(X). The additive shift b vanishes entirely — variance sees only spread, never location.

Variance shortcut formula

Theorem

\mathrm{Var}(X) = E[X^2] - \big(E[X]\big)^2.

Proof: expand (X-E[X])^2 = X^2-2XE[X]+(E[X])^2 and take E[\cdot], using linearity (note E[X] is a constant): E\big[(X-E[X])^2\big] = E[X^2]-2E[X]E[X]+(E[X])^2 = E[X^2]-(E[X])^2. This is almost always the easiest way to compute a variance.

Compute \mathrm{Var}(X) for each named distribution

Bernoulli(p): X^2=X on \{0,1\}, so E[X^2]=E[X]=p, giving \mathrm{Var}(X)=p-p^2=p(1-p).

Rademacher: X^2=1 always, so E[X^2]=1, E[X]=0, giving \mathrm{Var}(X)=1.

In-class activity

Using \mathrm{Var}(X)=E[X^2]-(E[X])^2, compute \mathrm{Var}(X) for Geometric and Poisson.

Geometric(p): Using the derivative series trick twice gives E[X^2]=\dfrac{2-p}{p^2}, so \mathrm{Var}(X) = \frac{2-p}{p^2}-\frac1{p^2} = \frac{1-p}{p^2}.

Poisson(\lambda): E[X^2]=\lambda^2+\lambda, which gives \mathrm{Var}(X) = \lambda^2+\lambda-\lambda^2 = \lambda. (Poisson is the rare distribution with mean = variance!)

Variance of a sum — independence is required this time

Contrast with linearity of expectation

E[X+Y]=E[X]+E[Y] always. But in general \mathrm{Var}(X+Y)\ne\mathrm{Var}(X)+\mathrm{Var}(Y) — cross terms from how X and Y move together get in the way.

Theorem

If X_1,\dots,X_n are independent, then \mathrm{Var}\!\left(\sum_{i=1}^n X_i\right) = \sum_{i=1}^n \mathrm{Var}(X_i).

Payoff: writing Binomial(n,p) as a sum of n i.i.d. Bernoulli(p) indicators: \mathrm{Var}(X) = \sum_{i=1}^n \mathrm{Var}(X_i) = np(1-p).

Central moments

Definition: Central moments

The k-th central moment of X is E[(X-E[X])^k] — an ordinary moment, but centered at the mean first.

The 1st central moment is always 0 (by definition of the mean).
The 2nd central moment is exactly the variance — nothing new, just a name.
Higher central moments capture other shape features: the 3rd relates to skewness (asymmetry), the 4th to kurtosis (tail heaviness). We won’t need these in this course, but variance is really just the first member of this whole family.