Lecture 6: Continuous Random Variables, PDFs, CDFs

Why do we need a new kind of random variable?

Every random variable so far has been discrete: a finite or countably infinite list of possible values, each with its own probability mass. But plenty of quantities we’d like to model — a waiting time, a position, a measurement error — can take any real value in an interval. There’s no way to list them out one at a time.

Definition: Continuous random variable, PDF

A random variable \(X\) is continuous if there is a function \(f_X(x)\ge 0\), called its probability density function (PDF), such that for every interval \([a,b]\), \[ P(a \le X \le b) = \int_a^b f_X(x)\, dx, \qquad \int_{-\infty}^{\infty} f_X(x)\,dx = 1. \]

---

Instead of bars of probability sitting on individual points (a PMF), probability is now area under a curve.

Example: the continuous uniform random variable

Definition: Uniform(\(a,b\))

\[ f_X(x) = \begin{cases}\dfrac{1}{b-a} & a\le x\le b\\[4pt] 0 & \text{otherwise.}\end{cases} \]

Every value in \([a,b]\) is “equally likely,” in the sense that the density is flat — probability of a sub-interval depends only on its length, not its location.

---

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

a_u, b_u = 2, 5

def f_uniform(x):
    x = np.asarray(x, dtype=float)
    return np.where((x >= a_u) & (x <= b_u), 1 / (b_u - a_u), 0.0)

xs = np.linspace(a_u - 1, b_u + 1, 400)
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(xs, f_uniform(xs), color='black')
ax.set_ylim(0, 1 / (b_u - a_u) * 1.3)
ax.set_xlabel('$x$')
ax.set_ylabel('$f_X(x)$')
ax.set_title(f'Uniform({a_u},{b_u}) density')
plt.tight_layout()
plt.show()

total, _ = quad(f_uniform, a_u - 1, b_u + 1)
p_3_4, _ = quad(f_uniform, 3, 4)
print(f"total area:        {total:.4f}")
print(f"P(3 <= X <= 4) =   {p_3_4:.4f}   (should be 1/3 = {1/3:.4f}, since [3,4] is 1/3 of [2,5])")

total area:        1.0000
P(3 <= X <= 4) =   0.3333   (should be 1/3 = 0.3333, since [3,4] is 1/3 of [2,5])

Example: a PDF can take arbitrarily large values

A density is not a probability

Nothing in the definition requires \(f_X(x)\le 1\) — only that the total area equals \(1\). Squeeze the same unit of probability into a short enough interval and the height of the density can be as large as you like.

Fact: \(P(X=x_0)=0\) for continuous \(X\)

For any single point \(x_0\), \[ P(X=x_0) = \int_{x_0}^{x_0} f_X(x)\,dx = 0, \] since we’re integrating over an interval of zero width. Continuous random variables put no probability mass on individual points — only on intervals of positive length.

Expectation of a continuous random variable

Definition: Expectation

\[ E[X] = \int_{-\infty}^{\infty} x\, f_X(x)\, dx. \]

The rules we’ve seen previously generalize once we replace summations with integrals!

Theorem: LOTUS for continuous random variables

For any function \(g\), \[ E[g(X)] = \int_{-\infty}^{\infty} g(x)\, f_X(x)\, dx. \]

---

Definitions: moments and variance

The \(k\)-th moment of \(X\) is \(E[X^k]\). Variance is defined exactly as in the discrete case, \[ \mathrm{Var}(X) = E\big[(X-E[X])^2\big] = E[X^2]-(E[X])^2, \] and the shortcut formula’s proof used only linearity — nothing about summing versus integrating — so it carries over unchanged.

In-class exercise

Prove that for constants \(a,b\) and a continuous random variable \(X\): \[ E[aX+b] = aE[X]+b \qquad\text{and}\qquad \mathrm{Var}(aX+b) = a^2\,\mathrm{Var}(X). \]

Solution: identical to L4’s discrete proof, with sums replaced by integrals. \[ E[aX+b] = \int (ax+b)\, f_X(x)\,dx = a\int x\, f_X(x)\,dx + b\int f_X(x)\,dx = aE[X]+b. \] \[ \mathrm{Var}(aX+b) = E\big[(aX+b-aE[X]-b)^2\big] = E\big[a^2(X-E[X])^2\big] = a^2\,\mathrm{Var}(X). \]

Every fact about expectation and variance built on linearity in L4 carries over to continuous RVs with zero extra work.

The Exponential random variable

The continuous “cousin” of the Geometric — instead of counting trials until an event, it measures the waiting time until an event that can happen at any instant.

