Chapter 34: Probability Distributions in Depth

What Is a Random Variable?

A random variable $X$ is a measurable function from a sample space $\Omega$ to the real line $\mathbb{R}$. This formal definition is less important than the intuition: $X$ is a quantity whose value is determined by a random experiment. Rolling a fair die gives $X \in \{1, 2, 3, 4, 5, 6\}$. Recording a person's height gives $X \in (0, \infty)$. The distinction between these two examples — countable vs. uncountable range — defines the two major classes.

Discrete Random Variables and PMFs

A discrete random variable takes values in a finite or countably infinite set. Its behavior is completely described by the probability mass function (PMF):

The PMF must satisfy $p_X(x) \geq 0$ for all $x$ and $\sum_{x} p_X(x) = 1$.

Example — Poisson distribution. The number of events in a fixed interval when events occur independently at rate $\lambda$ follows a Poisson distribution:

This models website hits per second, customer arrivals per hour, mutations per genome segment, and countless other count phenomena.

Continuous Random Variables and PDFs

A continuous random variable $X$ has a probability density function (PDF) $f_X(x)$ such that:

The PDF itself is not a probability — it can exceed 1. The probability lives in the area under the curve. Requirements: $f_X(x) \geq 0$ and $\int_{-\infty}^{\infty} f_X(x)\, dx = 1$. Crucially, $P(X = x) = 0$ for any single point; only intervals carry positive probability for continuous variables.

Example — Exponential distribution. The time between successive Poisson events has PDF:

This is the unique continuous memoryless distribution: $P(X > s + t \mid X > s) = P(X > t)$.

The Cumulative Distribution Function

The CDF unifies discrete and continuous cases:

For discrete $X$: $F_X(x) = \sum_{k \leq x} p_X(k)$ (a right-continuous step function). For continuous $X$: $F_X(x) = \int_{-\infty}^x f_X(t)\, dt$, and differentiating recovers the PDF. The CDF always satisfies: non-decreasing, right-continuous, $\lim_{x \to -\infty} F_X(x) = 0$, $\lim_{x \to \infty} F_X(x) = 1$.

Why the CDF matters in practice. Computing $P(a < X \leq b) = F_X(b) - F_X(a)$ requires only two CDF evaluations regardless of distribution complexity. The empirical CDF (ECDF) from data converges uniformly to the true CDF by the Glivenko-Cantelli theorem — the theoretical basis for the Kolmogorov-Smirnov test.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Compare PMF and CDF for Poisson(lambda=3)
lam = 3
k = np.arange(0, 15)
pmf = stats.poisson.pmf(k, mu=lam)
cdf = stats.poisson.cdf(k, mu=lam)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ax1.stem(k, pmf, basefmt=' ')
ax1.set_title('Poisson PMF (λ=3)')
ax1.set_xlabel('k'); ax1.set_ylabel('P(X=k)')

ax2.step(k, cdf, where='post')
ax2.set_title('Poisson CDF (λ=3)')
ax2.set_xlabel('k'); ax2.set_ylabel('P(X≤k)')
plt.tight_layout()
plt.show()

# Compute P(2 <= X <= 5) two ways
print('Direct sum:', sum(stats.poisson.pmf(range(2,6), mu=lam)))
print('CDF method:', stats.poisson.cdf(5, mu=lam) - stats.poisson.cdf(1, mu=lam))

Common Pitfalls

• Confusing PDF value with probability. A PDF value of 2.5 at a point is valid; it means probability density is high there, but no single point has nonzero probability.
• Using pmf when you need cdf or vice versa. Always check: <!--MATHBLOCK28--> is a CDF call; <!--MATHBLOCK29--> is a PMF call.
• Off-by-one in discrete CDFs. <!--MATHBLOCK30--> for integer-valued <!--MATHBLOCK31-->, not <!--MATHBLOCK32-->.

Expectation: The Center of Mass

The expectation (mean) of a random variable summarizes where its distribution is centered. For discrete $X$:

For continuous $X$:

Expectation is linear: $E[aX + bY] = aE[X] + bE[Y]$, regardless of dependence between $X$ and $Y$. This is the most-used property in statistics.

Law of the unconscious statistician (LOTUS). To compute $E[g(X)]$, you do not need the distribution of $g(X)$ — you integrate $g(x)$ against the original PDF:

Variance and Standard Deviation

Variance measures spread around the mean:

