Probability Theory

Every probability model starts with three ingredients:

1. Sample space $\Omega$: the set of all possible outcomes of an experiment. For a fair coin, $\Omega = {H, T}$. For two dice, $\Omega = {(i,j) : 1 \le i,j \le 6}$, a set of 36 ordered pairs. For a random real number, $\Omega = \mathbb{R}$.

2. Event: any subset $A \subseteq \Omega$. The event "sum of two dice equals 7" is ${(1,6),(2,5),(3,4),(4,3),(5,2),(6,1)}$.

3. Probability measure $P$: a function from events to $[0,1]$ satisfying Kolmogorov's three axioms:

   • Non-negativity: $P(A) \ge 0$ for all events $A$.
   • Normalization: $P(\Omega) = 1$.
   • Countable additivity: if $A_1, A_2, \ldots$ are pairwise disjoint events, then $P(\bigcup_{i=1}^\infty A_i) = \sum_{i=1}^\infty P(A_i)$.

From these three axioms, everything else follows. For example: $P(\emptyset) = 0$ (the impossible event has probability zero). $P(A^c) = 1 - P(A)$ (complement rule). $P(A \cup B) = P(A) + P(B) - P(A \cap B)$ (inclusion-exclusion). These are all theorems, not additional axioms.

The countable additivity axiom (also called $\sigma$-additivity) is the subtle one — it extends finite additivity to infinite disjoint unions. This is what makes the theory work for continuous distributions, where individual outcomes have probability zero but ranges have positive probability.

Analogy: Think of a probability measure as a way of distributing a total "mass" of 1 across the sample space. Events are regions of the space, and their probability is their share of the total mass.

import itertools

# Sample space: two dice
omega = list(itertools.product(range(1,7), range(1,7)))

def P(event):
    """Uniform probability: |event| / |omega|"""
    return len(event) / len(omega)

# Event: sum == 7
sum7 = [(i,j) for i,j in omega if i+j == 7]
print(P(sum7))   # 6/36 = 0.1667

# Complement rule: P(not sum7) = 1 - P(sum7)
print(1 - P(sum7))   # 0.8333

# Inclusion-exclusion: P(A or B)
A = [(i,j) for i,j in omega if i==6]    # first die is 6
B = [(i,j) for i,j in omega if j==6]    # second die is 6
AorB = list(set(A) | set(B))
print(P(AorB))           # = P(A)+P(B)-P(A∩B) = 11/36

Conditional probability is how we update beliefs when we learn new information. Given that event $B$ has occurred (with $P(B) > 0$), the probability of $A$ is: $$P(A \mid B) = \frac{P(A \cap B)}{P(B)}$$

Intuition: we restrict the sample space to $B$ and renormalize.

Bayes' Theorem follows immediately from symmetry of the definition: $$P(B \mid A) = \frac{P(A \mid B), P(B)}{P(A)}$$ This is arguably the most important equation in all of data science. It tells you how to update a prior belief $P(B)$ in light of evidence $A$, producing a posterior $P(B \mid A)$.

Example (medical test): A disease affects 1% of the population ($P(D) = 0.01$). A test has 95% sensitivity ($P(\text{pos} \mid D) = 0.95$) and 90% specificity ($P(\text{neg} \mid \text{no D}) = 0.90$, so false-positive rate $= 0.10$). What is $P(D \mid \text{pos})$? $$P(D \mid +) = \frac{0.95 \times 0.01}{0.95 \times 0.01 + 0.10 \times 0.99} \approx 0.0875$$ Only about 9%! Most positive tests are false positives because the disease is rare. This base-rate neglect surprises almost everyone.

Independence: $A$ and $B$ are independent if $P(A \cap B) = P(A),P(B)$, equivalently $P(A \mid B) = P(A)$ — learning $B$ tells us nothing about $A$. Events can be pairwise independent yet not mutually independent (a subtlety that bites students and practitioners alike).

The chain rule: For any events $A_1, \ldots, A_n$: $$P(A_1 \cap A_2 \cap \cdots \cap A_n) = P(A_1),P(A_2 \mid A_1),P(A_3 \mid A_1,A_2)\cdots$$ This is the foundation of Bayesian networks and probabilistic graphical models.

