Roll two dice. Chapter 29 would ask "what’s the probability of a total of 8?". But often you care about the total itself as a quantity: what total should you expect on average? How variable is it? You want to treat "the total" as a single object you can compute with, the way a program treats a variable.

A random variable is a rule that turns each outcome of an experiment into a number. For two dice, the rule "add the faces" turns the outcome \((3, 5)\) into 8 and \((6, 6)\) into 12. Before the dice are rolled you don’t know the value; afterwards you do. Different rules on the same experiment give different random variables: "the larger face", "the number of sixes", "1 if doubles else 0".

For a discrete random variable, the distribution is the list of possible values with their probabilities. Let \(S\) be the total of two dice. Counting cells in the \(6 \times 6\) grid from Chapter 29 gives:

A game costs £1 to play. You roll a die; on a 6 you win £5, otherwise nothing. Should you play? Each game is a gamble, but over a thousand games, what happens to your money?

Over 6,000 games you’d expect about 1,000 sixes, winning \(1000 \times 5 = 5000\) pounds, while paying 6,000. You lose about £1,000, or about 17p per game. That per-game average is the expected value. It’s what the average of many repetitions settles down to (the frequency interpretation from Chapter 29, applied to quantities).

Computed directly: you win £5 with probability \(\frac16\) and £0 with probability \(\frac56\), so the average winnings per game are \(5 \times \frac16 + 0 \times \frac56 = \frac56 \approx 0.83\); subtract the £1 stake and the expected profit is \(-\frac16\) pound.

Let \(X\) be the face of a fair die.

You’ll never roll 3.5. The expected value is a balance point, not a "most likely value". Picture the PMF as weights on a seesaw: \(E[X]\) is where it balances. For the two-dice total, the symmetric PMF balances at 7, so \(E[S] = 7\).

The same idea works for any function of \(X\): \(E[g(X)] = \sum_x g(x)\,P(X = x)\). For example \(E[X^2] = \frac{1 + 4 + 9 + 16 + 25 + 36}{6} = \frac{91}{6}\). Note that \(E[X^2] = 15.17 \ne (E[X])^2 = 12.25\). In general you can’t push \(E\) through a nonlinear function.

The most useful fact about expectation is almost embarrassingly simple:

for any random variables \(X\) and \(Y\), whether or not they are independent.

Why it works. Think of the sample space as a finite list of outcomes \(\omega\), each with probability \(P(\omega)\). Then \(E[X] = \sum_\omega X(\omega)P(\omega)\) (group the outcomes by value and you get the earlier formula back). Now

Nothing about independence was used. A sum is a sum.

So the total of two dice has \(E[S] = E[X_1] + E[X_2] = 3.5 + 3.5 = 7\), without the PMF table. Ten dice: 35. A hundred: 350.

Linearity becomes a superpower when combined with indicators.

To find the expected number of things that happen, write the count as a sum of indicators, one per thing, and add their probabilities. You never need the distribution of the count.

\(n\) people check their hats at a restaurant. The attendant loses the tickets and hands the hats back in a uniformly random order. On average, how many people get their own hat back?

The events "person \(i\) gets their own hat" are heavily dependent (if the first \(n-1\) people all get their own hats, so does the last). Computing the distribution of the number of matches is a genuinely hard counting problem. Linearity ignores all that.

Exactly one match on average, whether there are 3 people or 3 million. (A harder inclusion–exclusion argument shows that the probability of no matches tends to \(1/e \approx 0.368\).)

Two servers both average 100 ms per request. One always takes 95–105 ms. The other takes 10 ms most of the time and occasionally 2 seconds. The means agree; the experiences don’t. You need a number for spread.

Measure how far \(X\) typically lands from its mean \(\mu\). The deviation \(X - \mu\) averages to zero (positives and negatives cancel: \(E[X - \mu] = \mu - \mu = 0\)), so square it first, making every deviation count positively, then average.

Expanding the square and using linearity gives a formula that’s usually easier:

"Mean of the square minus square of the mean."

From before, \(E[X] = 3.5\) and \(E[X^2] = 91/6\).

A die typically lands about 1.7 away from 3.5, which matches the eye: the values 1 to 6 are between 0.5 and 2.5 from the centre.

Adding a constant \(b\) shifts everything and doesn’t change spread; scaling by \(a\) scales deviations by \(a\), so squared deviations by \(a^2\). The sum rule needs independence (in general there’s an extra covariance term, Section 30.8). Note the consequence: for the average of \(n\) independent copies of \(X\),

Averaging \(n\) measurements shrinks the spread by a factor of \(\sqrt n\). To halve the noise, take four times as many measurements. Hold on to this: it’s the engine of the law of large numbers, the central limit theorem and nearly all of Chapter 31.

