Discrete Probability

A probability space is a triple (Ω, F, P) where:

• Ω (omega) is the sample space — the set of all possible outcomes of an experiment. It can be finite ({H,T} for a coin flip), countably infinite (ℕ for the number of coin flips until the first head), or uncountable (but we stay discrete here).
• F is a collection of events — subsets of Ω we care about. For finite Ω, F is just the power set 2^Ω.
• P is a probability measure: a function from events to [0,1] satisfying Kolmogorov's three axioms:
  1. Non-negativity: P(A) ≥ 0 for all events A.
  2. Normalization: P(Ω) = 1.
  3. Countable additivity: If A₁, A₂, … are pairwise disjoint events, then P(⋃ᵢ Aᵢ) = Σᵢ P(Aᵢ).

These three axioms are the entire foundation. Everything else — conditional probability, expectation, variance — is a theorem derived from them.

Uniform probability space: When Ω is finite and every outcome is equally likely, P(A) = |A| / |Ω|. This connects probability directly to counting (combinatorics).

Analogy: Think of Ω as a dart board where you throw darts uniformly at random. An event A is a painted region on the board. Its probability is the fraction of the board's area it covers.

Why axiom 3 matters: It lets us compute the probability of infinite unions by summing, and it implies all the useful rules below (complement rule, inclusion-exclusion, etc.) as theorems, not assumptions.

import random

# Sample space for rolling two dice
omega = [(i, j) for i in range(1, 7) for j in range(1, 7)]

# Event: sum equals 7
event_sum7 = [outcome for outcome in omega if sum(outcome) == 7]

# Probability (uniform space)
prob = len(event_sum7) / len(omega)
print(f"P(sum=7) = {len(event_sum7)}/{len(omega)} = {prob:.4f}")
# P(sum=7) = 6/36 = 0.1667

# Verify: empirical estimate via simulation
hits = sum(1 for _ in range(100_000) if sum(random.choices(range(1,7), k=2)) == 7)
print(f"Empirical estimate: {hits/100_000:.4f}")

Events are just sets, so all set operations have direct probabilistic interpretations:

Complement rule (from axioms 2 and 3): P(Aᶜ) = 1 − P(A). This is one of the most useful identities: it is often much easier to compute P("at least one X") as 1 − P("no X").

Inclusion-Exclusion (the general form): P(A₁ ∪ A₂ ∪ … ∪ Aₙ) = Σ P(Aᵢ) − Σ P(Aᵢ ∩ Aⱼ) + Σ P(Aᵢ ∩ Aⱼ ∩ Aₖ) − …

The alternating sum with subsets of increasing size. This is a theorem derived from axiom 3 by induction. It is powerful but expensive for large n — which is why the union bound (covered later) is often preferred.

DeMorgan in probability: Since (A ∪ B)ᶜ = Aᶜ ∩ Bᶜ and (A ∩ B)ᶜ = Aᶜ ∪ Bᶜ, you can always move complementation through unions/intersections, flipping them. This lets you convert hard union problems to intersection problems and vice versa.

# Verify inclusion-exclusion on rolling two dice
omega = [(i, j) for i in range(1, 7) for j in range(1, 7)]

A = {o for o in omega if o[0] == 6}   # first die is 6
B = {o for o in omega if o[1] == 6}   # second die is 6

p_A = len(A) / len(omega)             # 6/36
p_B = len(B) / len(omega)             # 6/36
p_AB = len(A & B) / len(omega)        # 1/36
p_AuB_incl_excl = p_A + p_B - p_AB   # 11/36

# Direct count
p_AuB_direct = len(A | B) / len(omega)

assert abs(p_AuB_incl_excl - p_AuB_direct) < 1e-12
print(f"P(A∪B) = {p_AuB_direct:.6f}  ({len(A|B)}/36)")

Conditional probability formalizes the act of updating beliefs given new information. The conditional probability of A given B is:

P(A | B) = P(A ∩ B) / P(B), provided P(B) > 0.

