4. Random Variables and Distributions

A random variable turns the outcome of a random experiment into a number, so that probability can be combined with arithmetic, algebra and calculus. This chapter defines discrete and continuous random variables and the functions that describe their distributions (PMF, PDF, CDF and quantile function), summarizes distributions with expectation, variance, moments and generating functions, and shows how a distribution changes when its random variable is transformed. Every named distribution in Chapter 5, and every estimator and test later in the guide, is described in this vocabulary.

<a id="41-discrete-random-variables"></a>

4.1 Discrete Random Variables

Discrete random variables take finitely or countably many values, which lets us describe their distributions with simple lists and sums before moving to the calculus that continuous variables require.

Definition of a Random Variable

Why it matters. Outcomes of experiments are often not numbers: "HTH", "the user clicked the second ad", "the server returned 503". We cannot average them, compare them or plot them. What we usually care about is a number computed from the outcome: the count of heads, the revenue of a session, the latency of a request. A random variable formalizes that computation, and once we have it, all of algebra and calculus become available for reasoning about uncertainty.

Definition. A random variable is a function that assigns a real number to every outcome of a random experiment:

$$X: \Omega \to \mathbb{R}, \qquad \omega \mapsto X(\omega).$$

Section 3.5 added the technical requirement that $X$ be measurable, so that $P(X \in B)$ exists for every Borel set $B$. In practice every function you will write down satisfies it.

The name is misleading on both counts: a random variable is neither random nor a variable. It is a fixed, deterministic function. All the randomness lives in which outcome $\omega$ the experiment produces; $X$ just reads a number off it. A useful analogy is a sensor: nature picks the outcome, and the sensor reports a measurement of it. Different sensors on the same experiment are different random variables.

Notation.

Symbol	Meaning
$X, Y, Z$ (capital)	random variables (functions)
$x, y, z$ (lowercase)	particular real values
${X = x}$	the event ${\omega : X(\omega) = x}$
${X \le x}$, ${X \in B}$	the events ${\omega : X(\omega) \le x}$, ${\omega : X(\omega) \in B}$
$P(X = x)$	the probability of the event ${X = x}$
$X \sim \text{Distribution}$	"$X$ has this distribution"
$X \overset{d}{=} Y$	$X$ and $Y$ have the same distribution

Many random variables, one experiment. Toss a coin three times. The sample space has 8 equally likely outcomes, and we can define as many random variables on it as we like:

   outcome w    X(w) = #heads    Y(w) = 1{first is H}    R(w) = longest run
   ---------    -------------    --------------------    ------------------
     HHH              3                   1                      3
     HHT              2                   1                      2
     HTH              2                   1                      1
     ...             ...                 ...                    ...
     TTT              0                   0                      3

$Y$ is an indicator random variable: it equals 1 if an event happens and 0 otherwise. For an event $A$, the indicator is written $\mathbf{1}_A$. Indicators are the bridge between events and random variables, and Section 4.3 will turn them into one of the most powerful counting tools in probability.

The distribution of a random variable. The distribution (or law) of $X$ is the collection of probabilities $P(X \in B)$ for all sets $B$ of values. It describes which values $X$ takes and how often, and it forgets the underlying sample space. Two random variables can have the same distribution without being the same function: with three tosses, the number of heads $X$ and the number of tails $T = 3 - X$ have identical distributions, yet $X$ and $T$ are never equal on any outcome.

Types of random variables.

Type	Values	Described by	Examples
Discrete	finite or countable set	probability mass function (Section 4.1)	number of defects, die roll, clicks per hour
Continuous	a continuum, $P(X = x) = 0$ for every $x$	probability density function (Section 4.2)	height, latency, temperature
Mixed	point masses plus a continuous part	CDF, or a mass function plus a density	rainfall (exactly 0 on dry days), insurance claims

from collections import Counter
from fractions import Fraction
from itertools import groupby, product

# Sample space: three tosses of a fair coin, 8 equally likely outcomes
omega = ["".join(t) for t in product("HT", repeat=3)]

# Several random variables on the SAME sample space: each is just a function of the outcome
X = lambda w: w.count("H")                                  # number of heads
T = lambda w: w.count("T")                                  # number of tails, T = 3 - X
Y = lambda w: int(w[0] == "H")                              # indicator of "first toss is heads"
R = lambda w: max(len(list(run)) for _, run in groupby(w))  # length of the longest run

print("outcome  X  T  Y  R")
for w in omega:
    print(f"  {w}    {X(w)}  {T(w)}  {Y(w)}  {R(w)}")

def distribution(V):
    """Push the uniform probability on omega forward to the values of V."""
    d = Counter()
    for w in omega:
        d[V(w)] += Fraction(1, len(omega))
    return {v: str(p) for v, p in sorted(d.items())}

for name, V in [("X", X), ("T", T), ("Y", Y), ("R", R)]:
    print(f"distribution of {name}: {distribution(V)}")

# Same distribution does not mean same random variable
print("P(X = T) =", Fraction(sum(X(w) == T(w) for w in omega), len(omega)))

outcome  X  T  Y  R
  HHH    3  0  1  3
  HHT    2  1  1  2
  HTH    2  1  1  1
  HTT    1  2  1  2
  THH    2  1  0  2
  THT    1  2  0  1
  TTH    1  2  0  2
  TTT    0  3  0  3
distribution of X: {0: '1/8', 1: '3/8', 2: '3/8', 3: '1/8'}
distribution of T: {0: '1/8', 1: '3/8', 2: '3/8', 3: '1/8'}
distribution of Y: {0: '1/2', 1: '1/2'}
distribution of R: {1: '1/4', 2: '1/2', 3: '1/4'}
P(X = T) = 0

Key takeaways

• A random variable is a deterministic function from outcomes to real numbers; the randomness comes from the outcome.
• Events about a random variable, such as ${X \le 3}$, are ordinary events in $\Omega$, so all the rules of Chapter 3 apply.
• The distribution of $X$ records how probability is spread over its values and forgets $\Omega$.
• Equal in distribution ($X \overset{d}{=} Y$) is much weaker than equal ($X = Y$).
• Random variables are discrete, continuous or mixed; each type has its own descriptive tools.

🧪 Practice

1. Roll two dice. Define random variables for (a) the larger of the two values, (b) the absolute difference, (c) the indicator that the dice show the same number. List the event ${\text{(b)} = 0}$ and compare it with the event ${\text{(c)} = 1}$.
2. Classify each random variable as discrete, continuous or mixed: the number of emails received today; the time until the next email; the amount of money a visitor spends on a website (most visitors spend nothing); the proportion of heads in 10 tosses.
3. Extend the code to four tosses. Find two different random variables (other than heads and tails) that have the same distribution.
4. Interview: A teammate says "a random variable is just a column in a dataframe". In what sense is this a good mental model, and where does it break down? Hint: think about what one row represents (one outcome $\omega$), what the column does (maps that outcome to a number), and what is missing (the probability measure that says how likely each row is, and the outcomes that were never observed).

Probability Mass Function

Why it matters. For a discrete random variable, the most direct description of its distribution is a list: each possible value, with its probability. That list is the probability mass function. Every probability about $X$ is a sum over it, every discrete model in Chapter 5 is specified by one, and the bar chart of a PMF is often the first thing to look at when reasoning about a count.

Definition. The probability mass function (PMF) of a discrete random variable $X$ is

$$p_X(x) = P(X = x), \qquad x \in \mathbb{R}.$$

It is zero except at the countably many values $X$ can take. A function $p$ is a valid PMF if and only if:

1. $p(x) \ge 0$ for every $x$;
2. $\sum_x p(x) = 1$, where the sum runs over the values with $p(x) > 0$.

Then, for any set $B$ of values,

$$P(X \in B) = \sum_{x \in B} p_X(x).$$

Analogy. Think of one kilogram of sand poured onto a number line, landing only in a few piles. The PMF tells you how much sand is in each pile. Probability "mass" is a literal image: the total is fixed at 1, and a probability is the amount of mass inside a region.

Where PMFs come from.

Source	Example	How $p(x)$ is obtained
Counting equally likely outcomes	sum of two dice	(outcomes with value $x$) / (total outcomes)
A model formula	geometric: tosses until the first head	$p(k) = (1 - q)^{k-1} q$ for $k = 1, 2, \dots$
Data	observed goals per football match	relative frequency of $x$ in the sample (the empirical PMF)

A model formula must be normalized. For example, $p(k) = c / k^2$ for $k = 1, 2, \dots$ is a valid PMF only for $c = 1 / \sum_k k^{-2} = 6 / \pi^2$. Normalizing constants like $c$ appear throughout statistics.

The empirical PMF is the discrete version of the relative frequency table from Section 2.1. It estimates the true PMF and converges to it as the sample grows, but on a small sample it is noisy and assigns probability 0 to every value that has not been observed yet.

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

# Exact PMF of S = sum of two fair dice, by counting the 36 equally likely outcomes
pmf = Counter()
for d1, d2 in product(range(1, 7), repeat=2):
    pmf[d1 + d2] += Fraction(1, 36)
assert sum(pmf.values()) == 1                         # a valid PMF sums to exactly 1

# Empirical PMF: relative frequencies in 600 simulated rolls
rng = np.random.default_rng(7)
rolls = rng.integers(1, 7, size=(600, 2)).sum(axis=1)
values, counts = np.unique(rolls, return_counts=True)
emp = dict(zip(values.tolist(), (counts / counts.sum()).tolist()))

print(" s   p(s)   exact  empirical")
for s in range(2, 13):
    bar = "#" * (3 * int(pmf[s] * 36))                # 3 characters per 1/36 of mass
    print(f"{s:2d}   {str(pmf[s]):5s}  {float(pmf[s]):.3f}  {emp.get(s, 0):.3f}    {bar}")

# Probabilities of events are sums of the PMF
print("P(S is even)  =", sum(p for s, p in pmf.items() if s % 2 == 0))
print("P(S >= 10)    =", sum(p for s, p in pmf.items() if s >= 10))

 s   p(s)   exact  empirical
 2   1/36   0.028  0.033    ###
 3   1/18   0.056  0.055    ######
 4   1/12   0.083  0.082    #########
 5   1/9    0.111  0.102    ############
 6   5/36   0.139  0.143    ###############
 7   1/6    0.167  0.140    ##################
 8   5/36   0.139  0.140    ###############
 9   1/9    0.111  0.117    ############
10   1/12   0.083  0.098    #########
11   1/18   0.056  0.062    ######
12   1/36   0.028  0.028    ###
P(S is even)  = 1/2
P(S >= 10)    = 1/6

Key takeaways

• The PMF $p_X(x) = P(X = x)$ lists the probability of each value of a discrete random variable.
• A valid PMF is non-negative and sums to 1; model formulas often need a normalizing constant.
• The probability of any event about $X$ is the sum of the PMF over the values in that event.
• The empirical PMF (relative frequencies) estimates the true PMF; it is noisy for small samples and gives 0 to unseen values.

🧪 Practice

1. A loaded die has $p(k) = c \cdot k$ for $k = 1, \dots, 6$. Find $c$, then $P(\text{roll} \ge 5)$.
2. Let $X$ be the number of tosses of a fair coin until the first head. Write down $p_X(k)$, verify that it sums to 1, and compute $P(X \text{ is even})$.
3. Change the simulation to 60 rolls and then 60,000 rolls. Which values of $s$ have the largest relative error between the empirical and the exact PMF for small samples, and why?
4. Interview: You have rand5(), which returns a uniform integer from 1 to 5. Implement rand7(), uniform on 1 to 7. Hint: two calls give 25 equally likely pairs; work out the PMF of a value computed from the pair, keep only a range of outcomes whose size is a multiple of 7, and retry otherwise (rejection).

Cumulative Distribution Function

Why it matters. Many questions are about thresholds rather than exact values: "at most 3 failures", "a latency below 200 ms", "at least 10 points". The cumulative distribution function answers all of them, and, unlike the PMF, it is defined for every kind of random variable: discrete, continuous or mixed. It is the universal way to describe a distribution (Section 3.5 showed that a CDF determines the distribution completely).

Definition. The cumulative distribution function (CDF) of a random variable $X$ is

$$F_X(x) = P(X \le x), \qquad x \in \mathbb{R}.$$

For a discrete random variable, it accumulates the PMF from the left:

$$F_X(x) = \sum_{k \le x} p_X(k).$$

Properties. Every CDF:

1. is non-decreasing: $x < y \Rightarrow F(x) \le F(y)$;
2. tends to 0 as $x \to -\infty$ and to 1 as $x \to +\infty$;
3. is right-continuous: $F(x) = \lim_{h \downarrow 0} F(x + h)$.

Conversely, any function with these three properties is the CDF of some random variable. Right-continuity comes from the "$\le$" in the definition: the value at a jump point belongs to the upper step.

For a discrete random variable, the CDF is a step function that jumps by $p(x)$ at each possible value $x$ and is flat in between. For a fair die:

  F(x)
   1 |                                    *=======
 5/6 |                              *=====o
 4/6 |                        *=====o
 3/6 |                  *=====o
 2/6 |            *=====o
 1/6 |      *=====o
   0 |======o-----+-----+-----+-----+-----+-------> x
            1     2     3     4     5     6
     * = value included at the jump (right-continuity), o = value excluded

Computing probabilities from the CDF.

Probability	Formula	Discrete pitfall
$P(X \le b)$	$F(b)$
$P(X > a)$	$1 - F(a)$
$P(a < X \le b)$	$F(b) - F(a)$
$P(X = x)$	$F(x) - F(x^-)$, the size of the jump at $x$	0 at points where $F$ is continuous
$P(X \ge a)$ for integer-valued $X$	$1 - F(a - 1)$	$1 - F(a)$ would wrongly exclude $a$ itself
$P(a \le X \le b)$ for integer-valued $X$	$F(b) - F(a - 1)$	$F(b) - F(a)$ would wrongly exclude $a$

Here $F(x^-) = \lim_{h \downarrow 0} F(x - h)$ is the left limit. The function $S(x) = 1 - F(x) = P(X > x)$ is the survival function; it is the natural object in reliability and survival analysis (Section 12.3).

import numpy as np

# PMF of the sum of two dice: p(s) = (6 - |s - 7|) / 36 for s = 2, ..., 12
support = np.arange(2, 13)
pmf = (6 - np.abs(support - 7)) / 36
cdf = np.cumsum(pmf)                                  # F evaluated at each support point

def F(x):
    """CDF of the sum, defined for EVERY real x, not only at integers."""
    i = np.searchsorted(support, x, side="right")     # number of support points <= x
    return 0.0 if i == 0 else float(cdf[i - 1])

for x in [1.9, 2, 6.5, 7, 7.99, 12, 100]:
    print(f"F({x:>4}) = {F(x):.4f}")

# Interval probabilities: < and <= are NOT interchangeable for a discrete variable
print(f"P(4 <  S <= 7) = F(7) - F(4)       = {F(7) - F(4):.4f}")
print(f"P(4 <= S <= 7) = F(7) - F(3)       = {F(7) - F(3):.4f}")
print(f"P(S >= 10)     = 1 - F(9)          = {1 - F(9):.4f}")
print(f"P(S >= 10)     = 1 - F(10) (wrong) = {1 - F(10):.4f}")
print(f"P(S = 7)       = F(7) - F(7^-)     = {F(7) - F(7 - 1e-9):.4f}")   # size of the jump at 7

F( 1.9) = 0.0000
F(   2) = 0.0278
F( 6.5) = 0.4167
F(   7) = 0.5833
F(7.99) = 0.5833
F(  12) = 1.0000
F( 100) = 1.0000
P(4 <  S <= 7) = F(7) - F(4)       = 0.4167
P(4 <= S <= 7) = F(7) - F(3)       = 0.5000
P(S >= 10)     = 1 - F(9)          = 0.1667
P(S >= 10)     = 1 - F(10) (wrong) = 0.0833
P(S = 7)       = F(7) - F(7^-)     = 0.1667

Key takeaways

• $F_X(x) = P(X \le x)$ is defined for every real $x$ and every type of random variable.
• A CDF is non-decreasing, right-continuous, and goes from 0 to 1.
• For a discrete variable, the CDF is a step function whose jump at $x$ equals $p(x)$.
• $P(a < X \le b) = F(b) - F(a)$; for integer-valued variables, watch the endpoints: $P(X \ge a) = 1 - F(a - 1)$.

🧪 Practice