Story: one trial that succeeds with probability \(p\) (a coin flip, a click or not, a packet lost or not). \(X = 1\) on success, 0 on failure.

The variance: \(E[X^2] = p\) (since \(X^2 = X\) when \(X\) is 0 or 1), so \(\operatorname{Var} = p - p^2 = p(1 - p)\). It’s largest at \(p = 1/2\), the most uncertain coin.

Story: \(n\) independent Bernoulli trials, each with success probability \(p\). \(X\) = number of successes. Write \(X \sim \text{Bin}(n, p)\), read "\(X\) has the binomial distribution with parameters \(n\) and \(p\)".

Deriving the PMF by counting. Take \(n = 5\), and ask for exactly \(k = 2\) successes. One particular sequence, say SSFFF, has probability \(p \cdot p \cdot (1-p)(1-p)(1-p) = p^2(1-p)^3\) by independence. Every sequence with exactly 2 S’s and 3 F’s has the same probability, just in a different order. How many such sequences? Choose which 2 of the 5 positions are S: \(\binom52 = 10\) (Chapter 24). The sequences are disjoint, so add: \(P(X = 2) = 10\,p^2(1-p)^3\). In general:

The binomial theorem confirms these sum to \((p + (1 - p))^n = 1^n = 1\); that’s where the distribution gets its name.

Mean and variance, without algebra. \(X = X_1 + \cdots + X_n\), a sum of independent Bernoullis. By linearity \(E[X] = np\); by independence \(\operatorname{Var}(X) = np(1 - p)\).

Worked example. Ten fair coin flips. Probability of exactly 7 heads? Try it: Binomial(10, 0.5) at 7.

Probability of 8 or more heads: \(\frac{\binom{10}{8} + \binom{10}{9} + \binom{10}{10}}{1024} = \frac{45 + 10 + 1}{1024} = \frac{56}{1024} \approx 0.055\).

The brute-force mean and variance match \(np = 5\) and \(np(1-p) = 2.5\).

Story: repeat independent trials with success probability \(p\) until the first success. \(X\) = number of trials needed (including the success).

The first success is on trial \(k\) exactly when the first \(k - 1\) trials fail and trial \(k\) succeeds:

The mean, by a conditioning trick. Let \(\mu = E[X]\). Look at the first trial. With probability \(p\) it succeeds, and you’re done after 1 trial. With probability \(1 - p\) it fails, you’ve used 1 trial, and, because trials are independent, you’re back where you started, facing another \(\mu\) trials on average. So

If a request fails with probability 0.25, so succeeds with probability \(p = 0.75\), a retry loop needs \(1/0.75 \approx 1.33\) attempts on average. Rolling a die until you get a six takes 6 rolls on average. (The variance is \((1 - p)/p^2\).)

The geometric distribution is memoryless: having already failed 10 times doesn’t make success any more likely on trial 11. That’s the gambler’s fallacy from Chapter 29 in distribution form.

The problem. A web server receives on average 3 requests per second, from a huge population of users each of whom is very unlikely to send a request in any particular second. What’s the probability of 0 requests in a given second? Of 6 or more?

Story: count events that happen independently, at a constant average rate, where each tiny time slice has only a tiny chance of an event. \(X \sim \text{Poisson}(\lambda)\), where \(\lambda\) (lambda) is the average count per interval:

Where it comes from. Chop the second into \(n\) tiny slices, each containing a request with probability \(p = \lambda/n\). The count is binomial, \(\text{Bin}(n, \lambda/n)\). As \(n \to \infty\) with \(np = \lambda\) fixed,

by the compound-interest limit from Chapter 12, and similar limits give the other terms. So the Poisson is the binomial for "many trials, each rare". Numerically: \((1 - 3/1000)^{1000} \approx 0.0496\), already close to \(e^{-3} \approx 0.0498\).

Worked example. With \(\lambda = 3\): Try it: Poisson with λ = 3.

About 1 second in 20 sees no traffic, and about 1 in 12 sees double the average or more. If the server can handle 5 requests per second, it will be overloaded 8% of the time, even though the average load is only 60% of capacity. Randomness creates queues.

How long until the next request arrives? It could be 0.3 s, 0.31 s, 0.3141 s: any real number. Now ask "what’s the probability it’s exactly 0.3 seconds, to infinitely many decimal places?". There are infinitely many possible values; if each had positive probability, the total would be infinite. So each individual value must have probability 0. Yet "between 0.2 and 0.4 seconds" clearly has positive probability. Lists of probabilities can’t work; we need a different description.