Intuition: When we learn that B occurred, we restrict the sample space to B and renormalize. P(A | B) asks: among all outcomes where B happened, what fraction also satisfy A?

The multiplication rule follows directly: P(A ∩ B) = P(A | B) · P(B) = P(B | A) · P(A). This chain extends to n events: P(A₁ ∩ A₂ ∩ … ∩ Aₙ) = P(A₁) · P(A₂|A₁) · P(A₃|A₁∩A₂) · … · P(Aₙ|A₁∩…∩Aₙ₋₁)

Law of Total Probability: If B₁, B₂, …, Bₙ is a partition of Ω (pairwise disjoint, union = Ω), then: P(A) = Σᵢ P(A | Bᵢ) · P(Bᵢ)

This is extraordinarily useful — it lets you compute a hard unconditional probability by conditioning on a convenient partition and summing weighted pieces.

Classic example (Monty Hall): Three doors, one car, two goats. You pick door 1. Host opens door 3 (always a goat, and host knows). Should you switch to door 2? Conditional probability: P(car behind 2 | host opened 3) = 2/3. Switching wins with probability 2/3. This surprises most people, which is why working it out carefully with the definitions matters.

# Monty Hall via simulation to build intuition
import random

def monty_hall(switch: bool) -> bool:
    doors = [0, 0, 1]          # 1 = car
    random.shuffle(doors)
    pick = random.randrange(3)
    # Host reveals a goat door (not picked, not the car)
    goat_doors = [d for d in range(3) if d != pick and doors[d] == 0]
    revealed = random.choice(goat_doors)
    if switch:
        # Switch to the remaining unopened door
        pick = next(d for d in range(3) if d != pick and d != revealed)
    return doors[pick] == 1

N = 200_000
win_stay   = sum(monty_hall(switch=False) for _ in range(N)) / N
win_switch = sum(monty_hall(switch=True)  for _ in range(N)) / N
print(f"Win rate (stay):   {win_stay:.4f}  (expected 1/3 ≈ 0.3333)")
print(f"Win rate (switch): {win_switch:.4f}  (expected 2/3 ≈ 0.6667)")

Two events A and B are independent if knowing one gives you no information about the other:

P(A ∩ B) = P(A) · P(B)

Equivalently (when P(B) > 0): P(A | B) = P(A). The conditional probability equals the unconditional — learning B happened tells you nothing new about A.

Mutual independence for a collection A₁, …, Aₙ: every subset S ⊆ {1,…,n} must satisfy P(⋂ᵢ∈S Aᵢ) = Πᵢ∈S P(Aᵢ). Pairwise independence does not imply mutual independence — this is a subtle but critical distinction. Example: flip two fair coins (H₁, H₂). Let A = {H₁}, B = {H₂}, C = {H₁ ≠ H₂} (exactly one head). Any two of these are pairwise independent, but all three together are not mutually independent (A ∩ B ∩ C = ∅, yet P(A)·P(B)·P(C) = 1/8 ≠ 0).

Independence in CS: Randomized algorithms depend heavily on independence. When we say "hash functions behave randomly," we usually mean k-wise independence: any k keys have independently uniform hash values. True independence would require exponentially large hash tables; k-wise independent families are much cheaper to construct and sufficient for most algorithms.

Conditional independence: A and B are conditionally independent given C if P(A ∩ B | C) = P(A | C) · P(B | C). This is the structural backbone of Bayesian networks: a child is conditionally independent of all non-descendants given its parents.

# Verify pairwise-but-not-mutually-independent example
from itertools import product

omega = list(product([0,1],[0,1]))  # (coin1, coin2)
P = {o: 0.25 for o in omega}       # uniform

A = {o for o in omega if o[0] == 1}       # coin1 = H
B = {o for o in omega if o[1] == 1}       # coin2 = H
C = {o for o in omega if o[0] != o[1]}    # exactly one H

def prob(event): return sum(P[o] for o in event)