The second form is computationally convenient. Variance is not linear: $\text{Var}(aX + b) = a^2 \text{Var}(X)$. For independent $X, Y$: $\text{Var}(X + Y) = \text{Var}(X) + \text{Var}(Y)$.

Standard deviation $\sigma = \sqrt{\text{Var}(X)}$ restores the original units.

Higher Moments and Skewness/Kurtosis

The $k$-th raw moment is $E[X^k]$. The $k$-th central moment is $E[(X - \mu)^k]$.

Skewness measures asymmetry:

Positive skew means a long right tail (income, wait times). Negative skew means a long left tail.

Excess kurtosis measures tail heaviness relative to the normal:

The $-3$ makes the normal distribution have excess kurtosis of 0. Student's $t$ with low degrees of freedom has high kurtosis — its tails are heavier than normal, which matters for financial risk and robust statistics.

Moment-Generating Functions

The moment-generating function (MGF) $M_X(t) = E[e^{tX}]$ encodes all moments: $E[X^k] = M_X^{(k)}(0)$. More importantly, if two distributions have the same MGF, they are identical — MGFs uniquely determine distributions (when they exist). They also make it easy to find the distribution of sums of independent variables: if $X$ and $Y$ are independent, $M_{X+Y}(t) = M_X(t) \cdot M_Y(t)$.

Example. For $X \sim \text{Poisson}(\lambda)$: $M_X(t) = e^{\lambda(e^t - 1)}$. If $X_1 \sim \text{Poisson}(\lambda_1)$ and $X_2 \sim \text{Poisson}(\lambda_2)$ independently, then $M_{X_1+X_2}(t) = e^{(\lambda_1+\lambda_2)(e^t-1)}$, proving $X_1 + X_2 \sim \text{Poisson}(\lambda_1 + \lambda_2)$.

import numpy as np
from scipy import stats

# Verify moment formulas numerically for Gamma(alpha=3, beta=2)
alpha, beta = 3, 2  # shape, rate; scipy uses scale=1/rate
scale = 1 / beta
dist = stats.gamma(a=alpha, scale=scale)

# Analytic moments for Gamma(alpha, beta):
# E[X] = alpha/beta, Var(X) = alpha/beta^2
print(f'Analytic mean:  {alpha/beta:.4f}')
print(f'Analytic var:   {alpha/beta**2:.4f}')
print(f'scipy mean:     {dist.mean():.4f}')
print(f'scipy var:      {dist.var():.4f}')
print(f'scipy skew:     {dist.stats(moments="s"):.4f}')
print(f'Analytic skew:  {2/np.sqrt(alpha):.4f}')  # 2/sqrt(alpha)

# Monte Carlo verification
rng = np.random.default_rng(0)
samples = dist.rvs(size=1_000_000, random_state=rng)
print(f'\nMC mean:  {samples.mean():.4f}')
print(f'MC var:   {samples.var():.4f}')
print(f'MC skew:  {stats.skew(samples):.4f}')
print(f'MC kurt:  {stats.kurtosis(samples):.4f}')  # excess kurtosis
print(f'Analytic kurt: {6/alpha:.4f}')              # 6/alpha for Gamma

Standard Distributions: Key Moments at a Glance

• Binomial <!--MATHBLOCK35-->: mean <!--MATHBLOCK36-->, variance <!--MATHBLOCK37-->, skew <!--MATHBLOCK38-->.
• Poisson <!--MATHBLOCK39-->: mean <!--MATHBLOCK40-->, variance <!--MATHBLOCK41-->. Mean equals variance — a key diagnostic.
• Exponential <!--MATHBLOCK42-->: mean <!--MATHBLOCK43-->, variance <!--MATHBLOCK44-->, skew <!--MATHBLOCK45-->.
• Normal <!--MATHBLOCK46-->: mean <!--MATHBLOCK47-->, variance <!--MATHBLOCK48-->, skew <!--MATHBLOCK49-->, kurtosis <!--MATHBLOCK50-->.
• Uniform <!--MATHBLOCK51-->: mean <!--MATHBLOCK52-->, variance <!--MATHBLOCK53-->.

Pitfalls

Cauchy distribution. The ratio of two independent standard normals has a Cauchy distribution with PDF $f(x) = 1/(\pi(1+x^2))$. Its mean does not exist (the integral diverges). Averages of Cauchy samples do not converge — a dramatic illustration that moments can fail to exist and the CLT does not always apply.

Joint Distributions