Think of a histogram of a million measured waiting times, with very narrow bins, scaled so the total area of the bars is 1. The probability of landing in any range is the area of the bars over that range. As the bins get narrower, the bar tops trace out a smooth curve: the probability density function. Probability is area under it.

One description works for discrete and continuous variables alike.

Story: \(X\) is equally likely to land anywhere in \(\lbrack a, b]\). The density is flat: \(f(x) = \frac{1}{b - a}\) on the interval (so the rectangle has area 1), and 0 outside. Probabilities are proportional to length: \(P(c \le X \le d) = \frac{d - c}{b - a}\). The mean is the midpoint \(\frac{a + b}{2}\), and the variance is \(\frac{(b - a)^2}{12}\).

The standard uniform, on \(\lbrack 0, 1]\), is what random.random() approximates, and as you’ll see in [sec-dist-sampling], every other distribution can be built from it.

Story: the waiting time until the next event in a Poisson process with rate \(\lambda\) events per unit time. The time until the next request, the next radioactive decay, the next failure of a component that doesn’t age.

Deriving it from Poisson. The wait \(T\) exceeds \(t\) exactly when there are zero events in \(\lbrack 0, t]\). In time \(t\) the count is Poisson with mean \(\lambda t\), so

The density is the slope of \(F\): \(f(t) = \lambda e^{-\lambda t}\) for \(t \ge 0\). The mean is \(1/\lambda\) (at 3 requests per second, the average gap is a third of a second), and the standard deviation is also \(1/\lambda\).

Worked example. Buses arrive as a Poisson process at rate \(\lambda = 0.5\) per minute. What’s the probability your wait is between 1 and 3 minutes? Try it: Exponential waiting time.

Memorylessness. \(P(T > s + t \mid T > s) = \frac{e^{-\lambda(s + t)}}{e^{-\lambda s}} = e^{-\lambda t} = P(T > t)\). If you’ve already waited 10 minutes, the remaining wait has the same distribution as if you’d just arrived. That’s right for radioactive atoms and for well-behaved request streams, and wrong for things that wear out (hard disks, people), which need other models.

Measure the heights of 10,000 adults, the errors of 10,000 repeated measurements of a length, or the average response time across 10,000 batches of requests. Plot histograms, and the same symmetric bell shape keeps appearing. Why one shape, and how do you compute with it?

The "why" is the central limit theorem (Section 30.7): any quantity that is the sum of many small, independent contributions comes out approximately bell-shaped. First, the "how".

You don’t need to work with the formula directly. Two facts do almost all the work.

Standardizing. If \(X \sim N(\mu, \sigma^2)\), then

\(Z\) (the z-score) counts how many standard deviations \(X\) is from its mean. So every normal question becomes a question about one curve.

The 68–95–99.7 rule. For any normal distribution, about 68% of the probability lies within \(1\sigma\) of the mean, 95% within \(2\sigma\), and 99.7% within \(3\sigma\).

Worked example. IQ scores are scaled to be approximately \(N(100, 15^2)\). What fraction score above 130? Try it: Normal tail above 130.

Quick check with the rule: 95% lie within \(\pm 2\sigma\), so 5% lie outside, half of it above: 2.5%. ✓

(The function erf, the error function, is in Python’s math module; \(\Phi(z) = \tfrac12\bigl(1 + \operatorname{erf}(z/\sqrt2)\bigr)\). The last line computes \(P(150 \le X \le 250)\) for \(X \sim N(200, 40^2)\), i.e. \(P(-1.25 \le Z \le 1.25)\).)

A single die roll is unpredictable. The average of 10,000 rolls is extremely predictable: it will be very close to 3.5. A single user’s session length is unpredictable; the mean session length across a million users barely moves day to day. Why does averaging tame randomness, and how fast?

Let \(X_1, X_2, \ldots, X_n\) be independent random variables, all with the same distribution (independent and identically distributed, or i.i.d.), with mean \(\mu\) and standard deviation \(\sigma\). Their sample mean is

(\(\bar X\) is read "X-bar".) We computed in Section 30.3 that \(E[\bar X_n] = \mu\) and \(\operatorname{SD}(\bar X_n) = \sigma/\sqrt n\). As \(n\) grows, the spread of the average around \(\mu\) shrinks to zero. That is the law of large numbers: the sample mean converges to the true mean.

To make "converges" precise, Chebyshev’s inequality says that any random variable is unlikely to be many standard deviations from its mean: \(P(|Y - E[Y]| \ge k\sigma_Y) \le 1/k^2\). Applied to \(\bar X_n\), with a fixed tolerance \(\varepsilon\) (epsilon, a small positive number), Try it: Run the law of large numbers.