print(f"P(A)={prob(A)}, P(B)={prob(B)}, P(C)={prob(C)}")
print(f"P(A∩B)={prob(A&B)} vs P(A)P(B)={prob(A)*prob(B)}  → pairwise indep")
print(f"P(A∩C)={prob(A&C)} vs P(A)P(C)={prob(A)*prob(C)}  → pairwise indep")
print(f"P(B∩C)={prob(B&C)} vs P(B)P(C)={prob(B)*prob(C)}  → pairwise indep")
print(f"P(A∩B∩C)={prob(A&B&C)} vs P(A)P(B)P(C)={prob(A)*prob(B)*prob(C)} → NOT mutually indep")

Bayes' theorem is perhaps the single most important formula in all of applied probability. It tells you how to update a prior belief in light of new evidence:

P(H | E) = P(E | H) · P(H) / P(E)

where:

• P(H) is the prior probability of hypothesis H before seeing evidence.
• P(E | H) is the likelihood — how probable is the evidence if H is true?
• P(E) is the marginal likelihood (normalizing constant), computed via total probability: P(E) = P(E|H)·P(H) + P(E|Hᶜ)·P(Hᶜ).
• P(H | E) is the posterior — our updated belief after seeing E.

Derivation: Pure algebra from the definition of conditional probability: P(H | E) = P(H ∩ E)/P(E) = [P(E | H)·P(H)] / P(E). That is all.

Worked example (medical test): A disease affects 1% of the population. A test is 99% sensitive (P(positive | sick) = 0.99) and 99% specific (P(negative | healthy) = 0.99). You test positive. P(sick | positive) = ?

P(positive) = 0.99 × 0.01 + 0.01 × 0.99 = 0.0198 P(sick | positive) = (0.99 × 0.01) / 0.0198 ≈ 0.5

You are only 50% likely to be sick! The low prior (1%) dominates. This surprises virtually everyone and is clinically important.

Sequential Bayes: The posterior from one update becomes the prior for the next. This is how a spam filter improves with each new email it classifies — and exactly how the Kalman filter works.

Naive Bayes classifier: Assume features are conditionally independent given the class. Then P(class | features) ∝ P(class) · Π P(featureᵢ | class). The "naive" independence assumption is almost always wrong, yet the classifier often works well in practice — a beautiful demonstration that a wrong model can still be useful.

def bayes_update(prior: float, likelihood: float, likelihood_complement: float) -> float:
    """Return P(H|E) given prior P(H), P(E|H), P(E|H^c)."""
    p_evidence = likelihood * prior + likelihood_complement * (1 - prior)
    return (likelihood * prior) / p_evidence

# Medical test example
prior   = 0.01    # 1% prevalence
lik_h   = 0.99    # sensitivity: P(pos | sick)
lik_hc  = 0.01    # 1-specificity: P(pos | healthy)

posterior = bayes_update(prior, lik_h, lik_hc)
print(f"P(sick | positive test) = {posterior:.4f}")  # ≈ 0.5

# Sequential updating: test positive twice (independent tests)
posterior2 = bayes_update(posterior, lik_h, lik_hc)
print(f"P(sick | two positive tests) = {posterior2:.4f}")  # ≈ 0.98

A random variable (RV) X is a function X: Ω → ℝ that assigns a real number to each outcome. It is a way to summarize outcomes numerically without specifying the full sample space.

Notation: X(ω) is the value of X at outcome ω. We write P(X = k) as shorthand for P({ω ∈ Ω : X(ω) = k}).

Probability mass function (PMF): For a discrete RV, pₓ(k) = P(X = k). It must satisfy:

• pₓ(k) ≥ 0 for all k.
• Σₖ pₓ(k) = 1.

Key discrete distributions:

Poisson as limit of Binomial: As n → ∞ and p → 0 with np = λ fixed, Binomial(n,p) → Poisson(λ). This is why the Poisson models "rare events in large populations": hash collisions, typos per page, radioactive decays per second, HTTP requests per millisecond.

Geometric is memoryless: P(X > m+n | X > m) = P(X > n). If you have been flipping and not seen a head yet, the coin has "no memory" of its previous failures. The geometric is the only discrete memoryless distribution.