1. Let $X$ be the number of heads in 3 fair tosses. Write $F_X(x)$ as a piecewise function for all real $x$ and sketch it.
2. A function is defined by $F(x) = 0$ for $x < 0$, $F(x) = 0.3$ for $0 \le x < 2$, $F(x) = 0.8$ for $2 \le x < 5$, and $F(x) = 1$ for $x \ge 5$. Recover the PMF, and compute $P(0 < X < 5)$ and $P(X \ge 2)$.
3. Explain why $G(x) = 0.3$ for $x \le 0$ and $G(x) = 1$ for $x > 0$ is not a CDF, and which property fails. What would fix it?
4. Interview: A dashboard reports "P(at least 3 retries)" as 1 - cdf(3), using a CDF of the number of retries per request. What is wrong, and how do you test such a function to catch this class of bug? Hint: the number of retries is an integer; compare 1 - cdf(3) with $P(\text{retries} > 3)$, and test the edge values with a tiny hand-computable distribution.

Support of a Distribution

Why it matters. Before choosing a model, you should know which values are even possible: counts cannot be negative, proportions live in $[0, 1]$, a binomial count cannot exceed the number of trials. The set of possible values is the support. A model whose support does not cover the data assigns probability zero to something that happened, which makes its likelihood zero and breaks estimation and model comparison. Matching the support to the domain is the first step of every modeling decision.

Definition. For a discrete random variable, the support is the set of values with positive probability:

$$\operatorname{supp}(X) = {x : p_X(x) > 0}.$$

For a continuous random variable with density $f$, the support is (the closure of) the set where $f(x) > 0$. In general, the support is the smallest closed set $S$ with $P(X \in S) = 1$. The values outside the support have no probability mass and no density.

Supports of common distributions (Chapter 5 covers each of them).

Distribution	Support	Typical use
Bernoulli$(p)$	${0, 1}$	a single success or failure
Binomial$(n, p)$	${0, 1, \dots, n}$	successes in $n$ trials
Geometric$(p)$	${1, 2, 3, \dots}$	trials until the first success
Poisson$(\lambda)$	${0, 1, 2, \dots}$	counts with no fixed upper limit
Uniform$(a, b)$	$[a, b]$	a value equally likely anywhere in a range
Exponential$(\lambda)$	$[0, \infty)$	waiting times
Beta$(\alpha, \beta)$	$[0, 1]$	proportions and probabilities
Normal$(\mu, \sigma^2)$	$(-\infty, \infty)$	measurement errors, sums of many effects

Support is a modeling choice. The support of a model encodes hard constraints:

• Too narrow: a Binomial$(10, p)$ model for "support tickets per hour" says 11 tickets is impossible. If 11 ever occurs, the model's likelihood is exactly 0 no matter what $p$ is.
• Too wide: a normal model for delivery times gives positive probability to negative times. This is often harmless if that probability is negligible, but it can produce nonsense predictions in the tails.
• Bounded vs. unbounded: whether the largest possible value is fixed (number of seats in a plane) or open-ended (number of visitors) decides between binomial-type and Poisson-type models.

A subtlety for discrete variables. The support of a discrete variable can be a countable set that is dense. Enumerate the rationals in $[0, 1]$ as $q_1, q_2, \dots$ and set $P(X = q_n) = 2^{-n}$. Every value has positive probability only at rationals, yet the smallest closed set of probability 1 is the whole interval $[0, 1]$. For discrete variables, it is usually more useful to work with the set of atoms ${x : p(x) > 0}$.

import math

def binom_logpmf(k, n=10, p=0.3):
    if not (k == int(k) and 0 <= k <= n):            # support {0, 1, ..., n}
        return -math.inf                              # log(0): an impossible value
    return math.log(math.comb(n, k)) + k * math.log(p) + (n - k) * math.log(1 - p)

def poisson_logpmf(k, lam=3.0):
    if not (k == int(k) and k >= 0):                  # support {0, 1, 2, ...}
        return -math.inf
    return -lam + k * math.log(lam) - math.lgamma(k + 1)

# Support tickets per hour; one busy hour had 11 tickets
data = [2, 4, 3, 1, 5, 3, 2, 11]

for name, logpmf in [("Binomial(10, 0.3)", binom_logpmf), ("Poisson(3)", poisson_logpmf)]:
    ll = sum(logpmf(k) for k in data)                 # log-likelihood of the whole sample
    print(f"{name:18s} log-likelihood = {ll:8.2f}")

# No choice of p rescues the binomial: 11 is outside its support for every p
print("best binomial over a grid of p:", max(sum(binom_logpmf(k, 10, p) for k in data)
                                              for p in [i / 100 for i in range(1, 100)]))

Binomial(10, 0.3)  log-likelihood =     -inf
Poisson(3)         log-likelihood =   -20.38
best binomial over a grid of p: -inf

Key takeaways

• The support is the set of values a random variable can actually take (positive mass or positive density).
• Choose a model whose support matches the domain: counts, bounded counts, positive reals, proportions, all reals.
• A single observation outside a model's support makes its likelihood 0 (log-likelihood $-\infty$), regardless of the parameters.
• Supports that are too wide produce impossible predictions, such as negative durations.

🧪 Practice

1. Give the support of: the number of sixes in 20 die rolls; the number of rolls until the first six; the proportion of sixes in 20 rolls; the waiting time until a bus that arrives uniformly within the next 10 minutes.
2. For each quantity, propose a distribution from the support table and justify it from the support alone: daily rainfall in mm, number of goals in a football match, click-through rate of an ad, the error of a scale's reading.
3. A model for exam scores out of 100 is Normal$(70, 15^2)$. Using statistics.NormalDist, compute the probability it assigns to scores below 0 and above 100. Is the support mismatch a practical problem here?
4. Interview: A forecasting model for daily order counts outputs a normal predictive distribution, and the operations team sees "negative orders" in the 5% scenarios for small stores. How would you fix the model? Hint: the support should be ${0, 1, 2, \dots}$; consider count distributions (Poisson, negative binomial), modeling on a log scale, or truncation, and explain why rounding and clipping at 0 distorts the probabilities.

Functions of a Discrete Random Variable

Why it matters. We rarely stop at the raw random variable. A store cares about sales $= \min(\text{demand}, \text{stock})$, a casino about the payoff of a bet, an engineer about the squared error. Each is a function $Y = g(X)$ of a random variable, and it is itself a random variable whose distribution we need. For discrete variables, finding it requires nothing more than bookkeeping.

The rule. If $X$ is discrete and $Y = g(X)$, then $Y$ is discrete and

$$p_Y(y) = P(g(X) = y) = \sum_{x ,:, g(x) = y} p_X(x).$$

In words: send every value $x$ to $g(x)$, and add up the probabilities of all $x$ that land on the same $y$.

• If $g$ is one-to-one, the probabilities are just relabeled: $p_Y(g(x)) = p_X(x)$.
• If $g$ is many-to-one, the probabilities of the merged values add up.

   X uniform on {-2, -1, 0, 1, 2}           Y = X^2

     x     p_X(x)                            y     p_Y(y)
    -2      1/5   ----------------\
                                   >-------> 4      2/5
     2      1/5   ----------------/
    -1      1/5   ----------------\
                                   >-------> 1      2/5
     1      1/5   ----------------/
     0      1/5   -------------------------> 0      1/5

Applied example: the newsvendor. A bakery stocks 3 units of a perishable cake each day. Daily demand $D$ has a known PMF. The quantities the bakery cares about are all functions of $D$:

• sales $= \min(D, 3)$;
• leftover units $= \max(3 - D, 0)$, thrown away at closing;
• lost sales $= \max(D - 3, 0)$, customers who were turned away.

Each one merges several demand values into a single value (every $D \ge 3$ gives sales of exactly 3), which is precisely the many-to-one rule.