Definition: Exponential(\(\lambda\))

\[ f_X(x) = \begin{cases}\lambda e^{-\lambda x} & x \ge 0\\ 0 & x<0,\end{cases} \qquad \lambda>0. \]

\[ E[X] = \frac1\lambda, \qquad \mathrm{Var}(X) = \frac1{\lambda^2} \] (both derivable by integration by parts) — the same shape of answer as the Geometric(\(p\)) mean \(1/p\), with \(\lambda\) playing the role of a continuous rate in place of \(p\).

---

Example: customers arrive at a coffee shop at an average rate of \(\lambda=3\) per hour, memorylessly (the time since the last arrival tells you nothing about when the next one comes). The waiting time \(X\) until the next customer is Exponential(\(3\)), measured in hours, so \(E[X]=\tfrac13\) hour \(=20\) minutes.

def f_exp(x, lam):
    x = np.asarray(x, dtype=float)
    return np.where(x >= 0, lam * np.exp(-lam * x), 0.0)

xs = np.linspace(0, 3, 400)
fig, ax = plt.subplots(figsize=(6, 3.5))
for lam in [0.5, 1, 3]:
    ax.plot(xs, f_exp(xs, lam), label=f'$\\lambda={lam}$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f_X(x)$')
ax.legend()
plt.tight_layout()
plt.show()

lam = 3.0  # the coffee shop: 3 customers/hour
E_exp_num, _ = quad(lambda x: x * f_exp(x, lam), 0, 50)
E_exp2_num, _ = quad(lambda x: x ** 2 * f_exp(x, lam), 0, 50)
Var_exp_num = E_exp2_num - E_exp_num ** 2

rng = np.random.default_rng(151)
wait_times = rng.exponential(1 / lam, size=1_000_000)  # numpy's scale parameter is 1/lambda

print(f"lambda = {lam}")
print(f"E[X]:   formula 1/lambda   = {1/lam:.4f}   quad = {E_exp_num:.4f}   simulated = {wait_times.mean():.4f}")
print(f"Var(X): formula 1/lambda^2 = {1/lam**2:.4f}   quad = {Var_exp_num:.4f}   simulated = {wait_times.var():.4f}")

lambda = 3.0
E[X]:   formula 1/lambda   = 0.3333   quad = 0.3333   simulated = 0.3327
Var(X): formula 1/lambda^2 = 0.1111   quad = 0.1111   simulated = 0.1112

Cumulative Distribution Functions

Definition: CDF

For any random variable \(X\) — discrete or continuous — \[ F_X(x) = P(X \le x). \]

Unlike the PMF or PDF, the CDF is defined the same way regardless of type. And each of PMF, PDF, or CDF — whichever applies to \(X\) — unambiguously pins down the entire distribution: from any one of them, you can recover the others.

Properties of the CDF

Every CDF, discrete or continuous, satisfies:

• monotonically nondecreasing: \(x_1 \le x_2 \implies F_X(x_1)\le F_X(x_2)\),
• \(F_X(x) \to 0\) as \(x\to-\infty\), and \(F_X(x)\to 1\) as \(x\to\infty\),
• for discrete \(X\): \(F_X\) is a piecewise-constant step function — flat except for a jump of size \(p_X(x)\) at every \(x\) in the range of \(X\),
• for continuous \(X\): \(F_X\) is continuous everywhere (no jumps at all — consistent with \(P(X=x_0)=0\) for every \(x_0\)).

Relation between PDF and CDF, for \(X\) continuous: \[ F_X(x) = \int_{-\infty}^x f_X(t)\,dt, \qquad\text{equivalently}\qquad f_X(x) = F_X'(x) \] wherever \(F_X\) is differentiable — the Fundamental Theorem of Calculus!

---

from scipy.stats import binom

fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))

# discrete: Binomial(5, 0.4) CDF -- a step function
n_b, p_b = 5, 0.4
xs_d = np.linspace(-1, 6, 700)
cdf_d = binom.cdf(np.floor(xs_d), n_b, p_b)
axes[0].step(xs_d, cdf_d, where='post', color='black')
axes[0].set_title('discrete CDF\n(Binomial(5,0.4)) -- a step function')
axes[0].set_xlabel('$x$')
axes[0].set_ylabel('$F_X(x)$')

# continuous: Uniform(2,5) CDF -- continuous, and (in this special case) piecewise linear
xs_c = np.linspace(0, 7, 400)
cdf_c = np.clip((xs_c - a_u) / (b_u - a_u), 0, 1)
axes[1].plot(xs_c, cdf_c, color='black')
axes[1].set_title('continuous CDF\n(Uniform(2,5)) -- continuous, no jumps')
axes[1].set_xlabel('$x$')