# Bayes' theorem: medical test example
P_D   = 0.01           # prior: disease prevalence
P_pos_given_D      = 0.95   # sensitivity
P_pos_given_noD    = 0.10   # 1 - specificity

# Law of total probability for P(positive)
P_pos = P_pos_given_D * P_D + P_pos_given_noD * (1 - P_D)

# Posterior
P_D_given_pos = (P_pos_given_D * P_D) / P_pos
print(f"P(disease | positive test) = {P_D_given_pos:.4f}")
# 0.0875 — only 8.75%! Most positives are false alarms.

A random variable $X$ is a function $X: \Omega \to \mathbb{R}$ that assigns a real number to each outcome. It is not "variable" in the algebra sense — it is a deterministic function; the randomness comes from the random draw of $\omega \in \Omega$.

The distribution of $X$ is the complete description of the probabilities of all values $X$ can take:

• For a discrete RV: the probability mass function (PMF) $p(x) = P(X = x)$, which must satisfy $p(x) \ge 0$ and $\sum_x p(x) = 1$.
• For a continuous RV: the probability density function (PDF) $f(x) \ge 0$ with $\int_{-\infty}^{\infty} f(x),dx = 1$. Here $P(a \le X \le b) = \int_a^b f(x),dx$. Note: $P(X = x) = 0$ for any single point $x$.

The cumulative distribution function (CDF) $F(x) = P(X \le x)$ works for both discrete and continuous RVs. It is always non-decreasing, right-continuous, with $F(-\infty) = 0$ and $F(\infty) = 1$.

Key insight: Two different random variables on different sample spaces can have the same distribution. The distribution (PMF or PDF) is what we usually care about in practice — not the underlying sample space.

Analogy: The sample space $\Omega$ is like the raw machine code of an experiment. The random variable is a compiler that translates outcomes into numbers you can do arithmetic with. The distribution is the object-level description everyone shares.

import numpy as np
import matplotlib.pyplot as plt

# Discrete RV: roll a fair die
X_values = [1, 2, 3, 4, 5, 6]
pmf = [1/6] * 6                # PMF
cdf = np.cumsum(pmf)           # CDF

print("PMF:", pmf)
print("CDF:", cdf)

# Continuous RV: standard normal
from scipy.stats import norm
x = np.linspace(-4, 4, 200)
print(f"P(-1 <= Z <= 1) = {norm.cdf(1) - norm.cdf(-1):.4f}")  # 0.6827
print(f"P(-2 <= Z <= 2) = {norm.cdf(2) - norm.cdf(-2):.4f}")  # 0.9545

Bernoulli($p$): The simplest non-trivial distribution. A single binary trial with success probability $p \in [0,1]$. The RV $X$ takes value 1 (success) with probability $p$ and 0 (failure) with probability $1-p$. $$P(X = k) = p^k (1-p)^{1-k}, \quad k \in {0,1}$$

• $E[X] = p$, $\text{Var}(X) = p(1-p)$.
• Interpretation: a coin flip with bias $p$, a single click or no-click, a bit.

Binomial($n$, $p$): The sum of $n$ independent Bernoulli($p$) trials. $X$ counts the number of successes in $n$ trials. $$P(X = k) = \binom{n}{k} p^k (1-p)^{n-k}, \quad k = 0, 1, \ldots, n$$ The binomial coefficient $\binom{n}{k}$ counts the number of ways to place $k$ successes among $n$ positions; $p^k(1-p)^{n-k}$ is the probability of any specific arrangement.

• $E[X] = np$, $\text{Var}(X) = np(1-p)$.
• Example: 20 students each independently pass an exam with probability 0.7. The number who pass is Binomial(20, 0.7) with mean 14.

Normal approximation: For large $n$, Binomial($n$,$p$) $\approx$ Normal($np$, $np(1-p)$). This is the Central Limit Theorem in action — why proportions are approximately Gaussian in large samples.

from scipy.stats import binom
import numpy as np

# Binomial(20, 0.7): number passing an exam
n, p = 20, 0.7