A warning for later. The distribution of $g(X)$ is not obtained by applying $g$ to a summary of $X$. In particular, the average of $g(X)$ is generally not $g$ applied to the average of $X$. Section 4.3 makes this precise (the law of the unconscious statistician), and Section 6.5 quantifies it (Jensen's inequality).

from collections import defaultdict
from fractions import Fraction as F

def pmf_of(g, pmf):
    """PMF of Y = g(X): add up the probabilities of all x that map to the same y."""
    out = defaultdict(F)
    for x, p in pmf.items():
        out[g(x)] += p
    return {y: str(p) for y, p in sorted(out.items())}

# X uniform on {-2, ..., 2}: a many-to-one and a one-to-one function
pX = {x: F(1, 5) for x in range(-2, 3)}
print("Y = X^2    :", pmf_of(lambda x: x * x, pX))
print("Y = 2X + 1 :", pmf_of(lambda x: 2 * x + 1, pX))

# Newsvendor: daily demand D for a cake, 3 units stocked
pD = {0: F(1, 10), 1: F(2, 10), 2: F(3, 10), 3: F(2, 10), 4: F(1, 10), 5: F(1, 10)}
stock = 3
print("sales      :", pmf_of(lambda d: min(d, stock), pD))
print("leftover   :", pmf_of(lambda d: max(stock - d, 0), pD))
print("lost sales :", pmf_of(lambda d: max(d - stock, 0), pD))

Y = X^2    : {0: '1/5', 1: '2/5', 4: '2/5'}
Y = 2X + 1 : {-3: '1/5', -1: '1/5', 1: '1/5', 3: '1/5', 5: '1/5'}
sales      : {0: '1/10', 1: '1/5', 2: '3/10', 3: '2/5'}
leftover   : {0: '2/5', 1: '3/10', 2: '1/5', 3: '1/10'}
lost sales : {0: '4/5', 1: '1/10', 2: '1/10'}

Key takeaways

• A function of a random variable is a random variable.
• For discrete $X$: $p_Y(y) = \sum_{x : g(x) = y} p_X(x)$. One-to-one functions relabel values; many-to-one functions merge them and add their probabilities.
• Many business quantities are functions of an uncertain input (sales, leftovers, payoffs, errors), and their distributions follow from this rule.
• Summaries do not commute with functions in general: the mean of $g(X)$ is usually not $g(\text{mean of } X)$.

🧪 Practice

1. Let $X$ be a fair die roll. Find the PMF of $Y = (X - 3)^2$ and of $Z = X \bmod 3$.
2. Using the bakery's demand PMF, find the PMF of daily profit if each cake costs 3 to make, sells for 10, and leftover cakes are worth nothing. Which stock level, from 0 to 5, makes a loss least likely?
3. Show that if $g$ is one-to-one, then $X$ and $g(X)$ carry the same information: knowing the value of one determines the other. Give an example where $g(X)$ has strictly less information than $X$.
4. Interview: A game pays $$2^K$, where $K$ is the number of tosses of a fair coin up to and including the first head. Find the PMF of the payoff. How much would you pay to play once? Hint: the payoff PMF follows directly from that of $K$ (a one-to-one function); the answer to the second question is the St. Petersburg paradox, which Section 4.3 revisits through expected values.

<a id="42-continuous-random-variables"></a>

4.2 Continuous Random Variables

When a random variable can take any value in a continuum, each single value has probability zero, so probability must be described by densities and areas instead of lists of point masses.

Probability Density Function

Why it matters. Measure a request's latency with infinite precision. The probability that it is exactly 123.456789... ms is zero, and the same holds for every other value; Section 3.5 showed why this is not a contradiction. A PMF would be zero everywhere and say nothing. What does carry information is how concentrated the probability is around each value: latencies near 120 ms may be far more common than latencies near 900 ms. The probability density function captures that concentration, and it is how every continuous model (normal, exponential, gamma, beta, ...) is written.

Definition. A random variable $X$ is continuous (more precisely, absolutely continuous) if there is a function $f_X \ge 0$, its probability density function (PDF), such that

$$P(a \le X \le b) = \int_a^b f_X(x),dx \quad \text{for all } a \le b.$$

A function $f$ is a valid density if and only if $f(x) \ge 0$ everywhere and $\int_{-\infty}^{\infty} f(x),dx = 1$.

What a density value means. For a small width $h$,

$$P(x < X \le x + h) \approx f_X(x),h, \qquad \text{so} \qquad f_X(x) = \lim_{h \to 0} \frac{P(x < X \le x + h)}{h}.$$

A density is probability per unit of $x$, not a probability. Consequences:

• $f(x)$ can be larger than 1. Uniform$(0, 0.25)$ has density 4 on its support, because its probability is squeezed into an interval of length 0.25.
• $f(x)$ has units. If $x$ is in seconds, $f$ is in "per second"; converting $x$ to milliseconds divides the density by 1000.
• $P(X = x) = 0$ for every $x$, so $P(a \le X \le b) = P(a < X < b)$: endpoints never matter.
• Changing $f$ at finitely (or countably) many points changes no probability, so a density is defined only "almost everywhere" (Section 3.5).

Analogy. A thin metal rod has a mass density in grams per centimetre. A single point of the rod has no mass, but every segment has a mass equal to the integral of the density over it. A dense part of the rod is like a range of likely values. Probability density works the same way, with total mass 1.

Discrete vs. continuous at a glance.

	Discrete	Continuous
Describing function	PMF $p(x) = P(X = x)$	PDF $f(x)$, with $P(X = x) = 0$
Range of values of the function	$[0, 1]$	$[0, \infty)$, can exceed 1
Normalization	$\sum_x p(x) = 1$	$\int f(x),dx = 1$
$P(a \le X \le b)$	$\sum_{a \le x \le b} p(x)$	$\int_a^b f(x),dx$
Endpoints	matter	do not matter

Finding a normalizing constant. Suppose a model says $f(x) = c,x(1 - x)$ on $[0, 1]$ and 0 elsewhere. Since $\int_0^1 x(1 - x),dx = \tfrac12 - \tfrac13 = \tfrac16$, the only valid choice is $c = 6$. (This is the Beta$(2, 2)$ density of Section 5.2.)

The density as the limit of a histogram. A histogram drawn with density=True (Section 2.1) has bars whose areas sum to 1. Its height over a bin is the fraction of the sample in the bin divided by the bin width: an estimate of the average density over that bin. With more data and narrower bins, it approaches $f$.

import numpy as np

rng = np.random.default_rng(42)
lam = 0.5                                             # rate: 0.5 arrivals per minute
x = rng.exponential(scale=1 / lam, size=1_000_000)    # waiting times, mean 1/lam = 2 minutes
f = lambda t: lam * np.exp(-lam * t)                  # exponential density

# Density = probability per unit length: P(t < X <= t + h) / h approaches f(t) as h shrinks
t = 1.0
for h in [1.0, 0.1, 0.01]:
    p = np.mean((x > t) & (x <= t + h))
    print(f"h = {h:4.2f}: P({t} < X <= {t}+h) / h = {p / h:.4f}    f({t}) = {f(t):.4f}")

# A density histogram estimates the AVERAGE density over each bin
heights, edges = np.histogram(x, bins=[0, 1, 2, 4, 8, 16], density=True)
for a, b, hgt in zip(edges[:-1], edges[1:], heights):
    avg_f = (np.exp(-lam * a) - np.exp(-lam * b)) / (b - a)   # exact P(a < X <= b) / (b - a)
    print(f"bin [{a:4.1f}, {b:4.1f}): histogram height = {hgt:.4f}, exact average density = {avg_f:.4f}")

# Densities can exceed 1: probability squeezed into a short interval
u = rng.uniform(0, 0.25, 1_000_000)
print("Uniform(0, 0.25), estimated density on [0.1, 0.11]:", round(np.mean((u > 0.1) & (u <= 0.11)) / 0.01, 3))

h = 1.00: P(1.0 < X <= 1.0+h) / h = 0.2384    f(1.0) = 0.3033
h = 0.10: P(1.0 < X <= 1.0+h) / h = 0.2938    f(1.0) = 0.3033
h = 0.01: P(1.0 < X <= 1.0+h) / h = 0.2977    f(1.0) = 0.3033
bin [ 0.0,  1.0): histogram height = 0.3941, exact average density = 0.3935
bin [ 1.0,  2.0): histogram height = 0.2385, exact average density = 0.2387
bin [ 2.0,  4.0): histogram height = 0.1161, exact average density = 0.1163
bin [ 4.0,  8.0): histogram height = 0.0293, exact average density = 0.0293
bin [ 8.0, 16.0): histogram height = 0.0022, exact average density = 0.0022
Uniform(0, 0.25), estimated density on [0.1, 0.11]: 4.004

Key takeaways

• A continuous random variable has $P(X = x) = 0$ for every $x$; probabilities come from integrating a density.
• A valid density is non-negative and integrates to 1; it may exceed 1 and it has units of "probability per unit of $x$".
• $f(x),h$ approximates the probability of a short interval of width $h$ near $x$.
• A density histogram is a piecewise-constant estimate of $f$.

🧪 Practice

1. Find $c$ so that $f(x) = c,e^{-2x}$ for $x \ge 0$ (and 0 otherwise) is a density. Then compute $P(X > 1)$.
2. Which are valid densities? (a) $f(x) = 2x$ on $[0, 1]$; (b) $f(x) = 1.5$ on $[0, 1]$; (c) $f(x) = 3$ on $[2, 2\tfrac13]$; (d) $f(x) = \sin x$ on $[0, \pi]$.
3. In the code, replace the bins by 100 equal bins on $[0, 16]$ and compare the histogram heights with $f$ at the bin midpoints. Where is the agreement worst, and why?
4. Interview: A model evaluation script reports that a test example has "likelihood 3.7" under a continuous model, and a colleague says this must be a bug because probabilities cannot exceed 1. Who is right? Hint: separate densities from probabilities, and consider what happens to the reported number if the feature is rescaled from metres to millimetres.

CDF for Continuous Variables

Why it matters. The density describes where probability is concentrated, but most practical questions ask how much probability lies below a threshold: "what fraction of requests finish under 200 ms?", "what fraction of components survive 5 years?". The CDF answers these directly, it is easier to estimate from data than a density (no bin widths to choose), and it connects densities to quantiles and to sampling algorithms later in this chapter.

From density to CDF and back. For a continuous random variable,

$$F_X(x) = P(X \le x) = \int_{-\infty}^{x} f_X(t),dt, \qquad f_X(x) = F_X'(x) ;\text{ wherever } f_X \text{ is continuous}.$$

The CDF is the running area under the density; the density is the slope of the CDF. A steep CDF means a high density, a flat CDF means no probability in that region.

Properties beyond the general ones (Section 4.1):

• $F_X$ is continuous: it has no jumps, because no single value has positive probability.
• Hence $P(X < x) = P(X \le x) = F_X(x)$, and $P(a < X \le b) = F(b) - F(a)$ regardless of whether the endpoints are included.

Worked example: the exponential distribution. With $f(x) = \lambda e^{-\lambda x}$ for $x \ge 0$:

$$F(x) = \int_0^x \lambda e^{-\lambda t},dt = 1 - e^{-\lambda x} \quad (x \ge 0), \qquad F(x) = 0 \quad (x < 0).$$

With $\lambda = 0.5$ per minute, the probability of waiting at most 2 minutes is $1 - e^{-1} \approx 0.632$.

   density f(x) = 0.5 e^(-0.5x)                   CDF F(x) = 1 - e^(-0.5x)

 0.5 |*                                     1.0 |                        . . . *
     | *                                        |              . . * * *
     |  *                                       |        . * *
     |   **                                 0.5 |      *
     |     **                                   |    *
     |       ***                                |  *
     |          *****                           | *
     |               **********                 |*
   0 +---------------------------> x          0 +---------------------------> x
     0     2     4     6     8                  0     2     4     6     8

   slope of the CDF at x  =  height of the density at x
   F(2) = area under the density between 0 and 2 = 0.632

The empirical CDF. From a sample $x_1, \dots, x_n$, the empirical CDF is

$$\hat F_n(x) = \frac{1}{n} #{i : x_i \le x},$$

a step function that jumps by $1/n$ at every data point. The Glivenko-Cantelli theorem says that $\sup_x |\hat F_n(x) - F(x)| \to 0$ almost surely: the empirical CDF converges to the true CDF uniformly. The gap shrinks roughly like $1/\sqrt{n}$, and its distribution is the basis of the Kolmogorov-Smirnov test (Section 9.5).

import numpy as np

lam = 0.5
F = lambda t: 1 - np.exp(-lam * t)                    # CDF of Exponential(rate 0.5)
f = lambda t: lam * np.exp(-lam * t)                  # its density

# The density is the derivative of the CDF (central finite difference)
h = 1e-5
for t in [0.5, 2.0, 5.0]:
    slope = (F(t + h) - F(t - h)) / (2 * h)
    print(f"t = {t}: slope of F = {slope:.5f}, f(t) = {f(t):.5f}")

# The CDF is the running area under the density (trapezoid rule on a fine grid)
grid = np.linspace(0, 2, 20_001)
area = np.sum((f(grid[1:]) + f(grid[:-1])) / 2 * np.diff(grid))
print(f"area under f on [0, 2] = {area:.5f}, F(2) = {F(2):.5f}")

# The empirical CDF converges uniformly to F, with a gap of order 1/sqrt(n)
rng = np.random.default_rng(0)
for n in [10, 100, 1_000, 10_000, 100_000]:
    x = np.sort(rng.exponential(1 / lam, n))
    i = np.arange(1, n + 1)
    # largest vertical distance between the ECDF steps and F: check just before and at each jump
    gap = max(np.max(i / n - F(x)), np.max(F(x) - (i - 1) / n))
    print(f"n = {n:7,d}: sup |ECDF - F| = {gap:.4f}   sqrt(n) * gap = {np.sqrt(n) * gap:.2f}")

t = 0.5: slope of F = 0.38940, f(t) = 0.38940
t = 2.0: slope of F = 0.18394, f(t) = 0.18394
t = 5.0: slope of F = 0.04104, f(t) = 0.04104
area under f on [0, 2] = 0.63212, F(2) = 0.63212
n =      10: sup |ECDF - F| = 0.2232   sqrt(n) * gap = 0.71
n =     100: sup |ECDF - F| = 0.0848   sqrt(n) * gap = 0.85
n =   1,000: sup |ECDF - F| = 0.0277   sqrt(n) * gap = 0.88
n =  10,000: sup |ECDF - F| = 0.0129   sqrt(n) * gap = 1.29
n = 100,000: sup |ECDF - F| = 0.0036   sqrt(n) * gap = 1.15

Key takeaways

• $F(x) = \int_{-\infty}^x f(t),dt$ and $f = F'$: the CDF is accumulated density, the density is the slope of the CDF.
• A continuous CDF has no jumps, so strict and non-strict inequalities give the same probabilities.
• The exponential CDF is $1 - e^{-\lambda x}$ for $x \ge 0$.
• The empirical CDF converges uniformly to the true CDF (Glivenko-Cantelli), with error of order $1/\sqrt{n}$.

🧪 Practice

1. Derive the CDF of $f(x) = 6x(1 - x)$ on $[0, 1]$ and use it to compute $P(0.25 \le X \le 0.75)$.
2. A CDF is $F(x) = 1 - 1/x^2$ for $x \ge 1$ and 0 otherwise. Find the density and $P(X > 3)$.
3. Modify the code to compute the empirical CDF of 50 samples and print it at $t = 1, 2, 4$ next to $F(t)$. Repeat with five different seeds. How large are the typical deviations compared with $1/\sqrt{50}$?
4. Interview: When comparing latency distributions of two releases, why might you prefer plotting empirical CDFs over histograms? Hint: think about choices a histogram forces on you (bin width, bin origin), how each plot shows tail percentiles, and how easily two curves can be overlaid and read off at a threshold.

Probability as Area

Why it matters. For a continuous variable, every probability is an area: the area under the density over the region of interest. Seeing probabilities as areas makes many problems geometric, explains why single points carry no probability (a line has no area), and suggests three ways to compute any probability: exactly via the CDF, numerically by integration, or by simulation.

The principle.

$$P(a \le X \le b) = \int_a^b f(x),dx = F(b) - F(a).$$

      f(x)
       |              _..----.._
       |           .-'##########'-.
       |         .'################ '.
       |       .'##################   '.
       |     .'####################     '.
       |  _-'######################       '-._
       +--------+--------------------+---------------> x
                a                    b

       shaded area = P(a <= X <= b)        total area under the curve = 1

Three ways to compute an area.

Method	How	When to use
Exact	$F(b) - F(a)$ with a closed-form or library CDF	whenever a CDF is available
Numerical integration	approximate $\int_a^b f$ with the trapezoid or Simpson rule	density known, CDF not available
Monte Carlo	simulate many values, count the fraction in $[a, b]$	complicated regions, many dimensions, or only a sampler is available

The normal distribution is the classic case where the density has a simple formula but the CDF has no closed form; it is computed numerically (via the error function), which is what tables and libraries such as statistics.NormalDist do.

The same idea in two dimensions. When two independent uniform variables are involved, probabilities are areas in the unit square (formally, integrals of a joint density, Section 6.1). Two friends each arrive at a café at a uniformly random time between 12:00 and 13:00, and each waits 15 minutes for the other. They meet if $|X - Y| \le 15$. The two corner triangles where they miss each other have legs of $45$ minutes, so

$$P(\text{meet}) = 1 - \left(\frac{45}{60}\right)^2 = \frac{7}{16} = 0.4375.$$

  Y (min)
   60 +----------------------+
      |             .########|      ### : |X - Y| <= 15, they meet
      |           .##########|
      |         .###########.|      blank corners: they miss each other
      |       .###########.  |      each corner is a right triangle with legs 45
   15 |.    .###########.    |
      |##.############.      |
      |###########.          |
    0 +----------------------+--> X (min)
      0  15                  60

import numpy as np
from statistics import NormalDist

X = NormalDist(mu=170, sigma=8)                       # adult heights in cm (illustrative model)
a, b = 160, 185

# 1. Exact: difference of CDF values
exact = X.cdf(b) - X.cdf(a)

# 2. Numerical integration of the density (trapezoid rule)
grid = np.linspace(a, b, 1_001)
dens = np.array([X.pdf(v) for v in grid])
trapezoid = np.sum((dens[1:] + dens[:-1]) / 2 * np.diff(grid))

# 3. Monte Carlo: fraction of simulated heights that fall in [a, b]
rng = np.random.default_rng(1)
heights = rng.normal(170, 8, 1_000_000)
mc = np.mean((heights >= a) & (heights <= b))

print(f"P({a} <= X <= {b}): exact = {exact:.5f}, trapezoid = {trapezoid:.5f}, Monte Carlo = {mc:.5f}")
print(f"P(X = 170) = area of a line = {X.cdf(170) - X.cdf(170):.1f}")

# Two-dimensional area: the meeting problem, by simulation
x, y = rng.uniform(0, 60, (2, 1_000_000))
print(f"P(meet) simulated = {np.mean(np.abs(x - y) <= 15):.4f}, exact = {7 / 16:.4f}")

P(160 <= X <= 185): exact = 0.86395, trapezoid = 0.86395, Monte Carlo = 0.86448
P(X = 170) = area of a line = 0.0
P(meet) simulated = 0.4383, exact = 0.4375

Key takeaways

• For a continuous variable, a probability is the area under the density over the region of interest.
• Single points and other sets of zero length have zero area, hence zero probability.
• Compute areas exactly with a CDF, numerically by integrating the density, or approximately by Monte Carlo.
• With several independent uniform variables, probabilities become areas or volumes of regions in a square or cube.

🧪 Practice

1. For $X \sim$ Exponential$(\lambda = 0.5)$, compute $P(1 \le X \le 3)$ exactly, then by the trapezoid rule with 10 and 1,000 intervals. How does the error depend on the number of intervals?
2. Solve the meeting problem when each friend waits $w$ minutes. For which $w$ is the probability of meeting $1/2$?
3. A stick of length 1 is broken at a uniformly random point. What is the probability that the longer piece is at least twice the length of the shorter one? Draw the relevant region on $[0, 1]$.
4. Interview: Estimate $\pi$ using only a uniform random number generator, and say how many samples you need for 3 correct decimal places. Hint: the quarter disc $x^2 + y^2 \le 1$ has area $\pi/4$ inside the unit square; the fraction of points that land in it is a proportion whose standard error is $\sqrt{p(1 - p)/n}$.

Quantile Function

Why it matters. The CDF answers "what fraction of values lies below $x$?". Many decisions need the reverse: "below which value do 99% of latencies fall?" (a service-level objective), "how much stock covers demand on 95% of days?", "what loss is exceeded only 1% of the time?" (value at risk). These are quantiles. The quantile function is the inverse of the CDF, and it is also the engine of random number generation (Section 4.4).

Definition. For $p \in (0, 1)$, the quantile function of $X$ is

$$Q(p) = F^{-1}(p) = \inf{x \in \mathbb{R} : F(x) \ge p}.$$

If $F$ is continuous and strictly increasing, $Q$ is its ordinary inverse: $Q(p)$ is the unique $x$ with $F(x) = p$. The general definition (the generalized inverse) handles two cases where an ordinary inverse does not exist:

   F with a jump (discrete X)                  F with a flat part (a gap with no probability)

  F |             *---------                  F |                    .-----
    |             |                             |                   /
  p |- - - - - - -|  <- every p inside the    p |- - - - - *-------*      <- F equals p on the whole
    |             |     jump gives the same     |         /|                 flat part; Q(p) is
    |*------------o     Q(p)                    |   _____/ |                 its LEFT end
    +-------------+-----------> x               +----------+-----------> x
                 Q(p)                                     Q(p)

Named quantiles: the median is $Q(0.5)$, the quartiles are $Q(0.25)$ and $Q(0.75)$, and the percentiles of Section 2.4 are the sample versions of $Q(p)$.

Key properties.

• $Q$ is non-decreasing, and $Q(p) \le x \iff p \le F(x)$. This equivalence is what makes inverse transform sampling work (Section 4.4).
• Quantiles commute with increasing functions: if $g$ is increasing, $Q_{g(X)}(p) = g(Q_X(p))$. For example, the median of $e^X$ is $e^{\text{median}(X)}$, and the 99th percentile of latency in milliseconds is 1000 times the 99th percentile in seconds. Means do not have this property.
• Quantiles do not average: the 99th percentile of a combined population is not the average of the groups' 99th percentiles.

Worked example. For Exponential$(\lambda)$, solving $1 - e^{-\lambda x} = p$ gives

$$Q(p) = -\frac{\ln(1 - p)}{\lambda}.$$

With a mean latency of $1/\lambda = 120$ ms, the median is $120 \ln 2 \approx 83$ ms, below the mean because of the long right tail, and the 99th percentile is $120 \ln 100 \approx 553$ ms.

Quantiles in decisions. In the newsvendor problem of Section 4.1, the stock level that maximizes expected profit is a quantile of demand: $Q(r)$, where $r = (\text{price} - \text{cost}) / \text{price}$ is the critical ratio (Section 4.3 verifies this numerically).

import numpy as np

# 1. Closed-form quantile function: exponential latencies with mean 120 ms
mean = 120.0
Q = lambda p: -mean * np.log(1 - p)                   # inverse of F(x) = 1 - exp(-x / mean)
for p in [0.5, 0.9, 0.99, 0.999]:
    print(f"Q({p:<5}) = {Q(p):5.1f} ms")

# 2. Generalized inverse for a discrete variable (a fair die): smallest x with F(x) >= p
support = np.arange(1, 7)
cdf = support / 6                                     # F(1), ..., F(6)
Q_die = lambda p: int(support[np.searchsorted(cdf, p, side="left")])   # first index with cdf >= p
print("die:", "  ".join(f"Q({p:.3f}) = {Q_die(p)}" for p in [0.1, 1 / 6, 0.17, 0.5, 0.51, 0.99]))

# 3. Quantiles commute with increasing functions; means do not
rng = np.random.default_rng(5)
x = rng.normal(0, 1, 1_000_000)
y = np.exp(x)                                         # an increasing transformation
print(f"median(e^X) = {np.median(y):.3f}   e^median(X) = {np.exp(np.median(x)):.3f}")
print(f"mean(e^X)   = {np.mean(y):.3f}   e^mean(X)   = {np.exp(np.mean(x)):.3f}")

# 4. Quantiles do not average: two servers with different latency profiles
fast = rng.exponential(50, 900_000)                   # 90% of traffic
slow = rng.exponential(400, 100_000)                  # 10% of traffic
p99 = lambda v: np.quantile(v, 0.99)
print(f"p99 fast = {p99(fast):.0f}, p99 slow = {p99(slow):.0f}, "
      f"average of p99s = {(p99(fast) + p99(slow)) / 2:.0f}, "
      f"traffic-weighted average = {0.9 * p99(fast) + 0.1 * p99(slow):.0f}, "
      f"true p99 of all traffic = {p99(np.concatenate([fast, slow])):.0f}")

Q(0.5  ) =  83.2 ms
Q(0.9  ) = 276.3 ms
Q(0.99 ) = 552.6 ms
Q(0.999) = 828.9 ms
die: Q(0.100) = 1  Q(0.167) = 1  Q(0.170) = 2  Q(0.500) = 3  Q(0.510) = 4  Q(0.990) = 6
median(e^X) = 1.002   e^median(X) = 1.002
mean(e^X)   = 1.651   e^mean(X)   = 1.001
p99 fast = 230, p99 slow = 1833, average of p99s = 1031, traffic-weighted average = 390, true p99 of all traffic = 919

Key takeaways

• The quantile function $Q(p) = \inf{x : F(x) \ge p}$ inverts the CDF; it is the ordinary inverse when $F$ is continuous and strictly increasing.
• Median, quartiles and percentiles are values of $Q$; the sample percentiles of Section 2.4 estimate them.
• Quantiles commute with increasing transformations, but they cannot be averaged across groups.
• For the exponential distribution, $Q(p) = -\ln(1 - p)/\lambda$.

🧪 Practice

1. Find the quantile function of Uniform$(a, b)$ and of the distribution with CDF $F(x) = x^3$ on $[0, 1]$. Compute the median and the interquartile range of each.
2. For $X$ = the number of heads in 3 fair tosses, compute $Q(p)$ for $p = 0.1, 0.125, 0.5, 0.6, 0.9$ using the generalized inverse.
3. Show that for an exponential distribution, the ratio of the 99th percentile to the median does not depend on $\lambda$. What does this imply about comparing the tail behavior of two exponential services with different means?
4. Interview: A monitoring system stores the p99 latency of each server every minute, and the global dashboard shows the average of those p99 values. What is wrong, and how do production systems compute global percentiles? Hint: quantiles do not average; consider merging raw values, fixed-bucket histograms, or mergeable quantile sketches (t-digest, DDSketch).

Mixed and Singular Distributions

Why it matters. Real data often fit neither the discrete nor the continuous template. A request that times out is logged at exactly the timeout value; a day's rainfall is exactly 0 on dry days and continuous otherwise; an insurance claim is 0 for most policies. These are mixed distributions, and treating them as purely continuous leads to wrong means, wrong quantiles and broken likelihoods. A third, stranger type, the singular distribution, shows that "no jumps in the CDF" does not guarantee a density.

Every distribution is a mixture of three types. By the Lebesgue decomposition theorem, every CDF can be written uniquely as

$$F = \alpha_d,F_d + \alpha_c,F_{c} + \alpha_s,F_s, \qquad \alpha_d + \alpha_c + \alpha_s = 1,$$

where $F_d$ is discrete (pure jumps), $F_c$ is absolutely continuous (has a density) and $F_s$ is singular (continuous but without a density).

Type	CDF looks like	Example
Discrete	pure staircase	die roll, count of defects
Absolutely continuous	smooth, no jumps	latency, height
Mixed (discrete + continuous)	smooth parts with some jumps	timeout-capped latency, rainfall, claims
Singular	continuous, but all growth happens on a set of length 0	the Cantor distribution

Mixed distributions in practice. Section 3.5 built the traffic-light waiting time as a point mass at 0 plus a uniform part. A common source of mixed data is censoring: the recorded value is $W = \min(T, c)$, where $T$ is continuous and $c$ is a cap (a timeout, a sensor's maximum reading, the end of a study). Then

• $W$ has a point mass at $c$ of size $P(T \ge c)$, a jump in its CDF;
• below $c$, $W$ has the same density as $T$;
• the mean of $W$ underestimates the mean of $T$, because large values are replaced by $c$.

For a mixed variable, probabilities are computed piece by piece: sum the point masses in the region, and integrate the density part over it.

   CDF of W = min(T, 5), T exponential with mean 2

  1.0 |                                   *--------------
      |                                   |   jump = P(T >= 5) = 0.082
 0.92 |                           . . . . o
      |                 . . . . '
      |         . . . '
      |     . '
      |  .'
      | '
    0 +-----------------------------------+-------------> w
      0                                   5 (timeout)

Singular distributions: the Cantor distribution. Choose the base-3 digits of a number independently, each equal to 0 or 2 with probability $1/2$. The result $C$ always lies in the Cantor set of Section 3.5, which has total length 0, so $C$ has no density: any density would have to integrate to 1 over a set of length 0. Yet every single point has probability 0, so the CDF has no jumps either. Its CDF, the Cantor function or "devil's staircase", is continuous, rises from 0 to 1, and is flat on every removed middle third: for example $F(x) = 1/2$ on all of $[1/3, 2/3]$.

  F(x)
  1  |                                                     *
     |                                                    *
     |                                                 ***
     |                                                *
 3/4 |                                         *******
     |                                        *
     |                                     ***
     |                                    *
 1/2 |                  ******************        flat on [1/3, 2/3]: no probability there
     |                 *
     |              ***
     |             *
 1/4 |      *******
     |     *
     |  ***
     | *
  0  |*
     +------------------------------------------------------> x
      0                 1/3              2/3               1

Singular distributions rarely appear in applied work. They matter conceptually: a continuous CDF guarantees $P(X = x) = 0$, but only an absolutely continuous CDF guarantees a density.

import numpy as np

rng = np.random.default_rng(11)

# Mixed: request durations T (exponential, mean 2 s), logged with a 5 s timeout
T = rng.exponential(scale=2.0, size=1_000_000)
timeout = 5.0
W = np.minimum(T, timeout)                            # what the logs actually record

print(f"P(W = 5): exact e^(-2.5) = {np.exp(-2.5):.4f}, simulated = {np.mean(W == timeout):.4f}")
print(f"CDF just below / at the timeout: F(4.999) = {np.mean(W <= 4.999):.4f}, F(5) = {np.mean(W <= 5):.4f}")
print(f"mean of logged W = {W.mean():.3f} s (exact {2 * (1 - np.exp(-2.5)):.3f}), true mean of T = 2.000 s")

# Singular: the Cantor distribution, base-3 digits are 0 or 2 with probability 1/2 each
digits = 2 * rng.integers(0, 2, size=(1_000_000, 35))
C = digits @ (3.0 ** -np.arange(1, 36))
for x in [1 / 9, 1 / 4, 1 / 3, 1 / 2, 2 / 3, 3 / 4]:
    print(f"F({x:.4f}) = {np.mean(C <= x):.4f}")
print("fraction of samples in the removed middle third (1/3, 2/3):", np.mean((C > 1 / 3) & (C < 2 / 3)))
print("largest probability of any single value (estimated):", np.unique(C, return_counts=True)[1].max() / len(C))

P(W = 5): exact e^(-2.5) = 0.0821, simulated = 0.0822
CDF just below / at the timeout: F(4.999) = 0.9178, F(5) = 1.0000
mean of logged W = 1.838 s (exact 1.836), true mean of T = 2.000 s
F(0.1111) = 0.2500
F(0.2500) = 0.3333
F(0.3333) = 0.5001
F(0.5000) = 0.5001
F(0.6667) = 0.5001
F(0.7500) = 0.6664
fraction of samples in the removed middle third (1/3, 2/3): 0.0
largest probability of any single value (estimated): 2e-06

Key takeaways

• Every distribution decomposes into a discrete part, an absolutely continuous part and a singular part.
• Mixed distributions (point masses plus a density) are common: censored or capped measurements, zero-inflated amounts.
• Censoring at $c$ creates a point mass $P(T \ge c)$ at $c$ and biases the naive mean downward.
• A singular distribution, such as the Cantor distribution, has a continuous CDF but no density.

🧪 Practice

1. For the timeout example, compute $P(W > 4)$, $P(W \ge 5)$ and the median of $W$. Which of these are affected by censoring?
2. Daily rainfall $R$ is 0 with probability 0.6; on rainy days it is exponential with mean 8 mm. Write the CDF of $R$, compute $P(R \le 5)$, and compute $E[R]$ by combining the two parts.
3. Explain why the Cantor function satisfies $F(x/3) = F(x)/2$ for $x \in [0, 1]$, and use it to compute $F(1/27)$ and $F(1/4)$. Compare with the simulation.
4. Interview: Latency logs are capped at a 30 s timeout, and 2% of requests hit the cap. A report gives the mean latency computed from the logs. Is it biased, and how would you estimate the true mean or median? Hint: the data are right-censored at 30 s; the median is unaffected as long as fewer than half the values are censored, while the mean needs a model of the tail or survival-analysis methods (Section 12.3).

<a id="43-expectation-and-moments"></a>

4.3 Expectation and Moments

A full distribution is often more than we need or can estimate, so we summarize it with a few numbers (its center, its spread and its shape), and we use generating functions to handle all of these summaries at once.

Expected Value

Why it matters. A casino needs to know how much it earns per bet on average, an insurer how much it pays per policy, an engineer how many retries a request needs. The expected value is the probabilistic version of the average: the value that the sample mean of many independent repetitions settles on. It is the basis of pricing, of decision-making under uncertainty, and of nearly every estimator in statistics.

Definition.

$$E[X] = \sum_x x,p_X(x) \quad \text{(discrete)}, \qquad E[X] = \int_{-\infty}^{\infty} x,f_X(x),dx \quad \text{(continuous)}.$$

For a mixed distribution, add the two parts: the sum over the point masses plus the integral over the density part. The expected value is also written $\mu$, $\mu_X$ or $EX$, and called the mean or expectation of $X$.

The definition requires absolute convergence: $E[X]$ exists only if $\sum |x|,p(x) < \infty$ (or $\int |x| f(x),dx < \infty$). Otherwise the result would depend on the order of summation, and we say the mean does not exist.

Three ways to read it.

Reading	Meaning	Example (fair die, $E[X] = 3.5$)
Weighted average	each value weighted by its probability	$\tfrac16(1 + 2 + \dots + 6) = 3.5$
Center of mass	the balance point of the PMF or density on a seesaw	a seesaw with equal weights at 1 to 6 balances at 3.5
Long-run average	the limit of the sample mean of many independent copies (law of large numbers, Section 6.5)	the average of many rolls approaches 3.5

The expected value need not be a possible value (no die shows 3.5), and it is not the "most likely" value (that is the mode).

The tail formula. For a non-negative random variable,

$$E[X] = \int_0^{\infty} P(X > x),dx, \qquad \text{and for } X \in {0, 1, 2, \dots}: \quad E[X] = \sum_{k=1}^{\infty} P(X \ge k).$$

It often avoids computing the PMF or density at all. For the number of tosses $K$ until the first head, $P(K \ge k) = (1/2)^{k-1}$, so $E[K] = 1 + \tfrac12 + \tfrac14 + \dots = 2$.

When the mean fails to exist.

• St. Petersburg game: a payoff of $2^K$ with $P(K = k) = 2^{-k}$ has $E[2^K] = \sum_k 2^k \cdot 2^{-k} = 1 + 1 + \dots = \infty$. The mean is infinite, yet nobody would pay a fortune to play.
• Cauchy distribution: density $f(x) = \frac{1}{\pi(1 + x^2)}$. It is symmetric about 0, but $\int |x| f(x),dx = \infty$, so $E[X]$ does not exist, and the sample mean never settles down no matter how much data is collected.

import numpy as np
from fractions import Fraction

# Discrete: E[X] = sum of x p(x) for a fair die
print("E[die] =", sum(Fraction(x, 6) for x in range(1, 7)))

# Tail formula for X >= 0: E[X] = integral of P(X > x). Exponential with rate 0.5 (mean 2)
lam = 0.5
x = np.linspace(0, 80, 800_001)
survival = np.exp(-lam * x)                           # P(X > x)
print(f"E[exponential] via the tail integral = {np.sum((survival[1:] + survival[:-1]) / 2 * np.diff(x)):.4f}")

# Long-run averages: the die's settle at 3.5, the Cauchy's never settle (no mean exists)
rng = np.random.default_rng(3)
n_max = 1_000_000
die = rng.integers(1, 7, n_max)
cauchy = np.tan(np.pi * (rng.random(n_max) - 0.5))    # standard Cauchy samples (see Section 4.4)
for n in [10, 100, 1_000, 10_000, 100_000, 1_000_000]:
    print(f"n = {n:9,d}: mean of die rolls = {die[:n].mean():.4f}   mean of Cauchy draws = {cauchy[:n].mean():8.3f}")

E[die] = 7/2
E[exponential] via the tail integral = 2.0000
n =        10: mean of die rolls = 2.9000   mean of Cauchy draws =    0.570
n =       100: mean of die rolls = 3.4600   mean of Cauchy draws =   -0.054
n =     1,000: mean of die rolls = 3.4530   mean of Cauchy draws =    0.392
n =    10,000: mean of die rolls = 3.4883   mean of Cauchy draws =    0.866
n =   100,000: mean of die rolls = 3.4974   mean of Cauchy draws =   -7.214
n = 1,000,000: mean of die rolls = 3.5022   mean of Cauchy draws =   -1.091

Key takeaways

• $E[X] = \sum x,p(x)$ or $\int x f(x),dx$: a probability-weighted average, the center of mass of the distribution.
• The sample mean of independent copies converges to $E[X]$ when the mean exists.
• For non-negative variables, $E[X] = \int_0^\infty P(X > x),dx$; for counts, $E[X] = \sum_{k \ge 1} P(X \ge k)$.
• The mean may be infinite (St. Petersburg) or undefined (Cauchy); then sample means do not stabilize.

🧪 Practice

1. A European roulette wheel has 37 pockets (18 red, 18 black, 1 green). A 1-unit bet on red pays +1 if red comes up and -1 otherwise. Compute the expected gain per bet.
2. Use the tail formula to compute the expected number of rolls of a fair die until the first six.
3. Let $X$ have density $f(x) = 2/x^3$ for $x \ge 1$. Compute $E[X]$. Then show that $Y$ with density $1/x^2$ for $x \ge 1$ has $E[Y] = \infty$. Simulate both (Section 4.4 shows how) and compare the running means.
4. Interview: A lottery ticket costs 2 and has an expected payout of 2.50 because of a huge jackpot rollover. Should a rational person buy one ticket? Should a syndicate buy every combination? Hint: expectation is only one summary; consider the probability of winning anything, the variance, the possibility of sharing the jackpot, taxes, and how many tickets other people buy when the jackpot is large.

Law of the Unconscious Statistician

Why it matters. We usually need the average of a function of a random variable: expected sales $E[\min(D, s)]$, expected squared error $E[(X - c)^2]$, expected utility $E[u(X)]$. Section 4.1 showed how to find the distribution of $g(X)$ first, but that is often tedious. The law of the unconscious statistician (LOTUS) says it is unnecessary: you can average $g$ directly against the distribution of $X$.

Statement. For any function $g$ (for which the expectation exists):

$$E[g(X)] = \sum_x g(x),p_X(x) \quad \text{(discrete)}, \qquad E[g(X)] = \int_{-\infty}^{\infty} g(x),f_X(x),dx \quad \text{(continuous)}.$$

The name is a joke: the formula is so natural that people use it "unconsciously", without noticing that it is a theorem rather than the definition. The definition of $E[Y]$ for $Y = g(X)$ uses the PMF of $Y$; LOTUS says that you get the same number by weighting $g(x)$ by the probabilities of $X$.

Why it is true (discrete case). Group the terms by the value of $g$:

$$\sum_x g(x),p_X(x) = \sum_y \sum_{x : g(x) = y} y,p_X(x) = \sum_y y \underbrace{\sum_{x : g(x) = y} p_X(x)}_{p_Y(y)} = E[Y].$$

$E[g(X)]$ is not $g(E[X])$. Applying a function to the mean is a common mistake, known as the flaw of averages:

$g$	Shape	Relationship	Example
linear, $aX + b$	straight line	$E[g(X)] = g(E[X])$	unit conversion
convex, e.g. $x^2$, $e^x$	curves upward	$E[g(X)] \ge g(E[X])$	$E[X^2] \ge (E[X])^2$
concave, e.g. $\min(x, s)$, $\ln x$	curves downward	$E[g(X)] \le g(E[X])$	expected sales below sales at average demand

The inequalities are Jensen's inequality (Section 6.5).

Applied example: choosing a stock level. For the bakery of Section 4.1, with price 10 and unit cost 3, the profit for a stock level $s$ is $g_s(D) = 10 \min(D, s) - 3s$. LOTUS gives the expected profit for each $s$ without computing any profit distribution, and the best $s$ is the 0.7-quantile of demand, as the critical-ratio rule of Section 4.2 predicts ($r = (10 - 3)/10 = 0.7$).

LOTUS and simulation. The continuous form explains Monte Carlo integration (Section 11.3): to estimate $E[g(X)] = \int g(x) f(x),dx$, draw $x_1, \dots, x_n$ from $f$ and average $g(x_i)$.

import numpy as np
from fractions import Fraction as F

pD = {0: F(1, 10), 1: F(2, 10), 2: F(3, 10), 3: F(2, 10), 4: F(1, 10), 5: F(1, 10)}
E = lambda g: sum(g(d) * p for d, p in pD.items())    # LOTUS: E[g(D)] = sum of g(d) p(d)

# Route 1: LOTUS directly. Route 2: build the PMF of Y = min(D, 3) first (Section 4.1)
lotus = E(lambda d: min(d, 3))
pY = {0: F(1, 10), 1: F(2, 10), 2: F(3, 10), 3: F(4, 10)}
print("E[min(D, 3)] via LOTUS =", lotus, "  via the PMF of Y =", sum(y * p for y, p in pY.items()))
print("min(E[D], 3)           =", min(E(lambda d: d), 3), "  <- the flaw of averages")

# Expected profit for each stock level s: price 10, unit cost 3, unsold cakes are lost
price, cost = 10, 3
cdf = 0
for s in range(6):
    cdf += pD[s]
    profit = E(lambda d: price * min(d, s) - cost * s)
    print(f"stock {s}: expected profit = {float(profit):5.2f}   F(s) = {float(cdf):.1f}")

# Continuous LOTUS by Monte Carlo: E[g(U)] ~ average of g over samples of U ~ Uniform(0, 1)
rng = np.random.default_rng(8)
u = rng.random(1_000_000)
print(f"E[U^2]  ~ {np.mean(u**2):.4f} (exact 1/3)     (E[U])^2 = 0.25")
print(f"E[e^U]  ~ {np.mean(np.exp(u)):.4f} (exact e - 1 = {np.e - 1:.4f})   e^E[U] = {np.exp(0.5):.4f}")

E[min(D, 3)] via LOTUS = 2   via the PMF of Y = 2
min(E[D], 3)           = 23/10   <- the flaw of averages
stock 0: expected profit =  0.00   F(s) = 0.1
stock 1: expected profit =  6.00   F(s) = 0.3
stock 2: expected profit = 10.00   F(s) = 0.6
stock 3: expected profit = 11.00   F(s) = 0.8
stock 4: expected profit = 10.00   F(s) = 0.9
stock 5: expected profit =  8.00   F(s) = 1.0
E[U^2]  ~ 0.3336 (exact 1/3)     (E[U])^2 = 0.25
E[e^U]  ~ 1.7187 (exact e - 1 = 1.7183)   e^E[U] = 1.6487

Key takeaways

• LOTUS: $E[g(X)] = \sum g(x),p(x)$ or $\int g(x) f(x),dx$; no need to find the distribution of $g(X)$.
• In general $E[g(X)] \ne g(E[X])$; equality holds for linear $g$. Convex $g$ pushes the average up, concave $g$ pulls it down.
• Expected profit, loss, utility and error are all LOTUS computations.
• Monte Carlo estimation of $E[g(X)]$ averages $g$ over samples of $X$.

🧪 Practice

1. For a fair die, compute $E[X^2]$, $E[1/X]$ and $E[(X - 3.5)^2]$ with LOTUS. Compare $E[1/X]$ with $1/E[X]$.
2. $X \sim$ Uniform$(0, 2)$. Compute $E[X^3]$ and $E[e^{-X}]$ exactly, then estimate them by Monte Carlo.
3. In the bakery example, add a salvage value: unsold cakes are sold to a discount store for 1 each. Recompute the optimal stock and check the critical-ratio rule with $r = (\text{price} - \text{cost}) / (\text{price} - \text{salvage})$.
4. Interview: A capacity planner computes expected revenue as price * min(mean_demand, capacity). Is the forecast too high, too low, or correct? Hint: $\min(\cdot, \text{capacity})$ is concave, so Jensen's inequality says which way the error goes; describe how a simulation over the demand distribution fixes it.

Variance and Standard Deviation

Why it matters. Two suppliers that both deliver in 5 days on average can be very different: one always takes 4 to 6 days, the other takes either 1 day or 9 days. The mean cannot tell them apart; their spread does. Variance and standard deviation quantify how far a random variable typically falls from its mean. They measure risk in finance, noise in measurement, and uncertainty in every estimate (a standard error is a standard deviation).

Definition. With $\mu = E[X]$:

$$\operatorname{Var}(X) = E\left[(X - \mu)^2\right], \qquad \operatorname{SD}(X) = \sqrt{\operatorname{Var}(X)}.$$

Variance is often written $\sigma^2$ and the standard deviation $\sigma$. Variance is in squared units (days$^2$); the standard deviation is back in the units of $X$, which makes it the one to report.

The computational formula. Expanding the square and using linearity (next topic):

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

It is convenient for exact calculations, but on floating-point data with a large mean relative to the spread it suffers from catastrophic cancellation (Section 2.3); numerical software uses the definition or Welford's algorithm instead.

Population vs. sample. $\operatorname{Var}(X)$ is a property of a distribution, a fixed number. The sample variance $s^2$ of Section 2.3 (with $n - 1$ in the denominator) is a statistic computed from data that estimates it. Chapter 8 makes the connection precise.

Standard results.

Distribution	Mean	Variance
Bernoulli$(p)$	$p$	$p(1 - p)$, largest at $p = 1/2$
Fair die	$3.5$	$35/12 \approx 2.917$
Uniform$(a, b)$	$(a + b)/2$	$(b - a)^2 / 12$
Exponential$(\lambda)$	$1/\lambda$	$1/\lambda^2$

For Uniform$(a, b)$: $E[X^2] = \frac{1}{b - a}\int_a^b x^2,dx = \frac{a^2 + ab + b^2}{3}$, and subtracting $\left(\frac{a + b}{2}\right)^2$ leaves $\frac{(b - a)^2}{12}$.

Why squares? Absolute deviations $E|X - \mu|$ (Section 2.3) are more intuitive, but squares have decisive mathematical advantages: the variance of a sum of independent variables is the sum of the variances (next topic), $E[(X - c)^2]$ is minimized exactly at $c = \mu$, and the normal distribution is parameterized by it. The price is sensitivity to extreme values, and some heavy-tailed distributions have infinite variance even though their mean is finite (for example a Pareto distribution with tail index between 1 and 2, Section 5.2).

import numpy as np
from fractions import Fraction as F

def mean_var(pmf):
    mu = sum(x * p for x, p in pmf.items())
    var_def = sum((x - mu) ** 2 * p for x, p in pmf.items())       # E[(X - mu)^2]
    var_short = sum(x * x * p for x, p in pmf.items()) - mu ** 2    # E[X^2] - mu^2
    assert var_def == var_short                                     # identical in exact arithmetic
    return mu, var_def

suppliers = {
    "A": {4: F(1, 3), 5: F(1, 3), 6: F(1, 3)},      # delivery days, supplier A
    "B": {1: F(1, 2), 9: F(1, 2)},                  # delivery days, supplier B
}
for name, pmf in suppliers.items():
    mu, var = mean_var(pmf)
    print(f"supplier {name}: mean = {mu} days, variance = {var} days^2, SD = {float(var) ** 0.5:.3f} days")

# Bernoulli(p): variance p(1 - p), largest at p = 1/2
ps = [0.01, 0.1, 0.3, 0.5, 0.7, 0.99]
print("p        :", "  ".join(f"{p:5.2f}" for p in ps))
print("p(1 - p) :", "  ".join(f"{p * (1 - p):5.3f}" for p in ps))

# Uniform(a, b): variance (b - a)^2 / 12
rng = np.random.default_rng(12)
a, b = 2.0, 8.0
u = rng.uniform(a, b, 1_000_000)
print(f"Uniform({a:g}, {b:g}): simulated variance = {u.var():.4f}, formula (b - a)^2 / 12 = {(b - a) ** 2 / 12:.4f}")

supplier A: mean = 5 days, variance = 2/3 days^2, SD = 0.816 days
supplier B: mean = 5 days, variance = 16 days^2, SD = 4.000 days
p        :  0.01   0.10   0.30   0.50   0.70   0.99
p(1 - p) : 0.010  0.090  0.210  0.250  0.210  0.010
Uniform(2, 8): simulated variance = 2.9981, formula (b - a)^2 / 12 = 3.0000

Key takeaways

• $\operatorname{Var}(X) = E[(X - \mu)^2] = E[X^2] - \mu^2$ measures spread around the mean; $\operatorname{SD}(X)$ is its square root, in the units of $X$.
• Distributions with the same mean can have very different variances, and therefore very different risk.
• Distribution variance is a fixed parameter; the sample variance estimates it.
• Bernoulli: $p(1 - p)$; Uniform$(a, b)$: $(b - a)^2/12$; Exponential: $1/\lambda^2$.

🧪 Practice

1. Compute the variance and standard deviation of the bakery's daily demand $D$ and of its daily sales $\min(D, 3)$. Which is larger, and why?
2. Derive $\operatorname{Var}(X) = 35/12$ for a fair die using $E[X^2] - (E[X])^2$.
3. Show that $E[(X - c)^2] = \operatorname{Var}(X) + (\mu - c)^2$ for any constant $c$, and conclude that the mean is the best constant predictor of $X$ under squared error.
4. Interview: Two trading strategies have the same expected daily return; one has a daily standard deviation of 1%, the other 3%. Which do you prefer, and is the standard deviation enough to decide? Hint: consider risk-adjusted metrics (Sharpe ratio), the effect of volatility on compounded returns ($E[\ln(1 + R)]$ is roughly $\mu - \sigma^2/2$), and what SD misses: skewness, fat tails, drawdowns.

Properties of Expectation and Variance

Why it matters. Few useful random variables have simple distributions, but many are sums: the total of a shopping cart, the number of successes in a series of trials, the error of an average. A handful of algebraic rules let us find their means and variances without ever computing their distributions. Linearity of expectation in particular solves problems that look hopelessly combinatorial in a single line.

The rules. For random variables $X, Y$ and constants $a, b$:

Rule	Formula	Requires
Constant	$E[c] = c$, $\operatorname{Var}(c) = 0$	nothing
Linearity	$E[aX + bY] = a,E[X] + b,E[Y]$	nothing, not even independence
Monotonicity	$X \le Y \Rightarrow E[X] \le E[Y]$	nothing
Shift and scale	$\operatorname{Var}(aX + b) = a^2 \operatorname{Var}(X)$, $\operatorname{SD}(aX + b) = \lvert a \rvert \operatorname{SD}(X)$	nothing
Sum	$\operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y) + 2\operatorname{Cov}(X, Y)$	nothing
Sum of independent	$\operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y)$	uncorrelated (independence suffices)
Product of independent	$E[XY] = E[X],E[Y]$	independence (or uncorrelated)
Zero variance	$\operatorname{Var}(X) = 0 \iff X = \mu$ almost surely	nothing