This is why the frequentist interpretation of probability (Chapter 29) is coherent: the proportion of heads is the sample mean of Bernoulli indicators, so it converges to \(p\). Try it: Proportion of heads converging.

The law of large numbers says the sample mean concentrates around \(\mu\). The central limit theorem (CLT) says what shape its fluctuations have: approximately normal, no matter what distribution the individual \(X_i\) have (provided their variance is finite).

Watch it happen with an extremely un-normal starting distribution, the exponential (all mass on the positive side, a long right tail):

Why should sums become normal? A rough intuition: each extra term adds a small independent nudge; the fine details of each nudge’s distribution wash out, and only the mean and variance survive. A sum of many small independent nudges is the situation the normal distribution describes. (The full proof uses more machinery than this book develops.) Try it: Build the bell curve yourself; Exponential means, n = 1.

Read the output carefully. The standard deviation of the mean matches \(\sigma/\sqrt n\) exactly at every \(n\): that part is not an approximation. The shape approaches normal more slowly. For a normal distribution, 2.5% would fall in each tail beyond \(\pm 1.96\) standard errors; with skewed data, even at \(n = 30\) the tails are lopsided (1.4% versus 3.3%).

Try it: Watch the CLT fail (Cauchy); Cauchy running means never settle.

Do slower pages have higher bounce rates? Does CPU load rise with request count? These questions are about two random variables at once and whether they move together.

The joint distribution of discrete \(X\) and \(Y\) gives \(P(X = x, Y = y)\) (the comma means "and") for every pair of values. It’s a two-way table, exactly like the deploy table in Chapter 29.

Example. Roll a die. Let \(X = 1\) if the roll is even (else 0), and \(Y = 1\) if the roll is at least 4 (else 0). Rolls: 1 → (0,0), 2 → (1,0), 3 → (0,0), 4 → (1,1), 5 → (0,1), 6 → (1,1).

The row and column totals, written in the margins, are the marginal distributions of \(X\) and \(Y\) on their own. \(X\) and \(Y\) are independent if every cell equals the product of its marginals: \(P(X = x, Y = y) = P(X = x)P(Y = y)\). Here \(P(X = 1, Y = 1) = 2/6 \ne \frac12 \cdot \frac12\), so they’re dependent: knowing the roll is even makes "at least 4" more likely.

Why the product? If both deviations have the same sign, their product is positive; if opposite signs, negative. Averaging the products measures which tendency wins.

For the die example: \(E[XY] = P(X = 1, Y = 1) = \frac26 = \frac13\) (since \(XY\) is 1 only when both are 1), and \(\mu_X = \mu_Y = \frac12\). So \(\operatorname{Cov}(X, Y) = \frac13 - \frac14 = \frac1{12}\): positive, as expected.

Covariance has awkward units (milliseconds times bounces) and depends on scale. Dividing by both standard deviations gives a pure number between \(-1\) and 1:

\(\rho\) (rho) is the correlation coefficient. \(\rho = 1\) means \(Y\) is an exact increasing linear function of \(X\); \(\rho = -1\), an exact decreasing one; \(\rho = 0\), no linear relationship. For the die example, \(\sigma_X = \sigma_Y = \sqrt{\frac12 \cdot \frac12} = \frac12\), so \(\rho = \frac{1/12}{1/4} = \frac13\).

Linear search. Search an array of \(n\) distinct items for one that’s equally likely to be in any position. The number of comparisons is uniform on \(\{1, \ldots, n\}\), so the expected number is \(\frac{1 + 2 + \cdots + n}{n} = \frac{n + 1}{2}\).

Hash table chains. Insert \(n\) keys into \(m\) buckets with a good hash. For a given bucket, the number of keys in it is a sum of \(n\) indicators each with probability \(1/m\), so by linearity its expected length is \(n/m\), the load factor \(\alpha\). An unsuccessful lookup scans a chain of expected length \(\alpha\): constant time if you keep \(\alpha\) bounded by resizing.

Randomized quicksort (promised in Chapter 29). Let the sorted values be \(z_1 < z_2 < \cdots < z_n\). Quicksort only ever compares an element with the current pivot, and any pair is compared at most once. Let \(I_{ij} = 1\) if \(z_i\) and \(z_j\) (\(i < j\)) are ever compared. Total comparisons \(C = \sum_{i<j} I_{ij}\).

Focus on the block \(z_i, z_{i + 1}, \ldots, z_j\) (\(j - i + 1\) values). Until a pivot is chosen from this block, all of them stay together in the same subarray. When the first pivot from the block is chosen, if it’s \(z_i\) or \(z_j\), those two get compared; if it’s anything in between, they’re split into different subarrays and never meet. Every element of the block is equally likely to be that first pivot, so