import math
from collections import Counter
import random

# Poisson PMF
def poisson_pmf(lam: float, k: int) -> float:
    return math.exp(-lam) * (lam ** k) / math.factorial(k)

# Verify Binomial → Poisson convergence
lam = 3.0
for n in [10, 100, 1000, 10_000]:
    p = lam / n
    # Binomial P(X=3)
    binom = math.comb(n, 3) * p**3 * (1-p)**(n-3)
    pois  = poisson_pmf(lam, 3)
    print(f"n={n:6d}: Binom={binom:.6f}  Poisson={pois:.6f}  diff={abs(binom-pois):.2e}")

# Geometric memoryless check
def geometric_sample(p=0.3):
    k = 1
    while random.random() > p:
        k += 1
    return k

samples = [geometric_sample() for _ in range(200_000)]
# P(X>5) empirical
p_gt5 = sum(1 for s in samples if s > 5) / len(samples)
# P(X>5+3 | X>5) = P(X>3)
cond = sum(1 for s in samples if s > 8) / sum(1 for s in samples if s > 5)
print(f"P(X>3) = {sum(1<s<=5 for s in range(1,1000))}")
print(f"P(X>5)={p_gt5:.4f}   P(X>8|X>5)={cond:.4f}   (should be approx equal)")

The expectation (expected value, mean) of a discrete random variable X is:

E[X] = Σₖ k · P(X = k)

It is the probability-weighted average of all possible values. Concretely: if you repeat the experiment millions of times, E[X] is the long-run average value of X.

LOTUS (Law of the Unconscious Statistician): For a function g of X: E[g(X)] = Σₖ g(k) · P(X = k)

You do not need the distribution of g(X); just use the distribution of X.

Expected values of key distributions:

• E[Bernoulli(p)] = p
• E[Binomial(n,p)] = np
• E[Geometric(p)] = 1/p
• E[Poisson(λ)] = λ

Linearity of Expectation — the single most powerful tool in probabilistic analysis:

E[X + Y] = E[X] + E[Y], and E[cX] = c·E[X] for any constant c.

This holds for ALL random variables X and Y — even dependent ones.

This is startling. You never need to know the joint distribution of X and Y to compute E[X+Y]. This is not obvious, but follows immediately from the definition: E[X+Y] = Σ_ω (X(ω)+Y(ω))·P(ω) = Σ_ω X(ω)·P(ω) + Σ_ω Y(ω)·P(ω) = E[X] + E[Y].

Indicator random variables: Define Iₐ = 1 if event A occurs, 0 otherwise. Then E[Iₐ] = P(A). Combined with linearity, this is a killer technique: decompose a complicated count X = I₁ + I₂ + … + Iₙ into a sum of indicators, take expectations of each separately.

Example — expected number of hash collisions: n balls thrown uniformly into n bins. Let Iᵢⱼ = 1 if balls i and j collide. E[collisions] = E[Σᵢ<ⱼ Iᵢⱼ] = C(n,2) · (1/n) = (n−1)/2. No need to track dependencies at all.

# Linearity of expectation: count fixed points in a random permutation
# ("derangement" problem) — E[fixed points] = 1, by linearity!
import random

def count_fixed_points(n=10):
    perm = list(range(n))
    random.shuffle(perm)
    return sum(1 for i, x in enumerate(perm) if i == x)

n = 10
N = 300_000
fp_counts = [count_fixed_points(n) for _ in range(N)]
print(f"E[fixed points in perm of {n}] ≈ {sum(fp_counts)/N:.4f}  (theory: 1.0)")

# Proof via linearity: let I_i = 1 if element i is fixed
# E[I_i] = P(position i is fixed) = 1/n
# E[fixed points] = sum_{i=1}^{n} E[I_i] = n * (1/n) = 1
# No need to know distribution of the total count!