Here $\operatorname{Cov}(X, Y) = E[(X - \mu_X)(Y - \mu_Y)]$ is the covariance (Section 6.3). Note that $\operatorname{Var}(X - Y) = \operatorname{Var}(X) + \operatorname{Var}(Y)$ for independent $X, Y$: variances add even when the variables are subtracted.

Linearity with indicators. Write a count as a sum of indicators, $X = \sum_i \mathbf{1}_{A_i}$. Since $E[\mathbf{1}_A] = P(A)$,

$$E[X] = \sum_i P(A_i),$$

whether or not the events are independent. Example: $n$ letters are put into $n$ envelopes at random (the matching problem of Section 3.2). Letter $i$ is in its own envelope with probability $1/n$, so the expected number of correct letters is $n \cdot \frac1n = 1$, for every $n$. Computing the full distribution needed inclusion-exclusion; the mean needs one line. (The variance is also 1 for $n \ge 2$, using covariances of the indicators.)

The variance of an average. If $X_1, \dots, X_n$ are independent with variance $\sigma^2$, then

$$\operatorname{Var}(\bar X) = \operatorname{Var}\left(\frac1n \sum_i X_i\right) = \frac{1}{n^2} \cdot n\sigma^2 = \frac{\sigma^2}{n}, \qquad \operatorname{SD}(\bar X) = \frac{\sigma}{\sqrt n}.$$