where \(H_n = 1 + \frac12 + \frac13 + \cdots + \frac1n \approx \ln n + 0.577\) is the \(n\)-th harmonic number. (The exact value is \(2(n + 1)H_n - 4n\).) No input is bad; only the coin flips are random.

The input list(range(n)) is already sorted, the worst case for a first-element pivot (about 500,000 comparisons); with random pivots it costs about 11,000.

Try it: Collect 50 coupons. ==== The coupon collector

The problem. A cereal brand puts one of \(n\) equally likely toys in each box. How many boxes until you have all \(n\)? The same question: how many random requests until every one of \(n\) servers has been hit; how many random tests until every one of \(n\) code paths is covered.

The method. Split the process into phases. Phase \(k\) starts when you have \(k - 1\) distinct toys and ends when you get a new one. During phase \(k\), each box is new with probability \(p_k = \frac{n - k + 1}{n}\), so the phase length is geometric with mean \(\frac{n}{n - k + 1}\). By linearity,

For \(n = 50\): \(50 \times H_{50} \approx 50 \times 4.499 = 225\) boxes. The last toy alone takes 50 boxes on average; the first takes 1.

Try it: Balls into bins, two choices. ==== Balls into bins: load balancing

Throw \(n\) requests (balls) at \(n\) servers (bins) uniformly at random.

The power of two choices: for each ball, pick two random bins and put the ball in the less loaded one. This tiny change cuts the maximum load to about \(\frac{\ln \ln n}{\ln 2}\), an exponential improvement. It’s used in real load balancers (for example, NGINX’s random two method) and in hash tables (cuckoo hashing).

Try it: Inverse-transform sampling, animated. ==== Sampling from a distribution: inverse transform sampling

The problem. Your language gives you random(), uniform on \(\lbrack 0, 1)\). A simulation needs exponential inter-arrival times, or a weighted choice among options. How do you turn uniform randomness into any distribution?

Intuition. Look at the CDF \(F\). It maps values to probabilities in \(\lbrack 0, 1]\). Run it backwards: pick a height \(u\) uniformly on the vertical axis, move across to the curve, and drop down to the \(x\) axis. Where the CDF is steep (high density), a wide band of heights maps to a narrow band of \(x\), so those \(x\) values get picked often. Exactly right.

Why it works. Let \(U\) be uniform on \(\lbrack 0,1]\) and \(X = F^{-1}(U)\) (the inverse function, Chapter 11). Since \(F\) is increasing, \(X \le x\) exactly when \(U \le F(x)\), so \(P(X \le x) = P(U \le F(x)) = F(x)\): the uniform’s CDF is \(P(U \le u) = u\). So \(X\) has CDF \(F\), as required.

For the exponential, solve \(u = 1 - e^{-\lambda x}\) for \(x\):

For a discrete distribution, the CDF is a staircase, and "inverting" it means: walk up the staircase until the cumulative probability first exceeds \(u\). Try it: Walk up the CDF staircase.

Alternatives and choosing. Inverse transform is the default whenever \(F^{-1}\) has a formula. The normal distribution doesn’t, so libraries use special methods (the Box–Muller transform, or the faster ziggurat algorithm). For discrete distributions with many values sampled many times, the linear walk is slow; use binary search on the cumulative array (\(O(\log n)\), this is what Python’s random.choices does) or the alias method (\(O(1)\) per sample after setup). Rejection sampling handles densities you can evaluate but not invert.

Try it: Step through reservoir sampling. ==== Reservoir sampling

The problem. Choose \(k\) items uniformly at random from a stream (log lines, tweets, packets) whose length you don’t know in advance and which is too big to store. You get one pass.

The method (Algorithm R). Keep the first \(k\) items in a "reservoir". For the \(i\)-th item (\(i > k\)), keep it with probability \(k/i\), replacing a uniformly random reservoir slot.

Why it works (the case \(k = 1\)). Keep item \(i\) with probability \(1/i\). Item \(i\) survives to the end of a stream of length \(n\) if it’s chosen at step \(i\) and not replaced at any later step \(i + 1, \ldots, n\):

Every item ends up equally likely, whatever \(n\) turns out to be. The general-\(k\) argument is the same with \(k/i\) and a factor \((1 - \frac{1}{j})\) for "replaced at step \(j\) and that slot was mine", and gives \(k/n\) for every item (an induction proof, Chapter 22).

Each of the 10 items appears in the sample with probability \(3/10\), as promised.