plt.tight_layout()
plt.show()

Geometric and Exponential CDFs

\[ X\sim\text{Geometric}(p): \quad F_X(k) = 1-(1-p)^k, \qquad k=1,2,\dots \] \[ X\sim\text{Exponential}(\lambda): \quad F_X(x) = 1-e^{-\lambda x}, \qquad x\ge 0. \]

Same shape of formula — because the Exponential really is the continuous limit of the Geometric. Chop a unit of time into \(n\) tiny slices, let each slice independently “succeed” with probability \(p=\lambda/n\), and count time (not trials) until the first success. As \(n\to\infty\), the Geometric CDF — rescaled so trial number \(k\) becomes time \(k/n\) — converges to the Exponential(\(\lambda\)) CDF.

Aside

Both distributions are memoryless — the property we used to derive the Geometric’s mean carries over exactly to the Exponential, and we’ll put it to work again later this course.

---

lam = 2.0
t = np.linspace(0, 3, 400)
exp_cdf = 1 - np.exp(-lam * t)

fig, ax = plt.subplots(figsize=(6, 3.5))
ax.plot(t, exp_cdf, 'k--', lw=2, label=f'Exponential($\\lambda={lam}$) CDF')

for n in [2, 5, 20, 200]:
    p = lam / n
    k_vals = np.arange(0, int(3 * n) + 1)
    geo_cdf = 1 - (1 - p) ** k_vals
    ax.step(k_vals / n, geo_cdf, where='post', alpha=0.7, label=f'Geometric($\\lambda/n$), n={n}')

ax.set_xlabel('time $t = k/n$')
ax.set_ylabel('CDF')
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

# max pointwise CDF gap, shrinking as n grows
n_range = np.unique(np.geomspace(2, 5000, 25).astype(int))
max_gaps = []
for n in n_range:
    p = lam / n
    k_vals = np.arange(0, int(5 * n) + 1)
    geo_cdf = 1 - (1 - p) ** k_vals
    t_vals = k_vals / n
    exp_at_t = 1 - np.exp(-lam * t_vals)
    max_gaps.append(np.max(np.abs(geo_cdf - exp_at_t)))

fig2, ax2 = plt.subplots(figsize=(5, 3))
ax2.loglog(n_range, max_gaps, 'o-')
ax2.set_xlabel('$n$')
ax2.set_ylabel('max pointwise CDF gap')
plt.tight_layout()
plt.show()

Example: continuous RV in randomized trace estimation

Recall the randomized trace estimator: for \(A\in\mathbb R^{n\times n}\), draw \(z\in\mathbb R^n\) with i.i.d. mean-\(0\), variance-\(1\), pairwise independent entries and set \(\tau = z^\top A z\). The unbiasedness proof only ever used \[ E[z_i^2] = 1 \qquad\text{and}\qquad E[z_iz_j] = 0 \ (i\ne j), \] which followed from mean \(0\), variance \(1\), and pairwise independence — nothing about \(z_i\) being discrete. So a continuous \(z_i\) works exactly as well.

A continuous alternative

Instead of Rademacher entries \(z_i=\pm1\), draw \(z_i\) i.i.d. from a continuous distribution with mean \(0\) and variance \(1\) — for concreteness, the (bell-shaped) standard Gaussian, whose PDF we’ll meet formally in a later lecture. The same proof works exactly, just a different distribution!

---

rng = np.random.default_rng(151)
d = 50
A = rng.standard_normal((d, d)) + 5 * np.eye(d)
true_trace = np.trace(A)

num_samples = 20_000
Z_rad = rng.choice([-1.0, 1.0], size=(num_samples, d))  # discrete: Rademacher
Z_gauss = rng.standard_normal((num_samples, d))          # continuous: standard Gaussian

tau_rad = np.einsum('si,ij,sj->s', Z_rad, A, Z_rad)
tau_gauss = np.einsum('si,ij,sj->s', Z_gauss, A, Z_gauss)

print(f"true trace(A)                        = {true_trace:.4f}")
print(f"mean of tau, Rademacher z (discrete)  = {tau_rad.mean():.4f}")
print(f"mean of tau, Gaussian z (continuous)  = {tau_gauss.mean():.4f}")

true trace(A)                        = 250.0783
mean of tau, Rademacher z (discrete)  = 249.9342
mean of tau, Gaussian z (continuous)  = 249.5031