This $1/\sqrt n$ law is why averaging reduces noise, why diversification reduces risk, and why the standard error of Section 7.3 shrinks with the square root of the sample size: four times the data halves the uncertainty.

Dependence changes everything for variance. For $\operatorname{Var}(X) = \operatorname{Var}(Y) = 1$:

Relationship	$\operatorname{Var}(X + Y)$	Interpretation
$Y$ independent of $X$	2	risks partly cancel
$Y = X$	4	risks reinforce, no diversification
$Y = -X$	0	a perfect hedge

import numpy as np

rng = np.random.default_rng(14)

# Linearity with indicators: fixed points of a random permutation of n = 10 items
n, trials = 10, 200_000
perms = rng.permuted(np.tile(np.arange(n), (trials, 1)), axis=1)   # each row is a random shuffle
fixed = (perms == np.arange(n)).sum(axis=1)                         # sum of indicators 1{pi(i) = i}
print(f"fixed points: mean = {fixed.mean():.3f} (theory 1), variance = {fixed.var():.3f} (theory 1)")

# Variance of a sum depends on the dependence between the terms
x = rng.standard_normal(1_000_000)
for name, y in [("independent", rng.standard_normal(1_000_000)), ("Y = X", x), ("Y = -X", -x)]:
    print(f"{name:12s}: Var(X + Y) = {np.var(x + y):.3f}")

# The average of n independent die rolls has variance sigma^2 / n (sigma^2 = 35/12)
for k in [1, 4, 16, 64]:
    means = rng.integers(1, 7, size=(200_000, k)).mean(axis=1)
    print(f"n = {k:2d}: Var(mean) = {means.var():.4f}, sigma^2 / n = {35 / 12 / k:.4f}, SD = {means.std():.3f}")

fixed points: mean = 1.001 (theory 1), variance = 0.997 (theory 1)
independent : Var(X + Y) = 2.000
Y = X       : Var(X + Y) = 3.992
Y = -X      : Var(X + Y) = 0.000
n =  1: Var(mean) = 2.9218, sigma^2 / n = 2.9167, SD = 1.709
n =  4: Var(mean) = 0.7308, sigma^2 / n = 0.7292, SD = 0.855
n = 16: Var(mean) = 0.1829, sigma^2 / n = 0.1823, SD = 0.428
n = 64: Var(mean) = 0.0457, sigma^2 / n = 0.0456, SD = 0.214

Key takeaways

• Expectation is linear for any random variables, dependent or not.
• Writing a count as a sum of indicators gives its mean as a sum of probabilities.
• $\operatorname{Var}(aX + b) = a^2 \operatorname{Var}(X)$; variances add for independent (or uncorrelated) variables, otherwise a covariance term appears.
• The average of $n$ independent copies has variance $\sigma^2/n$ and standard deviation $\sigma/\sqrt n$.

🧪 Practice

1. A shopping cart contains 3 items with independent prices: mean 10 and SD 2, mean 25 and SD 5, mean 4 and SD 1. Find the mean and SD of the total. Then find the mean and SD of the total if every price is converted to another currency by multiplying by 1.1.
2. Roll a die 60 times. Use indicators to find the expected number of faces that never appear.
3. $X$ and $Y$ each have variance 1 and $\operatorname{Var}(X + Y) = 3$. Find $\operatorname{Cov}(X, Y)$ and $\operatorname{Var}(X - Y)$.
4. Interview: Shuffle $n$ distinct numbers. A position is a record if its number is larger than all numbers before it. What is the expected number of records? Hint: use one indicator per position; by symmetry, the $i$-th number is the largest of the first $i$ with probability $1/i$, and the answer grows like $\ln n$.

Raw and Central Moments

Why it matters. The mean locates a distribution and the variance measures its spread, but two distributions can agree on both and still differ in shape: one may have a long right tail, another may produce extreme values far more often. Moments are the systematic family of summaries that includes the mean and variance and continues to skewness and kurtosis. They connect to the descriptive statistics of Section 2.4 and to the method-of-moments estimation of Section 8.2.

Definitions. For $k = 1, 2, 3, \dots$:

Name	Definition	$k = 1$	$k = 2$
Raw moment	$m_k = E[X^k]$	the mean $\mu$	$E[X^2]$
Central moment	$\mu_k = E[(X - \mu)^k]$	0	the variance $\sigma^2$
Standardized moment	$\tilde\mu_k = E\left[\left(\frac{X - \mu}{\sigma}\right)^k\right] = \mu_k/\sigma^k$	0	1

The third and fourth standardized moments have names:

• Skewness $\tilde\mu_3$: 0 for symmetric distributions, positive for a long right tail, negative for a long left tail.
• Kurtosis $\tilde\mu_4$: measures how heavy the tails are. The normal distribution has kurtosis 3, so excess kurtosis $\tilde\mu_4 - 3$ compares tails with the normal's.

Mechanical analogy. Treat the distribution as a mass on a rod. The first moment locates the center of mass; the second central moment is the moment of inertia about it (how hard it is to spin); higher moments describe the asymmetry and the far-flung mass. The vocabulary of moments was borrowed from physics.

From raw to central moments. Expanding $(X - \mu)^k$ with the binomial theorem:

$$\mu_k = \sum_{j=0}^{k} \binom{k}{j},m_j,(-\mu)^{k - j}, \qquad \text{e.g.} \quad \mu_2 = m_2 - \mu^2, \quad \mu_3 = m_3 - 3\mu m_2 + 2\mu^3.$$

Existence. If $E|X|^k < \infty$, all lower moments are finite too. Heavy-tailed distributions have only finitely many moments: a Pareto distribution with tail index $\alpha$ has $E[X^k] < \infty$ only for $k < \alpha$, and Student's $t$ with $\nu$ degrees of freedom only for $k < \nu$ (Chapter 5). When a moment is infinite, the corresponding sample statistic does not converge; it is dominated by the few largest observations.