# PMF: P(X = k) for each k
for k in [10, 14, 15, 20]:
    print(f"P(X={k}) = {binom.pmf(k, n, p):.4f}")

# Expected value and variance
print(f"Mean = {binom.mean(n, p)}")   # 14.0
print(f"Var  = {binom.var(n, p)}")    # 4.2

# P(at least 15 pass)
print(f"P(X >= 15) = {1 - binom.cdf(14, n, p):.4f}")

The Gaussian distribution $\mathcal{N}(\mu, \sigma^2)$ is the most important distribution in all of probability and statistics. Its PDF is: $$f(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$ where $\mu$ is the mean and $\sigma^2$ is the variance ($\sigma$ is the standard deviation). The "bell curve" shape comes from the squared exponent: probability falls off exponentially fast as you move away from the mean.

Standard normal $Z = \mathcal{N}(0, 1)$: any Gaussian can be standardized as $Z = (X - \mu)/\sigma$. The 68-95-99.7 rule: $P(\mu - \sigma \le X \le \mu + \sigma) \approx 0.6827$; within $2\sigma$, about 95.45%; within $3\sigma$, about 99.73%.

Why Gaussian is everywhere — the Central Limit Theorem: the sum (or average) of many independent, identically distributed random variables with finite mean and variance converges in distribution to a Gaussian, regardless of the original distribution. This is why measurement errors (sums of many small independent perturbations) are Gaussian, why test scores aggregate this way, why returns on diversified portfolios trend Gaussian.

Properties:

• Closed under linear transformations: if $X \sim \mathcal{N}(\mu, \sigma^2)$, then $aX + b \sim \mathcal{N}(a\mu + b, a^2\sigma^2)$.
• Sum of independent Gaussians is Gaussian: $\mathcal{N}(\mu_1, \sigma_1^2) + \mathcal{N}(\mu_2, \sigma_2^2) = \mathcal{N}(\mu_1+\mu_2, \sigma_1^2+\sigma_2^2)$.
• Maximum entropy: among all distributions with given mean and variance, the Gaussian has the highest entropy — it makes the weakest possible assumptions beyond those two numbers.

Caution: Real-world data often has heavier tails than Gaussian (financial returns, insurance claims). The Gaussian is a model, not a law of nature.

from scipy.stats import norm
import numpy as np

mu, sigma = 170, 10   # heights in cm

# 68-95-99.7 rule
print(norm.cdf(mu+sigma, mu, sigma) - norm.cdf(mu-sigma, mu, sigma))   # ~0.6827
print(norm.cdf(mu+2*sigma, mu, sigma) - norm.cdf(mu-2*sigma, mu, sigma)) # ~0.9545

# Standardize: P(X > 185)
z = (185 - mu) / sigma   # z = 1.5
print(1 - norm.cdf(z))   # P(Z > 1.5) = 0.0668

# Inverse CDF (quantile): height at 90th percentile
print(norm.ppf(0.90, mu, sigma))   # ~182.8 cm

Poisson($\lambda$): Counts the number of events in a fixed interval when events arrive at a constant average rate $\lambda$ and are independent of each other. Classic examples: phone calls arriving at a call center per hour, mutations per genome per generation, clicks on a website per minute, radioactive decays per second. $$P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \ldots$$

• $E[X] = \lambda$, $\text{Var}(X) = \lambda$ (mean equals variance — a hallmark of Poisson).
• Poisson limit theorem: Binomial($n$,$p$) $\to$ Poisson($\lambda = np$) as $n \to \infty$, $p \to 0$ with $np$ fixed. So the Poisson models "rare events in large populations".

Exponential($\lambda$): The waiting time between Poisson events. If events arrive at rate $\lambda$, the time between consecutive events is Exponential($\lambda$): $$f(t) = \lambda e^{-\lambda t}, \quad t \ge 0$$

• $E[T] = 1/\lambda$, $\text{Var}(T) = 1/\lambda^2$.
• Memoryless property: $P(T > s + t \mid T > s) = P(T > t)$. The exponential distribution has no memory of how long you've already waited. The only continuous memoryless distribution. This makes it the building block of queueing theory and reliability analysis.

Geometric distribution (discrete analog): the number of Bernoulli trials until the first success; the only discrete memoryless distribution.

from scipy.stats import poisson, expon
import numpy as np

# Poisson: website gets avg 5 requests/second
lam = 5

print(f"P(0 requests)  = {poisson.pmf(0, lam):.4f}")  # 0.0067
print(f"P(5 requests)  = {poisson.pmf(5, lam):.4f}")  # 0.1755
print(f"P(>10 requests)= {1 - poisson.cdf(10, lam):.4f}")  # 0.0137

# Exponential: time (sec) between requests
scale = 1/lam  # mean = 1/lambda
print(f"Mean wait = {expon.mean(scale=scale):.2f} sec")  # 0.20
# Memoryless: P(T > 0.5 | T > 0.3) = P(T > 0.2)
print(f"P(T > 0.2) = {1 - expon.cdf(0.2, scale=scale):.4f}")

The expected value (expectation, mean) $E[X]$ is the probability-weighted average of all values $X$ can take: $$E[X] = \sum_x x \cdot P(X = x) \quad \text{(discrete)}$$ $$E[X] = \int_{-\infty}^\infty x \cdot f(x),dx \quad \text{(continuous)}$$

Think of $E[X]$ as the center of mass of the probability distribution.

Linearity of expectation: $E[aX + bY] = aE[X] + bE[Y]$, for any random variables $X$ and $Y$, regardless of dependence. This is enormously powerful — it lets you compute expectations of sums by summing expectations, even when the summands are correlated.

Classic application: Expected number of fixed points of a random permutation of $n$ elements. Define $X_i = 1$ if element $i$ is a fixed point. $E[X_i] = 1/n$. By linearity, $E[\sum_i X_i] = \sum_i 1/n = 1$. Elegant!

Variance: $\text{Var}(X) = E[(X - E[X])^2] = E[X^2] - (E[X])^2$. It measures spread around the mean. The standard deviation $\sigma = \sqrt{\text{Var}(X)}$ has the same units as $X$.

Key variance rules: $\text{Var}(aX + b) = a^2 \text{Var}(X)$. For independent $X$, $Y$: $\text{Var}(X+Y) = \text{Var}(X) + \text{Var}(Y)$. (This additivity for independent sums is what drives the $\sqrt{n}$ scaling in the CLT.)

Covariance and correlation: $\text{Cov}(X,Y) = E[(X - E[X])(Y - E[Y])] = E[XY] - E[X]E[Y]$. The Pearson correlation $\rho(X,Y) = \text{Cov}(X,Y) / (\sigma_X \sigma_Y) \in [-1, 1]$ is the normalized version. Independent $\Rightarrow$ $\text{Cov} = 0$, but the converse is false in general.

import numpy as np

# Linearity of expectation: fixed-point example
# Random permutation of n=5 elements
n = 5
trial_count = 100_000
fixed_point_counts = []
for _ in range(trial_count):
    perm = np.random.permutation(n)
    fixed = np.sum(perm == np.arange(n))
    fixed_point_counts.append(fixed)
print(f"Simulated E[#fixed points] = {np.mean(fixed_point_counts):.4f}")
# Should be close to 1.0 regardless of n

# Variance additivity (independent RVs)
X = np.random.normal(3, 2, size=100_000)   # N(3,4)
Y = np.random.normal(1, 3, size=100_000)   # N(1,9)
Z = X + Y
print(f"Var(Z) ≈ {np.var(Z):.2f}, expected 13")

When we have two or more random variables, we need to describe their joint behavior.

Joint PMF/PDF: For discrete $(X, Y)$, $p(x,y) = P(X=x, Y=y)$. For continuous, $f(x,y)$ with $\iint f(x,y),dx,dy = 1$.

Marginal distribution: "Summing out" (or integrating out) one variable to get the distribution of the other: $$p_X(x) = \sum_y p(x,y) \quad \text{or} \quad f_X(x) = \int f(x,y),dy$$ The marginal distribution of $X$ describes $X$ alone, ignoring $Y$.

Conditional distribution: Given $X = x$, the conditional distribution of $Y$ is: $$p_{Y|X}(y \mid x) = \frac{p(x,y)}{p_X(x)}$$ This is just Bayes' theorem applied to densities.

Independence of RVs: $X$ and $Y$ are independent iff $p(x,y) = p_X(x) \cdot p_Y(y)$ for all $(x,y)$. The joint distribution factors into the product of the marginals.

Law of total expectation: $E[X] = E[E[X \mid Y]]$. You can compute the expectation of $X$ by first conditioning on $Y$ and then taking the expectation over $Y$. This is a workhorse of computation — it lets you break a hard expectation into easier conditional expectations.

Multivariate Gaussian: The $d$-dimensional Gaussian $\mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ has mean vector $\boldsymbol{\mu} \in \mathbb{R}^d$ and covariance matrix $\boldsymbol{\Sigma} \in \mathbb{R}^{d \times d}$ (positive semi-definite). All marginals and conditionals of a multivariate Gaussian are themselves Gaussian — a fact that makes Gaussian processes and Kalman filters tractable.

import numpy as np

# Joint distribution of (X, Y): bivariate discrete
# X = first die, Y = second die (independent)
joint = np.full((6, 6), 1/36)   # each (i,j) equally likely

# Marginal of X: sum over Y
marginal_X = joint.sum(axis=1)   # axis=1 sums over Y
print("Marginal P(X=k):", marginal_X)   # all 1/6

# Conditional: P(Y=y | X=1)
x_val = 0   # X=1 (0-indexed)
cond_Y_given_X1 = joint[x_val] / marginal_X[x_val]
print("P(Y|X=1):", cond_Y_given_X1)    # all 1/6 (independent)

# Law of total expectation: E[X] = E[E[X | Y]]
# Since independent, E[X | Y=y] = E[X] = 3.5 for all y
print("E[X] =", np.sum((np.arange(1,7)) * marginal_X))  # 3.5

The Law of Large Numbers (LLN) formalizes the intuition that averages stabilize as samples grow. Let $X_1, X_2, \ldots$ be i.i.d. (independent, identically distributed) with mean $\mu = E[X_i]$. Define the sample mean $\bar{X}n = \frac{1}{n}\sum{i=1}^n X_i$. Then:

Weak LLN: $\bar{X}_n \xrightarrow{P} \mu$ as $n \to \infty$. That is, for any $\varepsilon > 0$, $P(|\bar{X}_n - \mu| > \varepsilon) \to 0$. (Convergence in probability.)

Strong LLN: $\bar{X}_n \to \mu$ almost surely — the sample path of averages converges with probability 1.

Proof sketch of Weak LLN via Chebyshev's inequality: Chebyshev says $P(|Y - E[Y]| \ge k) \le \text{Var}(Y)/k^2$. Apply with $Y = \bar{X}_n$: $E[\bar{X}_n] = \mu$ and $\text{Var}(\bar{X}_n) = \sigma^2/n$ (by variance additivity for independent sums). So $P(|\bar{X}_n - \mu| \ge \varepsilon) \le \sigma^2/(n\varepsilon^2) \to 0$.

What it means: The sample mean is a consistent estimator of the true mean. With enough data, our empirical estimate will be arbitrarily close to the truth. This is the mathematical foundation of statistics, polling, Monte Carlo simulation (estimate integrals by averaging random samples), and stochastic gradient descent.

What it does NOT say: It says nothing about the rate of convergence (that's the CLT) or about any individual sequence of outcomes (gambler's fallacy: past outcomes don't change future probabilities).

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

# Fair die: true mean = 3.5
rolls = rng.integers(1, 7, size=10_000)
cumulative_means = np.cumsum(rolls) / np.arange(1, 10_001)

# The cumulative mean converges to 3.5
print(f"After   10 rolls: mean = {cumulative_means[9]:.3f}")
print(f"After  100 rolls: mean = {cumulative_means[99]:.3f}")
print(f"After 1000 rolls: mean = {cumulative_means[999]:.3f}")
print(f"After 10000 rolls: mean = {cumulative_means[-1]:.3f}")
# Should see convergence toward 3.5

The Central Limit Theorem (CLT) is perhaps the most remarkable theorem in probability. It says that the standardized sum of many i.i.d. random variables converges to the standard normal, regardless of the original distribution (as long as mean $\mu$ and variance $\sigma^2 < \infty$ exist).

Formally: let $X_1, \ldots, X_n$ be i.i.d. with mean $\mu$ and variance $\sigma^2$. Then: $$\sqrt{n},\frac{\bar{X}_n - \mu}{\sigma} \xrightarrow{d} \mathcal{N}(0,1)$$

Equivalently, $\sum_{i=1}^n X_i \approx \mathcal{N}(n\mu,, n\sigma^2)$ for large $n$.

Interpretation: The LLN says the mean converges; the CLT describes how fast and what the fluctuations look like — they are Gaussian with standard deviation $\sigma/\sqrt{n}$. Doubling precision requires $4\times$ more data.

Why it works intuitively: The characteristic function (Fourier transform of the PDF) of a sum of i.i.d. variables is a product of characteristic functions. Near the origin, every characteristic function looks like $e^{i\mu t - \sigma^2 t^2/2 + \ldots}$. After $n$ copies and standardization, higher-order terms vanish and you get exactly the Gaussian characteristic function $e^{-t^2/2}$.

In practice:

• A/B testing: the difference in sample means is approximately Gaussian, justifying z-tests and t-tests.
• Neural network training: weight updates from minibatch SGD are noisy but approximately Gaussian.
• Monte Carlo: the error in estimating an integral by $n$ samples is $O(1/\sqrt{n})$ with Gaussian distribution.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Start from a HIGHLY non-Gaussian distribution: Exponential(1)
# Mean=1, Var=1. Sum of n such RVs.
n_per_sample = 30
num_samples = 10_000

samples = rng.exponential(scale=1.0, size=(num_samples, n_per_sample))
sample_means = samples.mean(axis=1)

# Standardize
std_means = (sample_means - 1.0) / (1.0 / np.sqrt(n_per_sample))

# Test for normality: KS test against N(0,1)
stat, p = stats.kstest(std_means, 'norm')
print(f"KS p-value = {p:.4f}")   # Large p -> not rejected -> looks Gaussian!
print(f"Skewness = {stats.skew(std_means):.4f}")  # Should be near 0

Probability inequalities give tail bounds — upper bounds on how far a random variable can stray from its mean. They are essential for algorithm analysis and proving convergence results.

Markov's Inequality: For any non-negative RV $X$ and $a > 0$: $$P(X \ge a) \le \frac{E[X]}{a}$$ Proof: $E[X] = \int_0^\infty x f(x),dx \ge \int_a^\infty x f(x),dx \ge a \int_a^\infty f(x),dx = a, P(X \ge a)$. Divide by $a$.

This is weak (uses only the mean) but universal.

Chebyshev's Inequality: For any RV $X$ with mean $\mu$ and variance $\sigma^2$, and $k > 0$: $$P(|X - \mu| \ge k\sigma) \le \frac{1}{k^2}$$ Proof: Apply Markov to $(X - \mu)^2$ with $a = k^2\sigma^2$.

So at least $1 - 1/k^2$ of the probability lies within $k$ standard deviations of the mean, for any distribution. For $k=2$: at least 75% within $2\sigma$; for $k=3$: at least 89%. The Gaussian beats these bounds substantially (95.45% and 99.73%), but Chebyshev holds universally.

Chernoff bounds (beyond this module): For sums of independent 0/1 random variables, exponentially tighter bounds are available. These are the tools for analyzing hash tables, randomized algorithms, and boosting in machine learning.

# Chebyshev illustration: compare bound vs. actual for Uniform[0,1]
import numpy as np
from scipy.stats import uniform

# Uniform[0,1]: mean=0.5, var=1/12, std=1/sqrt(12)
mu = 0.5
sigma = (1/12)**0.5

for k in [1, 2, 3]:
    cheb_bound = 1 / k**2
    # Actual probability for Uniform: P(|X-0.5| >= k*sigma)
    lo = mu - k*sigma
    hi = mu + k*sigma
    # P outside [lo, hi]
    outside = max(0, lo) + max(0, 1-hi)   # for Uniform[0,1]
    print(f"k={k}: Chebyshev <= {cheb_bound:.3f}, Actual = {outside:.4f}")