When we study two or more random variables simultaneously, we need their joint distribution. For discrete $(X, Y)$ the joint PMF is $p_{X,Y}(x, y) = P(X = x, Y = y)$. For continuous $(X, Y)$ the joint PDF $f_{X,Y}(x, y)$ satisfies:

The joint distribution contains the complete probabilistic description of how $X$ and $Y$ relate.

Marginal Distributions

To recover the distribution of $X$ alone, marginalize over $Y$:

This is the probabilistic analog of integrating out a nuisance variable. In Bayesian inference, computing the marginal likelihood $p(\text{data}) = \int p(\text{data}|\theta) p(\theta)\, d\theta$ is exactly this operation.

Conditional Distributions

The conditional PDF of $Y$ given $X = x$ is:

This is Bayes' theorem in disguise. The conditional distribution describes how our uncertainty about $Y$ updates once we observe $X$.

Example — Bivariate Normal. If $(X, Y)$ is bivariate normal with means $\mu_X, \mu_Y$, variances $\sigma_X^2, \sigma_Y^2$, and correlation $\rho$, then:

This is the theoretical backbone of linear regression: the conditional mean is linear in $x$, and the conditional variance $\sigma_Y^2(1-\rho^2)$ shrinks when $|\rho|$ is large.

Independence

$X$ and $Y$ are independent if and only if their joint distribution factors:

Independence implies $E[XY] = E[X]E[Y]$, $\text{Var}(X+Y) = \text{Var}(X) + \text{Var}(Y)$, and knowing $X$ provides no information about $Y$.

Critical warning: The converse — zero correlation implies independence — is false in general. It holds only for jointly normal variables.

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

# Construct uncorrelated but DEPENDENT variables
X = rng.uniform(-1, 1, size=2000)
Y = X**2 + rng.normal(0, 0.05, size=2000)  # Y determined by X

print(f'Pearson correlation: {np.corrcoef(X, Y)[0,1]:.4f}')  # near 0
print(f'Spearman correlation: {stats.spearmanr(X, Y).statistic:.4f}')  # detects monotone

# Yet clearly dependent:
fig, ax = plt.subplots(figsize=(6, 4))
ax.scatter(X, Y, alpha=0.3, s=5)
ax.set_title('Zero Pearson correlation, strong dependence')
ax.set_xlabel('X'); ax.set_ylabel('Y = X² + ε')
plt.tight_layout()
plt.show()

Conditional Expectation and the Tower Property

The conditional expectation $E[Y|X]$ is itself a random variable — a function of $X$. The tower property (law of iterated expectations) states:

This is used constantly in probability calculations and in the proof that ordinary least-squares regression gives the best linear predictor.

The law of total variance decomposes uncertainty:

The first term is average within-group variance; the second is between-group variance. This decomposition underlies ANOVA, random-effects models, and variance component estimation.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Bivariate normal: verify conditional distribution
mu = np.array([1.0, 2.0])
rho = 0.7
cov = np.array([[1.0, rho], [rho, 1.0]])

samples = rng.multivariate_normal(mu, cov, size=200_000)
X_samples, Y_samples = samples[:, 0], samples[:, 1]

# Condition on X in narrow band around x0=0
x0 = 0.0
mask = np.abs(X_samples - x0) < 0.05
Y_given_X = Y_samples[mask]

cond_mean_empirical = Y_given_X.mean()
cond_var_empirical = Y_given_X.var()

# Theoretical values:
cond_mean_theory = mu[1] + rho * (x0 - mu[0])  # = 2 + 0.7*(0-1) = 1.3
cond_var_theory = 1 - rho**2                     # = 1 - 0.49 = 0.51

print(f'Conditional mean  — empirical: {cond_mean_empirical:.3f}, theory: {cond_mean_theory:.3f}')
print(f'Conditional var   — empirical: {cond_var_empirical:.3f}, theory: {cond_var_theory:.3f}')
print(f'Number of points in band: {mask.sum()}')

Covariance: Measuring Linear Co-movement

Covariance quantifies the degree to which two variables move together linearly:

Positive covariance: when $X$ is above its mean, $Y$ tends to be above its mean. Negative covariance: they move in opposite directions. Zero covariance: no linear association.

Covariance has an awkward scale-dependence: it changes when we rescale $X$ or $Y$. This motivates the normalized version.

Pearson Correlation

The Pearson correlation coefficient standardizes covariance:

$\rho = 1$ means perfect positive linear relationship; $\rho = -1$ means perfect negative linear relationship; $\rho = 0$ means no linear relationship (but possibly nonlinear dependence).

The Cauchy-Schwarz inequality $|\text{Cov}(X,Y)|^2 \leq \text{Var}(X)\text{Var}(Y)$ guarantees $|\rho| \leq 1$.

The Covariance Matrix

For a random vector $\mathbf{X} = (X_1, \ldots, X_p)^\top$, the covariance matrix is:

where $\Sigma_{ij} = \text{Cov}(X_i, X_j)$. The diagonal entries are variances. $\Sigma$ is always symmetric positive semidefinite (PSD): $\mathbf{v}^\top \Sigma \mathbf{v} \geq 0$ for all $\mathbf{v}$. Violating PSD is a sign of a bug — for instance, estimating a correlation matrix from insufficient data often produces one that is not PSD.

For an affine transformation $\mathbf{Y} = A\mathbf{X} + \mathbf{b}$:

This is the key formula behind principal component analysis (PCA), portfolio variance, and propagation of uncertainty.

Rank-Based Correlations

Pearson correlation only captures linear dependence. Two rank-based alternatives are robust to outliers and detect monotone (not just linear) relationships:

Spearman's $\rho_s$ is the Pearson correlation applied to the ranks of $X$ and $Y$.

Kendall's $\tau$ counts concordant minus discordant pairs:

import numpy as np
from scipy import stats

rng = np.random.default_rng(5)

# Example: exponentially distributed data, known analytic Pearson correlation
# Generate bivariate log-normal: X=exp(U), Y=exp(V), (U,V)~BVN(0,0,1,1,rho=0.8)
rho_latent = 0.8
cov = np.array([[1., rho_latent], [rho_latent, 1.]])
uv = rng.multivariate_normal([0., 0.], cov, size=100_000)
X = np.exp(uv[:, 0])
Y = np.exp(uv[:, 1])

# Analytic Pearson for log-normal: (exp(rho)-1)/(exp(1)-1)
import math
rho_pearson_theory = (math.exp(rho_latent) - 1) / (math.exp(1) - 1)

pearson_r, _ = stats.pearsonr(X, Y)
spearman_r, _ = stats.spearmanr(X, Y)
kendall_tau, _ = stats.kendalltau(X, Y)

print(f'Latent Gaussian correlation:   {rho_latent:.3f}')
print(f'Theoretical Pearson (log-norm):{rho_pearson_theory:.3f}')
print(f'Empirical Pearson:             {pearson_r:.3f}')
print(f'Spearman rho:                  {spearman_r:.3f}')
print(f"Kendall tau:                   {kendall_tau:.3f}")

Common Pitfalls

Outlier inflation. A single outlier can drive Pearson correlation from near 0 to 0.8. Always plot your data and consider Spearman for skewed variables.

Correlation is not causation. This caveat is statistical, not just philosophical: two variables correlated through a common cause have the same joint distribution as two causally linked variables. You cannot distinguish them from observational data alone without additional assumptions (causal DAG structure, instrumental variables, etc.).

Spurious correlations in time series. Two independent random walks will exhibit high correlation because both trend over time. Always detrend or difference before computing correlation of time series.

Restriction of range. If you compute correlation only within a subgroup (e.g., SAT scores for admitted students), you attenuate the observed correlation relative to the full population — known as Berkson's paradox or range restriction bias.

The Law of Large Numbers

Let $X_1, X_2, \ldots$ be i.i.d. with finite mean $\mu = E[X_i]$. Define the sample mean $\bar{X}_n = n^{-1}\sum_{i=1}^n X_i$. The Weak Law of Large Numbers (WLLN) states:

($\xrightarrow{p}$ means convergence in probability: $P(|\bar{X}_n - \mu| > \epsilon) \to 0$ for every $\epsilon > 0$.) The Strong LLN strengthens this to almost-sure convergence.

The LLN is the theoretical basis for Monte Carlo integration: to estimate $E[g(X)]$, draw many samples $x_i \sim p$ and compute $(1/n)\sum g(x_i)$. With enough samples, this is as accurate as desired.

The Central Limit Theorem

The LLN tells us where $\bar{X}_n$ lands; the Central Limit Theorem (CLT) tells us the shape of its distribution. If $X_i$ are i.i.d. with mean $\mu$ and finite variance $\sigma^2$, then:

($\xrightarrow{d}$ means convergence in distribution.) Equivalently, $\bar{X}_n \approx N(\mu, \sigma^2/n)$ for large $n$.

Why this is profound. The CLT holds for any distribution with finite variance — Poisson, Bernoulli, exponential, uniform, whatever. The normal distribution emerges as a universal limit, which explains its ubiquity. The convergence rate is $O(1/\sqrt{n})$, meaning quadrupling sample size halves the spread.

Berry-Esseen theorem (quantitative CLT). The maximum deviation of the true CDF of $\bar{X}_n$ from the normal approximation is bounded:

where $C < 0.5$. Convergence is slower when the distribution is skewed or has heavy tails.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

rng = np.random.default_rng(42)
B = 5000  # bootstrap repetitions

fig, axes = plt.subplots(2, 3, figsize=(14, 8))

for row, dist_name, dist_fn, mu, sigma in [
    (0, 'Bernoulli(0.1)', lambda n: rng.binomial(1, 0.1, n), 0.1, np.sqrt(0.09)),
    (1, 'Exponential(1)', lambda n: rng.exponential(1, n), 1.0, 1.0),
]:
    for col, n in enumerate([5, 30, 300]):
        means = np.array([dist_fn(n).mean() for _ in range(B)])
        standardized = (means - mu) / (sigma / np.sqrt(n))

        x = np.linspace(-4, 4, 300)
        axes[row][col].hist(standardized, bins=60, density=True,
                            alpha=0.6, label='Simulated')
        axes[row][col].plot(x, stats.norm.pdf(x), 'r-', lw=2, label='N(0,1)')
        axes[row][col].set_title(f'{dist_name}\nn={n}')
        if col == 0:
            axes[row][col].set_ylabel('Density')
        if row == 1:
            axes[row][col].set_xlabel('Standardized mean')
        axes[row][col].legend(fontsize=7)

plt.suptitle('CLT Convergence: standardized sample means', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

When the CLT Fails

• Infinite variance. Student's <!--MATHBLOCK24--> (Cauchy) — variance undefined, CLT does not apply.
• Heavy tails. For <!--MATHBLOCK25--> with <!--MATHBLOCK26-->, variance exists but kurtosis is large; CLT convergence is slow.
• Dependence. Correlated data violates i.i.d. For time-series data, CLT generalizations require conditions like mixing.
• Small <!--MATHBLOCK27--> with skewed data. For Bernoulli(<!--MATHBLOCK28-->), <!--MATHBLOCK29--> is far too small for the normal approximation.

Transformations of Random Variables

Change of Variables

If $Y = g(X)$ where $g$ is monotone and differentiable, the PDF of $Y$ is:

The Jacobian term $|dg^{-1}/dy|$ accounts for how $g$ stretches or compresses probability mass.

Example. If $X \sim N(0,1)$ and $Y = e^X$, then $Y$ is log-normal with $f_Y(y) = (1/y)\phi(\ln y)$ for $y > 0$, where $\phi$ is the standard normal PDF.

The Delta Method

For large samples, if $\sqrt{n}(\hat{\theta} - \theta) \xrightarrow{d} N(0, \sigma^2)$ and $g$ is differentiable at $\theta$, then:

This is the workhorse for approximate variance propagation: if you know the variance of $\hat{\theta}$, the delta method gives you the approximate variance of any smooth function of $\hat{\theta}$.

Example. Proportion $\hat{p}$ estimates $p$ with variance $p(1-p)/n$. For the log-odds $g(p) = \ln(p/(1-p))$, $g'(p) = 1/(p(1-p))$, so $\text{Var}(g(\hat{p})) \approx 1/(np(1-p))$. This is the theoretical basis for the standard error of a logistic regression coefficient.

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
n = 500
p_true = 0.3
B = 50_000

p_hat = rng.binomial(n, p_true, size=B) / n  # sample proportions
logodds_hat = np.log(p_hat / (1 - p_hat))

# Delta method variance for log-odds
delta_var = 1 / (n * p_true * (1 - p_true))
print('Delta method SE of log-odds:', np.sqrt(delta_var))
print('Simulated SE of log-odds:   ', logodds_hat.std())

true_logodds = np.log(p_true / (1 - p_true))
standardized = (logodds_hat - true_logodds) / np.sqrt(delta_var)
print(f'\nKS test vs N(0,1): stat={stats.kstest(standardized, "norm").statistic:.4f}')