Do moments determine a distribution? Usually, but not always. If the moment generating function exists (next topic), the sequence of moments determines the distribution. Without it, two different distributions can share every moment; the lognormal distribution is the classic example.

from math import comb, factorial
import numpy as np

def central_from_raw(m, k):
    """k-th central moment from raw moments m[0..k] (with m[0] = 1), via the binomial expansion."""
    mu = m[1]
    return sum(comb(k, j) * m[j] * (-mu) ** (k - j) for j in range(k + 1))

# Exponential(rate 1): raw moments E[X^k] = k!
raw = [factorial(k) for k in range(5)]
central = [central_from_raw(raw, k) for k in range(5)]
print("raw moments m_0..m_4      :", raw)
print("central moments mu_0..mu_4:", central)
print(f"skewness = {central[3] / central[2] ** 1.5:.1f}, kurtosis = {central[4] / central[2] ** 2:.1f} (excess {central[4] / central[2] ** 2 - 3:.1f})")

rng = np.random.default_rng(15)
x = rng.exponential(1.0, 2_000_000)
print("sample raw moments        :", [round(float(np.mean(x ** k)), 3) for k in range(5)])

# Heavy tails: Pareto(alpha = 3) has finite mean and variance, but infinite 3rd and 4th moments
for seed in range(4):
    y = np.random.default_rng(seed).pareto(3, 1_000_000) + 1     # classical Pareto with minimum 1
    z = (y - y.mean()) / y.std()
    print(f"Pareto(3), seed {seed}: mean = {y.mean():.4f} (exact 1.5), "
          f"sample skewness = {np.mean(z ** 3):6.2f}, sample kurtosis = {np.mean(z ** 4):7.2f}")

raw moments m_0..m_4      : [1, 1, 2, 6, 24]
central moments mu_0..mu_4: [1, 0, 1, 2, 9]
skewness = 2.0, kurtosis = 9.0 (excess 6.0)
sample raw moments        : [1.0, 0.999, 1.993, 5.953, 23.675]
Pareto(3), seed 0: mean = 1.5002 (exact 1.5), sample skewness =  10.76, sample kurtosis =  345.13
Pareto(3), seed 1: mean = 1.4982 (exact 1.5), sample skewness = 146.68, sample kurtosis = 72792.55
Pareto(3), seed 2: mean = 1.5005 (exact 1.5), sample skewness =  15.17, sample kurtosis = 1085.39
Pareto(3), seed 3: mean = 1.5003 (exact 1.5), sample skewness =  37.78, sample kurtosis = 8228.22

Key takeaways

• Raw moments $E[X^k]$, central moments $E[(X - \mu)^k]$ and standardized moments $E[((X - \mu)/\sigma)^k]$ describe location, spread and shape.
• Skewness (3rd standardized moment) measures asymmetry; kurtosis (4th) measures tail heaviness, with 3 for the normal.
• Central moments follow from raw moments by the binomial expansion.
• Heavy-tailed distributions lack higher moments; their sample skewness and kurtosis are unstable and meaningless.

🧪 Practice

1. For $X \sim$ Bernoulli$(p)$, show that $E[X^k] = p$ for all $k \ge 1$, and derive the skewness. For which $p$ is it zero?
2. For $X \sim$ Uniform$(0, 1)$, compute $m_k = E[X^k]$ for $k = 1, \dots, 4$, convert to central moments and find the excess kurtosis.
3. Run the Pareto loop with $\alpha = 5$ instead of 3. Do the sample skewness and kurtosis stabilize across seeds? Relate the answer to which moments exist.
4. Interview: A risk report flags a return series as "approximately normal because its sample skewness is 0.1 and its excess kurtosis is 0.5". What would make you distrust this conclusion? Hint: the sampling variability of kurtosis depends on the 8th moment; consider a series with a few rare crashes, the sample size, and complementing moments with tail plots or Q-Q plots (Section 5.3).

Moment Generating Functions

Why it matters. Computing $E[X^k]$ one integral at a time is laborious, and finding the distribution of a sum of independent variables directly requires convolution (Section 6.4). The moment generating function packs all the moments of a distribution into one function, turns sums of independent variables into products, and identifies distributions uniquely. It also produces the sharpest elementary tail bounds (Chernoff bounds), widely used in the analysis of algorithms and machine learning.

Definition. The moment generating function (MGF) of $X$ is

$$M_X(t) = E\left[e^{tX}\right],$$

defined for the values of $t$ where the expectation is finite. The MGF is useful when it is finite on an open interval around $t = 0$.

Why "moment generating". Expanding the exponential,

$$M_X(t) = E\left[\sum_{k=0}^{\infty} \frac{t^k X^k}{k!}\right] = \sum_{k=0}^{\infty} \frac{E[X^k]}{k!},t^k, \qquad \text{so} \qquad E[X^k] = M_X^{(k)}(0).$$

The moments are the Taylor coefficients of $M_X$ at 0: differentiate $k$ times and evaluate at 0.

Properties.

Property	Statement
Normalization	$M_X(0) = 1$
Linear transformation	$M_{aX + b}(t) = e^{bt} M_X(at)$
Independent sums	$M_{X + Y}(t) = M_X(t),M_Y(t)$ if $X$ and $Y$ are independent
Uniqueness	if $M_X = M_Y$ on an interval around 0, then $X \overset{d}{=} Y$

MGFs of common distributions.

Distribution	$M_X(t)$	Domain
Bernoulli$(p)$	$1 - p + p e^t$	all $t$
Binomial$(n, p)$	$(1 - p + p e^t)^n$	all $t$
Poisson$(\lambda)$	$\exp(\lambda(e^t - 1))$	all $t$
Exponential$(\lambda)$	$\frac{\lambda}{\lambda - t}$	$t < \lambda$
Normal$(\mu, \sigma^2)$	$\exp(\mu t + \sigma^2 t^2/2)$	all $t$

Using the product rule and uniqueness together. The binomial MGF is the Bernoulli MGF to the $n$-th power, because a binomial count is a sum of $n$ independent Bernoulli trials. If $X \sim$ Poisson$(\lambda_1)$ and $Y \sim$ Poisson$(\lambda_2)$ are independent, then

$$M_{X + Y}(t) = e^{\lambda_1(e^t - 1)},e^{\lambda_2(e^t - 1)} = e^{(\lambda_1 + \lambda_2)(e^t - 1)},$$

which is the Poisson$(\lambda_1 + \lambda_2)$ MGF, so by uniqueness $X + Y \sim$ Poisson$(\lambda_1 + \lambda_2)$. No convolution sum was needed.

Chernoff bounds. For any $t > 0$, Markov's inequality (Section 6.5) applied to $e^{tX}$ gives

$$P(X \ge a) = P\left(e^{tX} \ge e^{ta}\right) \le e^{-ta} M_X(t), \qquad \text{so} \qquad P(X \ge a) \le \min_{t > 0} e^{-ta} M_X(t).$$

For sums of independent variables, the bound decays exponentially in the number of terms, far faster than Chebyshev's $1/k^2$ (Section 2.3).

Limitations. The MGF does not exist for heavy-tailed distributions: the Cauchy and Pareto distributions have $M_X(t) = \infty$ for every $t > 0$, and so does the lognormal, even though all its moments are finite. The characteristic function (next topic) fixes this.

import numpy as np
from math import comb

# Moments as derivatives at 0: Exponential(rate 1), M(t) = 1 / (1 - t) for t < 1
M = lambda t: 1 / (1 - t)
h = 1e-3
d1 = (M(h) - M(-h)) / (2 * h)                                     # central finite differences
d2 = (M(h) - 2 * M(0) + M(-h)) / h ** 2
d3 = (M(2 * h) - 2 * M(h) + 2 * M(-h) - M(-2 * h)) / (2 * h ** 3)
print(f"M'(0) = {d1:.4f}, M''(0) = {d2:.4f}, M'''(0) = {d3:.4f}   (exact E[X^k] = k!: 1, 2, 6)")

# Independent sums multiply MGFs: Poisson(2) + Poisson(3) has the Poisson(5) MGF
rng = np.random.default_rng(16)
s = rng.poisson(2, 1_000_000) + rng.poisson(3, 1_000_000)
for t in [-1.0, 0.2, 0.5]:
    print(f"t = {t:4.1f}: empirical E[e^(tS)] = {np.mean(np.exp(t * s)):8.4f}   "
          f"Poisson(5) MGF = {np.exp(5 * (np.exp(t) - 1)):8.4f}")

# Tail bounds for X ~ Binomial(100, 1/2): P(X >= 75)
n, a = 100, 75
exact = sum(comb(n, k) for k in range(a, n + 1)) / 2 ** n
chebyshev = (n / 4) / (a - n / 2) ** 2                            # Var / (a - mu)^2
ts = np.linspace(0.001, 3, 30_000)
chernoff = np.min(np.exp(-ts * a) * ((1 + np.exp(ts)) / 2) ** n)  # min over t of e^(-ta) M(t)
print(f"P(X >= 75): exact = {exact:.2e}, Chernoff bound = {chernoff:.2e}, Chebyshev bound = {chebyshev:.2e}")

M'(0) = 1.0000, M''(0) = 2.0000, M'''(0) = 6.0000   (exact E[X^k] = k!: 1, 2, 6)
t = -1.0: empirical E[e^(tS)] =   0.0424   Poisson(5) MGF =   0.0424
t =  0.2: empirical E[e^(tS)] =   3.0245   Poisson(5) MGF =   3.0253
t =  0.5: empirical E[e^(tS)] =  25.6416   Poisson(5) MGF =  25.6260
P(X >= 75): exact = 2.82e-07, Chernoff bound = 2.08e-06, Chebyshev bound = 4.00e-02

Key takeaways

• $M_X(t) = E[e^{tX}]$; its $k$-th derivative at 0 is $E[X^k]$.
• The MGF of a sum of independent variables is the product of their MGFs.
• An MGF that is finite near 0 determines the distribution uniquely, which proves closure results such as Poisson + Poisson = Poisson.
• Chernoff bounds $P(X \ge a) \le \min_t e^{-ta} M_X(t)$ give exponentially small tail bounds.
• Heavy-tailed distributions (Cauchy, Pareto, lognormal) have no MGF.

🧪 Practice