# Hash collision example
def expected_collisions_sim(n=100, trials=50_000):
    total = 0
    for _ in range(trials):
        bins = [random.randrange(n) for _ in range(n)]
        # count pairs that collide
        from collections import Counter
        c = Counter(bins)
        total += sum(v*(v-1)//2 for v in c.values())
    return total / trials

print(f"E[collisions, n=100]: sim={expected_collisions_sim():.2f}  theory={(100-1)/2:.2f}")

Expectation tells you the center of a distribution. Variance measures how spread out it is:

Var(X) = E[(X − E[X])²] = E[X²] − (E[X])²

The second form (computational formula) is usually easier. The standard deviation σ = √Var(X) has the same units as X.

Key variances:

• Var(Bernoulli(p)) = p(1−p)
• Var(Binomial(n,p)) = np(1−p)
• Var(Geometric(p)) = (1−p)/p²
• Var(Poisson(λ)) = λ

Variance is NOT linear in general — the crucial distinction from expectation: Var(X + Y) = Var(X) + Var(Y) + 2·Cov(X, Y)

where Cov(X,Y) = E[XY] − E[X]·E[Y].

For independent X and Y: Cov(X,Y) = 0, so Var(X+Y) = Var(X) + Var(Y).

Scaling: Var(cX) = c²·Var(X). Standard deviation scales linearly: σ(cX) = |c|·σ(X).

Why variance matters in CS: If X is a randomized algorithm's runtime, high variance means the algorithm is unpredictable even if its average is good. Quicksort with random pivots has E[comparisons] = O(n log n) and reasonably small variance — you can make strong guarantees. Variance analysis is also central to concentration inequalities (Chebyshev, Chernoff bounds), which give you tail probability bounds: how often will X deviate far from its mean?

Chebyshev's inequality: P(|X − E[X]| ≥ t) ≤ Var(X) / t². Weak but universal — it holds for any random variable with finite variance.

import statistics
import random

# Verify computational formula for variance on Binomial(n=20, p=0.3)
n, p = 20, 0.3
N = 200_000
samples = [sum(random.random() < p for _ in range(n)) for _ in range(N)]

emp_var = statistics.variance(samples)
theory_var = n * p * (1 - p)
print(f"Var(Bin(20,0.3)): empirical={emp_var:.4f}  theory={theory_var:.4f}")

# Computational formula: E[X^2] - (E[X])^2
emp_mean  = statistics.mean(samples)
emp_esq   = sum(x**2 for x in samples) / N
var_comp  = emp_esq - emp_mean**2
print(f"Computational formula check: {var_comp:.4f}")

# Chebyshev bound check
mean = n * p  # = 6
std  = theory_var**0.5
t = 2 * std
outside = sum(1 for s in samples if abs(s - mean) >= t) / N
chebyshev_bound = theory_var / t**2
print(f"P(|X-E[X]| >= 2σ): empirical={outside:.4f}  Chebyshev bound={chebyshev_bound:.4f}")

The union bound (Boole's inequality) is the simplest possible tool for bounding the probability of a bad event:

P(A₁ ∪ A₂ ∪ … ∪ Aₙ) ≤ Σᵢ P(Aᵢ)

Proof: Immediate from axiom 3. Countable additivity says equality holds for disjoint events. Union can only decrease probability (events overlap), so the sum over-counts, giving an upper bound.

This bound is often called the "first moment method" because P(⋃ Aᵢ) ≤ E[number of events that occur] by linearity (the RHS is exactly E[I₁ + I₂ + … + Iₙ] = Σ P(Aᵢ)).

When is it useful? When you want to show a bad event is unlikely, and the individual bad-event probabilities are small. It ignores all positive correlations between events, so it can be loose — but it is:

1. Always correct (it is always an upper bound).
2. Tight when events are nearly disjoint (almost no overlap).
3. Trivial to compute (just sum up individual probabilities).

Example — no collision in a hash table: Hash n elements into a table of size m. Let Aᵢⱼ be the event that elements i and j collide. P(Aᵢⱼ) = 1/m. By union bound: P(any collision) ≤ C(n,2)/m = n(n−1)/(2m). For m = n², this is at most (n−1)/(2n) < 1/2. So with a table of size n², collisions are unlikely!

Example — randomized algorithms: Suppose a Monte Carlo algorithm succeeds with probability ≥ 3/4 on each trial. Run it independently k times. P(all fail) ≤ (1/4)ᵏ (by independence, not union bound — but union bound over the k failure events gives the same order). For k = 100, this is essentially zero.

In complexity theory: The union bound is used constantly in probabilistic proofs: "there are at most M bad events, each with probability at most p, so the probability that any bad event occurs is at most Mp."

# Demonstrate union bound for birthday problem / hash collision
import random
import math

def collision_prob_empirical(n_elements: int, table_size: int, trials: int = 200_000) -> float:
    """Empirical P(at least one collision) when hashing n elements into m slots."""
    collisions = 0
    for _ in range(trials):
        slots = [random.randrange(table_size) for _ in range(n_elements)]
        if len(set(slots)) < n_elements:
            collisions += 1
    return collisions / trials

n, m = 20, 400  # n=20 elements, m=n^2 slots
emprical = collision_prob_empirical(n, m)
union_bound = math.comb(n, 2) / m   # C(n,2)/m
print(f"n={n}, m={m}=n^2")
print(f"Empirical P(collision) = {emprical:.4f}")
print(f"Union bound            = {union_bound:.4f}")
print(f"Bound is {'valid' if union_bound >= emprical else 'VIOLATED (bug!)'}")

# Tightness: how loose is the bound?
print(f"Ratio (bound/actual)   = {union_bound/emprical:.2f}x")

The conditional expectation E[X | Y = y] is the expected value of X, given that we know Y = y:

E[X | Y = y] = Σₓ x · P(X = x | Y = y)

This is itself a function of y. If we think of it as a random variable (before observing Y), we write E[X | Y], which is a random variable that equals E[X | Y = y] whenever Y = y.

Tower Property (Law of Total Expectation): E[X] = E[E[X | Y]] = Σᵧ E[X | Y = y] · P(Y = y)

This is the expectation analog of the law of total probability. It says: if you first compute the conditional expectation given each value of Y, then average those over Y, you get E[X] back. This is enormously useful when computing E[X] directly is hard but conditioning on Y simplifies it.

Example — expected length of a random walk: Let X = number of steps a random walk takes to return to the origin. Condition on the first step — by symmetry and the tower property, you can set up a recursion.

Wald's identity (a consequence): If X₁, X₂, … are i.i.d. with mean μ, and N is a stopping time (a positive integer-valued RV, independent of future Xᵢ), then: E[X₁ + X₂ + … + X_N] = E[N] · μ

This is why quicksort analysis is tractable: the expected number of comparisons at a subproblem of size k is E[comparisons | size k] = k − 1, and we can tower over the random partition.

CS application — randomized divide-and-conquer: In randomized Quicksort, let T(n) = expected comparisons on n elements. Conditioning on the rank of the chosen pivot: T(n) = (n−1) + (1/n) Σₖ₌₁ⁿ [T(k−1) + T(n−k)] Solving this gives T(n) = 2n ln n + O(n). The tower property is the analytical tool that makes this tractable.

# Verify Wald's identity: sum of N i.i.d. geometric RVs
# Each trial: flip until first head (p=0.5), then flip N more groups
# N itself is Geometric(0.3)

import random
import statistics

def geometric(p):
    k = 0
    while random.random() >= p:
        k += 1
    return k + 1

p_N = 0.3  # distribution of N
p_X = 0.5  # distribution of each X_i
E_N = 1 / p_N  # = 3.33...
E_X = 1 / p_X  # = 2.0

trials = 200_000
totals = []
for _ in range(trials):
    N = geometric(p_N)
    total = sum(geometric(p_X) for _ in range(N))
    totals.append(total)

emp = statistics.mean(totals)
wald = E_N * E_X
print(f"E[X_1+...+X_N]: empirical={emp:.4f}  Wald's identity={wald:.4f}")