1. Derive the MGF of Bernoulli$(p)$ and use it to get the binomial mean $np$ and variance $np(1 - p)$.
2. Derive $M_X(t) = \lambda/(\lambda - t)$ for the exponential distribution, and use it to find $E[X]$ and $\operatorname{Var}(X)$.
3. Use MGFs to show that if $X \sim$ Normal$(\mu_1, \sigma_1^2)$ and $Y \sim$ Normal$(\mu_2, \sigma_2^2)$ are independent, then $X + Y$ is normal. What are its parameters?
4. Interview: A randomized algorithm succeeds on each independent run with probability 0.6, and you take a majority vote over $n$ runs. How large should $n$ be for the majority to be wrong with probability below $10^{-6}$? Hint: the number of failures is binomial; apply a Chernoff bound (or Hoeffding's inequality, $P(\bar X - p \ge \varepsilon) \le e^{-2n\varepsilon^2}$) and notice that $n$ grows only logarithmically in the target error.

Characteristic Functions

Why it matters. The MGF is powerful but fails for heavy-tailed distributions. The characteristic function keeps every advantage of the MGF (sums become products, it determines the distribution) and always exists. It is the tool behind the standard proof of the central limit theorem (Section 6.5), it explains why averaging Cauchy data does not help, and through the fast Fourier transform it computes distributions of sums efficiently.

Definition. The characteristic function of $X$ is

$$\varphi_X(t) = E\left[e^{itX}\right] = E[\cos(tX)] + i,E[\sin(tX)], \qquad t \in \mathbb{R},$$

where $i = \sqrt{-1}$. Since $|e^{itX}| = 1$, the expectation is always finite. For a continuous variable, $\varphi_X$ is the Fourier transform of the density: it describes the distribution in terms of frequencies instead of locations. Formally, $\varphi_X(t) = M_X(it)$ whenever the MGF exists.

Properties.

Property	Statement
Always exists	$\lvert\varphi_X(t)\rvert \le 1$ and $\varphi_X(0) = 1$
Linear transformation	$\varphi_{aX + b}(t) = e^{ibt}\varphi_X(at)$
Independent sums	$\varphi_{X + Y}(t) = \varphi_X(t),\varphi_Y(t)$
Moments	if $E\lvert X\rvert^k < \infty$, then $\varphi_X^{(k)}(0) = i^k E[X^k]$
Symmetry	$\varphi_X$ is real-valued if and only if $X \overset{d}{=} -X$
Uniqueness and inversion	$\varphi_X$ determines the distribution, and there are explicit inversion formulas
Lévy continuity theorem	$X_n \to X$ in distribution if and only if $\varphi_{X_n}(t) \to \varphi_X(t)$ for every $t$

Examples.

Distribution	$\varphi_X(t)$
Normal$(\mu, \sigma^2)$	$\exp(i\mu t - \sigma^2 t^2/2)$
Exponential$(\lambda)$	$\frac{\lambda}{\lambda - it}$
Poisson$(\lambda)$	$\exp(\lambda(e^{it} - 1))$
Standard Cauchy	$e^{-\lvert t\rvert}$

Why averaging Cauchy data is useless. If $X_1, \dots, X_n$ are independent standard Cauchy, the mean $\bar X$ has characteristic function

$$\varphi_{\bar X}(t) = \left(\varphi_X(t/n)\right)^n = \left(e^{-|t|/n}\right)^n = e^{-|t|}.$$

By uniqueness, the average of $n$ Cauchy variables has the same distribution as a single one. This is exactly the erratic running mean seen in the expected value topic. Note also that $e^{-|t|}$ has a corner at $t = 0$: it is not differentiable there, consistent with the Cauchy distribution having no mean.

Characteristic functions and the FFT. For an integer-valued variable, evaluating $\varphi_X$ at the frequencies $t = 2\pi k/L$ is a discrete Fourier transform of its PMF. To get the exact distribution of a sum of independent variables, transform each PMF with the FFT, multiply, and transform back. This is the same idea as multiplying MGFs, but numerically stable and fast: $O(L \log L)$ instead of repeated convolutions.

import numpy as np

rng = np.random.default_rng(17)

def ecf(sample, t):
    """Empirical characteristic function: the average of exp(i t X) over the sample."""
    return np.mean(np.exp(1j * np.outer(t, sample)), axis=1)

t = np.array([0.5, 1.0, 2.0, 4.0])
normal = rng.standard_normal(200_000)
cauchy = np.tan(np.pi * (rng.random(200_000) - 0.5))
cauchy_mean50 = np.tan(np.pi * (rng.random((200_000, 50)) - 0.5)).mean(axis=1)   # averages of 50 draws

print("t                    :", "  ".join(f"{v:6.2f}" for v in t))
print("normal   ECF (real)  :", "  ".join(f"{v.real:6.3f}" for v in ecf(normal, t)), "  exp(-t^2/2):",
      " ".join(f"{v:.3f}" for v in np.exp(-t ** 2 / 2)))
print("Cauchy   ECF (real)  :", "  ".join(f"{v.real:6.3f}" for v in ecf(cauchy, t)), "  exp(-|t|)  :",
      " ".join(f"{v:.3f}" for v in np.exp(-t)))
print("mean of 50 Cauchy ECF:", "  ".join(f"{v.real:6.3f}" for v in ecf(cauchy_mean50, t)), "  <- unchanged")

# Exact PMF of the sum of 10 dice: multiply characteristic functions (FFT), then invert
die = np.zeros(61)                                    # length L > 60 avoids wrap-around
die[1:7] = 1 / 6
pmf_fft = np.fft.ifft(np.fft.fft(die) ** 10).real     # phi_sum = phi_die^10, inverted

pmf_conv = np.array([1.0])                            # check: 9 repeated convolutions
for _ in range(10):
    pmf_conv = np.convolve(pmf_conv, die[:7])
print(f"P(sum of 10 dice = 35) = {pmf_fft[35]:.6f}, P(sum >= 50) = {pmf_fft[50:].sum():.6f}")
print(f"max |FFT - convolution| = {np.max(np.abs(pmf_fft - pmf_conv[:61])):.1e}")

t                    :   0.50    1.00    2.00    4.00
normal   ECF (real)  :  0.882   0.606   0.135   0.002   exp(-t^2/2): 0.882 0.607 0.135 0.000
Cauchy   ECF (real)  :  0.608   0.369   0.136   0.020   exp(-|t|)  : 0.607 0.368 0.135 0.018
mean of 50 Cauchy ECF:  0.607   0.368   0.134   0.018   <- unchanged
P(sum of 10 dice = 35) = 0.072693, P(sum >= 50) = 0.002890
max |FFT - convolution| = 4.3e-17

Key takeaways

• $\varphi_X(t) = E[e^{itX}]$ always exists and determines the distribution.
• Sums of independent variables multiply characteristic functions; moments are derivatives at 0 when they exist.
• Convergence of characteristic functions is equivalent to convergence in distribution (Lévy), the basis of the central limit theorem proof.
• The average of $n$ standard Cauchy variables is standard Cauchy: $\varphi_{\bar X}(t) = e^{-|t|}$.
• For integer-valued variables, FFTs of PMFs compute exact distributions of sums efficiently.

🧪 Practice

1. Compute the characteristic function of a fair coin $X \in {-1, 1}$ and of the sum of $n$ independent copies. Why is it real-valued?
2. Use $\varphi_X^{(k)}(0) = i^k E[X^k]$ on $\varphi(t) = e^{-t^2/2}$ to find the first four moments of the standard normal.
3. Modify the FFT code to compute the exact PMF of the sum of 100 dice, and compare $P(\text{sum} \le 330)$ with the normal approximation NormalDist(350, (100 * 35 / 12) ** 0.5).cdf(330.5).
4. Interview: How would you compute the exact probability distribution of the total of 10,000 independent integer-valued random variables, each with a different known PMF on ${0, \dots, 20}$? Hint: naive repeated convolution costs $O(\text{total support}^2)$; multiply transforms instead (FFT with enough zero padding, possibly in a divide-and-conquer tree), and watch floating-point error in tiny tail probabilities.

<a id="44-transformations-of-random-variables"></a>

4.4 Transformations of Random Variables

Applying a function to a random variable produces a new random variable, and this section gives the tools for finding its distribution, culminating in the result that turns uniform random numbers into samples from any distribution.

Linear Transformations

Why it matters. Changing units (Celsius to Fahrenheit, seconds to milliseconds), adding a fixed fee, and standardizing a feature before training a model are all linear transformations $Y = aX + b$. They are the simplest transformations and the most common, and they come with exact rules: the shape of the distribution is preserved, while its location and scale change in predictable ways.

Rules. For $Y = aX + b$ with $a \ne 0$:

Quantity	Rule
Mean	$E[Y] = a,E[X] + b$
Variance, SD	$\operatorname{Var}(Y) = a^2 \operatorname{Var}(X)$, $\operatorname{SD}(Y) = \lvert a\rvert \operatorname{SD}(X)$
Quantiles	$Q_Y(p) = a,Q_X(p) + b$ if $a > 0$; $Q_Y(p) = a,Q_X(1 - p) + b$ if $a < 0$
CDF	$F_Y(y) = F_X\left(\frac{y - b}{a}\right)$ if $a > 0$; $F_Y(y) = 1 - F_X\left(\frac{y - b}{a}\right)$ if $a < 0$ and $X$ is continuous
Density	$f_Y(y) = \frac{1}{\lvert a\rvert} f_X\left(\frac{y - b}{a}\right)$
Skewness	unchanged if $a > 0$; changes sign if $a < 0$
Kurtosis	unchanged

The density rule has a factor $1/|a|$: stretching the axis by a factor $a$ spreads the same probability over a region $|a|$ times longer, so the density must drop by that factor for the total area to remain 1.

   X ~ Uniform(0, 1)                     Y = 3X + 2 ~ Uniform(2, 5)

  f |                                   f |
  1 +--------+                            |
    |########|                       1/3  |       +-----------------+
    |########|                            |       |#################|
    +--------+-------> x                  +-------+-----------------+----> y
    0        1                            0       2                 5

     area = 1 x 1 = 1                      area = 3 x 1/3 = 1

Location-scale families. Many named families are closed under linear transformations: if $Z$ has a fixed "standard" distribution, then $X = \mu + \sigma Z$ ranges over the whole family. Normal, uniform, Cauchy and logistic are location-scale families; the exponential and gamma (for a fixed shape) are scale families. This is why a single standard normal table suffices for every normal distribution (Section 5.3).

Standardization. The linear transformation $Z = (X - \mu)/\sigma$ has mean 0 and standard deviation 1; it is the probabilistic version of the z-score of Section 2.4. Standardization changes location and scale only. It does not change the shape: a skewed variable remains exactly as skewed after standardization.

import numpy as np

rng = np.random.default_rng(18)

def summary(v):
    z = (v - v.mean()) / v.std()
    return f"mean = {v.mean():8.3f}, SD = {v.std():7.3f}, skewness = {np.mean(z ** 3):6.3f}, kurtosis = {np.mean(z ** 4):6.3f}"

x = rng.exponential(scale=10, size=1_000_000)         # processing time in minutes (mean 10)
print("X              :", summary(x))
print("Y = 60X + 30   :", summary(60 * x + 30))         # seconds, plus a 30 s setup time
print("-X             :", summary(-x))                  # a negative scale flips the skewness
print("Z = (X - 10)/10:", summary((x - 10) / 10))       # standardized: same shape as X

# Density rule: f_Y(y) = f_X((y - b) / a) / |a|, checked with a narrow-bin density histogram of Y = 60X + 30
f_X = lambda t: np.exp(-t / 10) / 10
y = 60 * x + 30
counts, edges = np.histogram(y, bins=np.arange(30, 1236, 6))   # bins of width 6 s
heights = counts / (len(y) * 6)                       # divide by ALL samples: some Y lie beyond the last bin
for v in [60, 210, 480, 930]:
    i = (v - 30) // 6                                 # index of the bin [v, v + 6)
    mid = edges[i] + 3
    print(f"y = {mid:5.0f}: histogram = {heights[i]:.6f}, f_X((y - 30)/60) / 60 = {f_X((mid - 30) / 60) / 60:.6f}")

X              : mean =   10.004, SD =  10.013, skewness =  2.018, kurtosis =  9.197
Y = 60X + 30   : mean =  630.268, SD = 600.755, skewness =  2.018, kurtosis =  9.197
-X             : mean =  -10.004, SD =  10.013, skewness = -2.018, kurtosis =  9.197
Z = (X - 10)/10: mean =    0.000, SD =   1.001, skewness =  2.018, kurtosis =  9.197
y =    63: histogram = 0.001567, f_X((y - 30)/60) / 60 = 0.001577
y =   213: histogram = 0.001227, f_X((y - 30)/60) / 60 = 0.001229
y =   483: histogram = 0.000786, f_X((y - 30)/60) / 60 = 0.000783
y =   933: histogram = 0.000374, f_X((y - 30)/60) / 60 = 0.000370

Key takeaways

• For $Y = aX + b$: the mean and quantiles transform like the values, the SD scales by $|a|$, and the variance by $a^2$.
• The density transforms as $f_Y(y) = f_X((y - b)/a)/|a|$; the factor $1/|a|$ keeps the total area at 1.
• Linear transformations preserve the shape of a distribution (skewness up to sign, and kurtosis).
• Standardizing produces mean 0 and SD 1 but does not make a variable normal.

🧪 Practice

1. Daily temperatures in Celsius have mean 18 and SD 6. Find the mean, SD and variance in Fahrenheit ($F = 1.8C + 32$). If the 90th percentile is 26 C, what is it in Fahrenheit?
2. $X \sim$ Exponential with mean 1. Find the CDF and density of $Y = 5X + 2$, and state its support.
3. Prove $F_Y(y) = 1 - F_X((y - b)/a)$ for $a < 0$ and continuous $X$. Where does the argument fail for discrete $X$?
4. Interview: A feature pipeline standardizes all inputs, and a junior engineer concludes that the model's normality assumptions now hold. What would you explain? Hint: standardization is linear, so it preserves skewness and kurtosis; normality needs a nonlinear transformation (log, Box-Cox, quantile transform), and the model's assumptions might concern the errors, not the features.

CDF Method

Why it matters. Most interesting transformations are not linear: the square of an error, the maximum of several latencies, the logarithm of a price. For a nonlinear $g$, there is no shortcut rule like $aX + b$, but there is one method that always works, for any $g$, monotone or not: compute the CDF of $Y = g(X)$ directly from its definition, then differentiate to get the density if needed.

The recipe. To find the distribution of $Y = g(X)$ for continuous $X$:

1. Write the CDF of $Y$ as a probability about $X$: $F_Y(y) = P(g(X) \le y)$.
2. Solve the inequality $g(X) \le y$ for $X$, giving a set $A_y = {x : g(x) \le y}$.
3. Compute $F_Y(y) = P(X \in A_y)$ using the CDF or density of $X$.
4. Differentiate: $f_Y(y) = F_Y'(y)$.
5. State the support of $Y$.

Example 1: a square. Let $X \sim$ Uniform$(-1, 1)$ and $Y = X^2$. The map is two-to-one, but the CDF method does not care:

$$F_Y(y) = P(X^2 \le y) = P(-\sqrt y \le X \le \sqrt y) = \frac{2\sqrt y}{2} = \sqrt y, \qquad f_Y(y) = \frac{1}{2\sqrt y}, \qquad 0 < y < 1.$$

     g(x) = x^2
        |
     1  |\                       /
        | \                     /
     y  |- *-------------------* - -     {X^2 <= y} = {-sqrt(y) <= X <= sqrt(y)}
        |   \                 /
        |    `-.           .-'
        |       `--.___.--'
        +---+----+-----+-----+---+-----> x
           -1  -sqrt(y)  0  sqrt(y)  1
                 |<--------->|
                  this interval of X maps below y

The density $1/(2\sqrt y)$ is unbounded near 0: squaring pushes values toward 0, piling probability up there.

Example 2: a squared standard normal. For $Z \sim$ Normal$(0, 1)$ and $W = Z^2$, the same steps give $F_W(w) = P(-\sqrt w \le Z \le \sqrt w) = 2\Phi(\sqrt w) - 1$, where $\Phi$ is the standard normal CDF. This is the chi-square distribution with 1 degree of freedom (Section 5.4), the building block of variance estimates and chi-square tests.

Example 3: a maximum. Let $M = \max(X_1, \dots, X_n)$ for independent $X_i$ with common CDF $F$. The maximum is at most $m$ exactly when every $X_i$ is:

$$F_M(m) = P(X_1 \le m, \dots, X_n \le m) = F(m)^n.$$

For uniforms, $F_M(m) = m^n$. This is the first result about order statistics (Section 6.4), and it explains tail latency at scale: if a page waits for 100 independent backend calls, each slower than 1 second with probability 1%, the page is slower than 1 second with probability $1 - 0.99^{100} \approx 63%$.

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(19)
N = 1_000_000
Phi = NormalDist().cdf

y = rng.uniform(-1, 1, N) ** 2                        # Y = X^2, X ~ Uniform(-1, 1):   F_Y(y) = sqrt(y)
w = rng.standard_normal(N) ** 2                       # W = Z^2, Z ~ Normal(0, 1):     F_W(w) = 2 Phi(sqrt w) - 1
m = rng.random((N, 10)).max(axis=1)                   # M = max of 10 Uniform(0, 1):   F_M(m) = m^10

print("   v    P(Y<=v)  sqrt(v)   P(W<=v)  2Phi(sqrt v)-1   P(M<=v)   v^10")
for v in [0.05, 0.25, 0.5, 0.9]:
    print(f"{v:5.2f}   {np.mean(y <= v):.4f}   {np.sqrt(v):.4f}    {np.mean(w <= v):.4f}   {2 * Phi(np.sqrt(v)) - 1:.4f}"
          f"           {np.mean(m <= v):.4f}   {v ** 10:.4f}")

# Tail at scale: a page waits for the slowest of k backend calls, each slow (> 1 s) with probability 1%
for k in [1, 10, 100]:
    print(f"{k:3d} calls: P(page slower than 1 s) = 1 - 0.99^{k} = {1 - 0.99 ** k:.3f}")

   v    P(Y<=v)  sqrt(v)   P(W<=v)  2Phi(sqrt v)-1   P(M<=v)   v^10
 0.05   0.2235   0.2236    0.1775   0.1769           0.0000   0.0000
 0.25   0.4997   0.5000    0.3834   0.3829           0.0000   0.0000
 0.50   0.7067   0.7071    0.5207   0.5205           0.0010   0.0010
 0.90   0.9486   0.9487    0.6571   0.6572           0.3486   0.3487
  1 calls: P(page slower than 1 s) = 1 - 0.99^1 = 0.010
 10 calls: P(page slower than 1 s) = 1 - 0.99^10 = 0.096
100 calls: P(page slower than 1 s) = 1 - 0.99^100 = 0.634

Key takeaways

• CDF method: $F_Y(y) = P(g(X) \le y)$, rewrite as an event about $X$, compute, then differentiate.
• It works for any transformation, including non-monotone ones like $X^2$.
• The square of a standard normal is chi-square with 1 degree of freedom: $F(w) = 2\Phi(\sqrt w) - 1$.
• The maximum of $n$ independent variables with CDF $F$ has CDF $F^n$, which drives tail latency in fan-out systems.

🧪 Practice

1. $X \sim$ Uniform$(0, 1)$. Use the CDF method to find the distribution of $Y = -\ln X$ and of $Y = \sqrt X$.
2. Find the CDF and density of the minimum of $n$ independent Uniform$(0, 1)$ variables. Hint: the minimum exceeds $m$ exactly when all of them do.
3. $X \sim$ Exponential$(\lambda)$. Find the distribution of $Y = \lfloor X \rfloor$ (a discrete variable) using the CDF of $X$, and identify it.
4. Interview: A service fans out each user request to 50 shards and returns when all of them respond. Each shard has a p99 latency of 200 ms. Estimate the fraction of user requests slower than 200 ms, and propose two engineering mitigations. Hint: use $F_M = F^n$ with $F(200) = 0.99$; look up "hedged requests" and the trade-offs of returning partial results.

Change of Variables Formula

Why it matters. When $g$ is monotone and differentiable, the CDF method always ends with the same calculation. Doing it once in general gives a formula for the density of $Y = g(X)$ that can be applied directly. It is how the lognormal and many other distributions are derived from simpler ones, and its multivariate version, with a Jacobian determinant, underlies the transformations of joint densities in Section 6.1 and generative models such as normalizing flows.

Statement. Let $X$ have density $f_X$, and let $g$ be strictly monotone and differentiable on the support of $X$, with inverse $h = g^{-1}$. Then $Y = g(X)$ has density

$$f_Y(y) = f_X\big(h(y)\big),\left|h'(y)\right| = \frac{f_X(x)}{|g'(x)|}\bigg|_{x = h(y)}.$$

Derivation for increasing $g$: $F_Y(y) = P(g(X) \le y) = P(X \le h(y)) = F_X(h(y))$; differentiate with the chain rule. For decreasing $g$, the inequality flips and the derivative picks up a minus sign, which the absolute value absorbs.

Intuition: probability is conserved. A short interval $[x, x + dx]$ is mapped to an interval $[y, y + dy]$ with $dy \approx |g'(x)|,dx$. Both intervals contain the same probability, so

$$f_Y(y),|dy| = f_X(x),|dx| \quad \Longrightarrow \quad f_Y(y) = f_X(x) \left|\frac{dx}{dy}\right|.$$

    x-axis    |--[x, x+dx]--|                       probability f_X(x) dx
                     |
                     |  g stretches by |g'(x)|
                     v
    y-axis    |------[ y, y+dy ]------|             same probability, spread over |g'(x)| dx
                                                    => density divided by |g'(x)|

Where $g$ stretches the axis, density thins out; where $g$ compresses it, density piles up. The factor $|h'(y)|$ is the one-dimensional Jacobian. Forgetting it is the most common error in transformations: the result does not integrate to 1.

Non-monotone $g$. If $g$ is many-to-one, sum over all preimages: $f_Y(y) = \sum_{x : g(x) = y} f_X(x)/|g'(x)|$. For $Y = X^2$, the two preimages $\pm\sqrt y$ each contribute, which reproduces the CDF method's result.

Example: the lognormal distribution. Let $X \sim$ Normal$(\mu, \sigma^2)$ and $Y = e^X$. Then $h(y) = \ln y$ and $h'(y) = 1/y$, so

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

This is the lognormal density (Section 5.2), a standard model for incomes, file sizes and prices. Its summaries show how differently they behave under a nonlinear map:

Summary of $Y = e^X$	Value	Equals $e^{(\text{same summary of } X)}$?
Median	$e^{\mu}$	yes: quantiles commute with increasing maps (Section 4.2)
Mean	$e^{\mu + \sigma^2/2}$	no: larger, by Jensen's inequality
Mode	$e^{\mu - \sigma^2}$	no: smaller, because the Jacobian $1/y$ shifts the peak left

import numpy as np
from statistics import NormalDist

mu, sigma = 1.0, 0.6
phi = lambda x: np.exp(-(x - mu) ** 2 / (2 * sigma ** 2)) / (sigma * np.sqrt(2 * np.pi))   # Normal(mu, sigma^2) density
f_Y = lambda y: phi(np.log(y)) / y                    # change of variables: f_X(ln y) * |d ln y / dy|
f_wrong = lambda y: phi(np.log(y))                    # the same formula WITHOUT the Jacobian 1/y

# Only the correct density integrates to 1 (trapezoid rule on a fine grid)
grid = np.linspace(1e-6, 200, 2_000_001)
for name, dens in [("with Jacobian   ", f_Y), ("without Jacobian", f_wrong)]:
    vals = dens(grid)
    print(f"{name}: integral over (0, 200] = {np.sum((vals[1:] + vals[:-1]) / 2 * np.diff(grid)):.4f}")

# Compare with simulation: density histogram, median and mean
rng = np.random.default_rng(20)
y = np.exp(rng.normal(mu, sigma, 2_000_000))
counts, edges = np.histogram(y, bins=np.arange(0, 12.25, 0.25))
heights = counts / (len(y) * 0.25)                    # normalize by ALL samples, including those above 12
for v in [1.0, 2.0, 4.0, 8.0]:
    i = np.searchsorted(edges, v) - 1                 # bin [v - 0.25, v)
    print(f"y = {v}: histogram = {heights[i]:.4f}, f_Y at the bin midpoint = {f_Y(v - 0.125):.4f}")

print(f"median: simulated {np.median(y):.3f}, e^mu = {np.exp(mu):.3f}")
print(f"mean  : simulated {y.mean():.3f}, e^(mu + sigma^2/2) = {np.exp(mu + sigma ** 2 / 2):.3f}")
print(f"mode  : argmax of f_Y = {grid[np.argmax(f_Y(grid))]:.3f}, e^(mu - sigma^2) = {np.exp(mu - sigma ** 2):.3f}, "
      f"e^(mode of X) = {np.exp(mu):.3f}")

with Jacobian   : integral over (0, 200] = 1.0000
without Jacobian: integral over (0, 200] = 3.2544
y = 1.0: histogram = 0.1276, f_Y at the bin midpoint = 0.1276
y = 2.0: histogram = 0.2929, f_Y at the bin midpoint = 0.2928
y = 4.0: histogram = 0.1451, f_Y at the bin midpoint = 0.1441
y = 8.0: histogram = 0.0174, f_Y at the bin midpoint = 0.0175
median: simulated 2.717, e^mu = 2.718
mean  : simulated 3.252, e^(mu + sigma^2/2) = 3.254
mode  : argmax of f_Y = 1.897, e^(mu - sigma^2) = 1.896, e^(mode of X) = 2.718

The "without Jacobian" integral is not arbitrary: substituting $y = e^x$ shows that it equals $E[e^X] = e^{\mu + \sigma^2/2} \approx 3.254$, the mean of $Y$, not the total probability.

Key takeaways

• For monotone differentiable $g$ with inverse $h$: $f_Y(y) = f_X(h(y)),|h'(y)|$.
• The factor $|h'(y)|$ (the Jacobian) accounts for stretching and compression of the axis; without it, the result is not a density.
• For many-to-one $g$, sum the contributions of all preimages.
• Exponentiating a normal gives the lognormal; its median transforms simply, but its mean and mode do not.

🧪 Practice

1. $X \sim$ Uniform$(0, 1)$. Use the change of variables formula to find the density of $Y = X^3$ and of $Y = 1/X$ (state the supports).
2. $X \sim$ Exponential$(1)$. Find the density of $Y = e^{-X}$ and identify the distribution.
3. Derive the density of $Y = X^2$ for $X \sim$ Normal$(0, 1)$ by summing over both preimages, and check that it matches the derivative of $2\Phi(\sqrt y) - 1$ from the CDF method.
4. Interview: In a normalizing flow, a neural network maps a simple latent variable $z$ (standard normal) to a data point $x = g(z)$ through invertible layers. The training loss contains a "log-determinant" term. Where does it come from, and why must the layers be designed so that it is cheap to compute? Hint: write the multivariate change of variables formula $\log p_X(x) = \log p_Z(h(x)) + \log|\det J_h(x)|$; a general $d \times d$ determinant costs $O(d^3)$, while triangular Jacobians cost $O(d)$.

Probability Integral Transform

Why it matters. Every continuous random variable, however skewed or heavy-tailed, becomes uniform when passed through its own CDF. This simple fact has wide-reaching uses: it is why p-values are uniformly distributed when the null hypothesis is true, how forecasters check whether predicted distributions are calibrated, how copulas separate dependence from marginal distributions (Section 6.3), and, run in reverse, how computers generate random numbers from any distribution (next topic).

Theorem. If $X$ is continuous with CDF $F$, then

$$U = F(X) \sim \text{Uniform}(0, 1).$$

Proof (for strictly increasing $F$). For $u \in (0, 1)$:

$$P(U \le u) = P(F(X) \le u) = P(X \le F^{-1}(u)) = F(F^{-1}(u)) = u,$$

which is the Uniform$(0, 1)$ CDF. (The general case uses the generalized inverse of Section 4.2.)

Intuition. $F(X)$ answers "what fraction of the distribution lies below the value we observed?" That is, it gives the percentile rank of the observation within its own distribution. A random draw is equally likely to land in any percentile: it is below the 10th percentile 10% of the time, below the median 50% of the time, and so on. So its percentile rank is uniform.

Discrete variables are an exception. If $X$ has atoms, $F(X)$ takes only finitely or countably many values and cannot be uniform. The randomized PIT fixes this by spreading each atom uniformly over its jump: $U = F(X^-) + V,p(X)$ with an independent $V \sim$ Uniform$(0, 1)$.

Calibration of probabilistic forecasts. Suppose a model issues, for each case, a predictive CDF $\hat F_i$, and the actual outcome is $y_i$. If the forecasts are correct, the values $u_i = \hat F_i(y_i)$ are uniform. A histogram of the $u_i$ (the PIT histogram) diagnoses miscalibration:

PIT histogram shape	Diagnosis
flat	calibrated
U-shaped (too many values near 0 and 1)	overconfident: predictive distributions too narrow
hump-shaped (too many values near 0.5)	underconfident: predictive distributions too wide
sloped (too many near 0 or near 1)	biased: forecasts systematically too high or too low

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

rng = np.random.default_rng(21)
N = 100_000

def pit_hist(u, bins=10):
    """Heights of a 10-bin density histogram of u on [0, 1]; all 1.0 means uniform."""
    counts, _ = np.histogram(u, bins=bins, range=(0, 1))
    return "  ".join(f"{c:.2f}" for c in counts / len(u) * bins)

# 1. Any continuous X through its own CDF is uniform: exponential with mean 2
x = rng.exponential(2.0, N)
print(f"{'exponential, F(X)':25s}:", pit_hist(1 - np.exp(-x / 2)))

# 2. Calibration check: outcomes are N(0, 2^2); the forecaster announces a normal predictive CDF
y = rng.normal(0, 2, N)
for name, mu, sd in [("calibrated, N(0, 2^2)", 0, 2), ("overconfident, N(0, 1)", 0, 1),
                     ("underconfident, N(0, 4^2)", 0, 4), ("biased, N(1, 2^2)", 1, 2)]:
    cdf = NormalDist(mu, sd).cdf
    print(f"{name:25s}:", pit_hist(np.array([cdf(v) for v in y])))

# 3. Discrete X: F(X) is not uniform; the randomized PIT is
lam = 2.0
k = rng.poisson(lam, N)
pmf = np.array([exp(-lam) * lam ** j / factorial(j) for j in range(40)])
cdf = np.cumsum(pmf)
print(f"{'Poisson, F(X)':25s}:", pit_hist(cdf[k]))
u_rand = cdf[k] - pmf[k] + rng.random(N) * pmf[k]     # F(X^-) + V p(X)
print(f"{'Poisson, randomized PIT':25s}:", pit_hist(u_rand))

exponential, F(X)        : 0.99  0.99  0.99  1.01  0.99  0.99  0.99  1.01  1.02  1.02
calibrated, N(0, 2^2)    : 1.00  1.00  0.99  1.00  1.00  1.00  1.01  0.99  1.01  1.00
overconfident, N(0, 1)   : 2.60  0.76  0.59  0.53  0.51  0.51  0.53  0.60  0.76  2.61
underconfident, N(0, 4^2): 0.05  0.41  1.01  1.58  1.94  1.95  1.59  1.01  0.41  0.05
biased, N(1, 2^2)        : 2.17  1.49  1.23  1.08  0.94  0.83  0.73  0.64  0.53  0.37
Poisson, F(X)            : 0.00  1.34  0.00  0.00  2.71  0.00  2.69  0.00  1.79  1.46
Poisson, randomized PIT  : 0.99  1.00  1.02  0.98  1.01  1.00  0.98  0.99  1.01  1.01

Key takeaways

• For continuous $X$ with CDF $F$, the value $F(X)$ is Uniform$(0, 1)$: every percentile rank is equally likely.
• It fails for discrete variables; the randomized PIT restores uniformity.
• PIT histograms diagnose probabilistic forecasts: flat means calibrated, U-shaped means overconfident, a hump means underconfident, a slope means bias.
• The same fact makes p-values uniform under a true null hypothesis (Chapter 9) and underlies copulas and goodness-of-fit tests.

🧪 Practice

1. $X$ has CDF $F(x) = x^2$ on $[0, 1]$. Simulate $X$ (for example as the maximum of two uniforms) and verify that $X^2$ is uniform.
2. Show that if $U = F(X)$ is uniform, then $1 - U$ is also uniform. What does this say about the PIT histogram of a forecaster who reports $1 - \hat F$ by mistake?
3. In the calibration example, use a forecaster with the correct mean but standard deviation 1.8. How visible is the miscalibration in the PIT histogram with $N = 1{,}000$ and with $N = 100{,}000$?
4. Interview: A weather model outputs a full predictive distribution for tomorrow's temperature at 500 stations every day. How would you check whether these distributions are reliable, beyond checking the accuracy of the mean forecast? Hint: compute $\hat F(y_{\text{observed}})$ for every forecast and look at the PIT histogram, overall and per station or season; complement it with a proper scoring rule such as CRPS, since a flat PIT histogram alone does not guarantee sharp forecasts.

Inverse Transform Sampling

Why it matters. Computers produce uniform random numbers natively, but simulations need samples from exponential, Pareto, empirical or custom distributions. Inverse transform sampling converts uniform numbers into samples from any distribution whose quantile function can be evaluated. It is the principle behind library samplers, behind weighted random choice (load balancers, recommendation sampling, A/B traffic splits), and behind the reproducible simulations used throughout this guide.

Theorem. If $U \sim$ Uniform$(0, 1)$ and $Q$ is the quantile function of a distribution with CDF $F$, then

$$X = Q(U) \quad \text{has CDF } F.$$

Proof. By the key property of the generalized inverse (Section 4.2), $Q(u) \le x \iff u \le F(x)$. Hence $P(X \le x) = P(Q(U) \le x) = P(U \le F(x)) = F(x)$. The result holds for every distribution: discrete, continuous or mixed.

This is the probability integral transform run backwards: $F$ maps the distribution to uniform, and $Q$ maps uniform back to the distribution.

   F(x)
    1 |                         ___-------
      |                     _.-'
    u |- - - - - - - - - -*'                1. draw u uniformly on the vertical axis
      |                 .'|                 2. move right until you hit the CDF
      |              .-'  |                 3. drop down: x = Q(u) is the sample
      |          _.-'     |
      |   __..--'         |                 Steep parts of F (high density) capture many u
    0 +-------------------+--------> x      values; flat parts capture few.
                        x = Q(u)

The algorithm.

1. Find the quantile function $Q$: solve $F(x) = u$ for $x$, or use the generalized inverse.
2. Draw $u \sim$ Uniform$(0, 1)$.
3. Return $Q(u)$.

Common cases.

Distribution	CDF $F(x)$	Sampler $Q(u)$
Uniform$(a, b)$	$\frac{x - a}{b - a}$	$a + (b - a)u$
Exponential$(\lambda)$	$1 - e^{-\lambda x}$	$-\ln(1 - u)/\lambda$
Pareto$(x_m, \alpha)$	$1 - (x_m/x)^{\alpha}$	$x_m (1 - u)^{-1/\alpha}$
Cauchy	$\tfrac12 + \tfrac1\pi \arctan x$	$\tan(\pi(u - \tfrac12))$ (used in Section 4.3)
Discrete with PMF $p_1, \dots, p_k$	cumulative sums $c_j$	smallest $j$ with $c_j \ge u$ (binary search)
Empirical distribution of data	ECDF	the sample quantile at $u$

When $Q$ is not available in closed form. The normal CDF cannot be inverted with elementary functions. Options: numerical inversion (root finding or accurate rational approximations, which is how NormalDist.inv_cdf works), special-purpose methods (Box-Muller, ziggurat), or rejection sampling and MCMC (Chapter 11), which do not need $F$ at all.

Discrete sampling cost. With cumulative sums and binary search, each sample costs $O(\log k)$ for $k$ categories after $O(k)$ setup. When millions of samples are drawn from a fixed PMF, the alias method reduces this to $O(1)$ per sample.

import numpy as np

rng = np.random.default_rng(22)
u = rng.random(1_000_000)                             # the only source of randomness below

# 1. Exponential(rate 0.5) via Q(u) = -ln(1 - u) / lambda
lam = 0.5
x = -np.log1p(-u) / lam                               # log1p(-u) = ln(1 - u), accurate for small u
print(f"exponential: mean = {x.mean():.4f} (exact 2), variance = {x.var():.4f} (exact 4)")

# 2. Pareto(x_m = 1, alpha = 2.5) via Q(u) = (1 - u)^(-1/alpha): a heavy tail from uniform numbers
alpha = 2.5
p = (1 - u) ** (-1 / alpha)
print(f"Pareto: P(X > 10) = {np.mean(p > 10):.5f} (exact {10 ** -alpha:.5f}), max of 1e6 draws = {p.max():.0f}")

# 3. Discrete: weighted random choice by binary search in the cumulative sums
weights = np.array([5, 1, 3, 1])                      # traffic split across four backends
cdf = np.cumsum(weights) / weights.sum()
cdf[-1] = 1.0                                         # guard against round-off leaving cdf[-1] < 1
idx = np.searchsorted(cdf, u, side="left")            # smallest j with cdf[j] >= u
print("backend frequencies:", np.round(np.bincount(idx) / len(idx), 4), " target:", weights / weights.sum())

# 4. Mixed: traffic-light wait of Section 3.5 (0 with prob 0.4, else Uniform(0, 60])
#    F(w) = 0.4 + 0.6 w / 60 on [0, 60], so the generalized inverse is 0 for u <= 0.4
w = np.where(u <= 0.4, 0.0, 60 * (u - 0.4) / 0.6)
print(f"traffic light: P(W = 0) = {np.mean(w == 0):.4f} (exact 0.4), mean = {w.mean():.3f} (exact 18)")

exponential: mean = 1.9977 (exact 2), variance = 3.9969 (exact 4)
Pareto: P(X > 10) = 0.00315 (exact 0.00316), max of 1e6 draws = 470
backend frequencies: [0.5    0.1003 0.3    0.0997]  target: [0.5 0.1 0.3 0.1]
traffic light: P(W = 0) = 0.4002 (exact 0.4), mean = 17.981 (exact 18)

Key takeaways

• If $U$ is uniform, $Q(U)$ has the distribution with quantile function $Q$; this works for discrete, continuous and mixed distributions.
• Closed-form samplers exist for the uniform, exponential, Pareto, Cauchy and others; discrete distributions use cumulative sums and binary search.
• When $Q$ has no closed form, use numerical inversion, special-purpose algorithms, or rejection and MCMC methods.
• Inverse transform sampling is the probability integral transform run backwards.

🧪 Practice

1. Derive an inverse transform sampler for the density $f(x) = 3x^2$ on $[0, 1]$ and verify it by comparing the sample mean with $E[X] = 3/4$.
2. Write a sampler for the geometric distribution $P(K = k) = (1 - p)^{k - 1} p$ using a closed-form generalized inverse (no loop or search). Hint: $P(K > k) = (1 - p)^k$.
3. Sample from the empirical distribution of a dataset by inverse transform, and explain why this is the same as drawing data points uniformly with replacement (the bootstrap of Section 8.4).
4. Interview: Implement weighted_choice(items, weights) returning one item with probability proportional to its weight. Then explain how to support millions of calls per second, and what changes if weights are updated frequently. Hint: prefix sums plus binary search give $O(\log n)$ per call; the alias method gives $O(1)$ for static weights; a Fenwick (binary indexed) tree supports $O(\log n)$ updates and sampling.

<a id="5-common-probability-distributions"></a>