7. Discrete Probability

Probability measures how likely an uncertain outcome is, and in the discrete setting it is counting (Chapter 5) with weights attached. This chapter builds the theory from sample spaces, axioms, conditional probability, and Bayes' theorem to random variables, expectation, and variance, and then to the standard distributions, the tail bounds of Markov and Chebyshev, the probabilistic method, and randomized algorithms. These tools explain why hash tables are fast, how far a Miller–Rabin (4.4) answer can be trusted, how spam filters and A/B tests reach conclusions, and why randomized quicksort runs in O(n log n) expected time on every input.

<a id="71-probability-fundamentals"></a>

7.1 Probability Fundamentals

This subchapter defines a probability precisely: a sample space of outcomes, events as subsets of it, and a function that obeys three axioms. From these it derives the rules for complements, unions, conditioning, and independence.

Sample Spaces and Events

Why. Before you can compute a probability, you have to say exactly what can happen. Most probability mistakes, including several famous paradoxes, come from a vague or wrong description of the possible outcomes. Treating the eleven sums of two dice as equally likely is one such mistake, and it gives wrong answers even though the arithmetic that follows is correct. The sample space is the data model of a probability problem. Get it right and the rest is calculation.

What. An experiment is any process whose result is uncertain: rolling dice, shuffling a deck, hashing a random key, sending a packet over a lossy link.

• The sample space Ω is the set of all possible outcomes. Outcomes must be mutually exclusive (exactly one happens) and exhaustive (one always happens).
• An event is a subset E ⊆ Ω. The event occurs if the actual outcome lies in E.
• {ω} for a single outcome ω is an elementary event. Ω itself is the certain event, and ∅ is the impossible event.

In this chapter, Ω is finite or countably infinite (2.3). That is what "discrete" means.

Because events are sets, the set operations of 2.1 express logical combinations of events:

In words	Set notation
A or B (or both)	A ∪ B
A and B	A ∩ B
not A	Ā = Ω − A
A but not B	A − B
A and B cannot both happen (mutually exclusive)	A ∩ B = ∅
whenever A happens, B happens	A ⊆ B

Example: two dice. Record the result as an ordered pair (first die, second die). Then Ω = {1, …, 6} × {1, …, 6} is a Cartesian product (2.1) with 36 outcomes. The event "the sum is 7" is the anti-diagonal of the grid:

                second die
              1   2   3   4   5   6
           1  .   .   .   .   .   7
  first    2  .   .   .   .   7   .
  die      3  .   .   .   7   .   .
           4  .   .   7   .   .   .
           5  .   7   .   .   .   .
           6  7   .   .   .   .   .

  "sum is 7" = {(1,6), (2,5), (3,4), (4,3), (5,2), (6,1)}      6 of 36 outcomes

Choosing the sample space. One experiment can be modeled by several sample spaces, which differ in how much detail they record:

Sample space for two dice	Size	Outcomes equally likely?	Can it answer "is the first die a 6"?
Ordered pairs (a, b)	36	Yes	Yes
Unordered pairs {a, b}	21	No: {1, 2} is twice as likely as {1, 1}	No
Sums {2, …, 12}	11	No: 7 is six times as likely as 2	No

A coarse model can always be computed from a fine one, but not the other way around. Start with the finest description that is still easy to count, usually one with equally likely outcomes. Some experiments need an infinite Ω: tossing a coin until the first head has Ω = {H, TH, TTH, TTTH, …}, which is countably infinite.

from itertools import product

omega = set(product(range(1, 7), repeat=2))           # 36 ordered pairs (first die, second die)
sum7    = {w for w in omega if w[0] + w[1] == 7}
doubles = {w for w in omega if w[0] == w[1]}
first6  = {w for w in omega if w[0] == 6}

print(len(omega), len(sum7), len(doubles))            # 36 6 6
print(sum7 & doubles)                                 # set(): "sum 7" and "doubles" are mutually exclusive
print(sorted(sum7 & first6))                          # [(6, 1)]
print(len(sum7 | doubles), len(omega - doubles))      # 12 30

Key takeaways

• A sample space Ω lists every possible outcome. Outcomes are mutually exclusive and exhaustive.
• Events are subsets of Ω, so "or", "and", and "not" are union, intersection, and complement.
• Mutually exclusive means disjoint: A ∩ B = ∅.
• Prefer a fine-grained sample space with equally likely outcomes, such as ordered pairs for dice rather than sums.

Practice

1. (Beginner) Three coins are tossed. Write out Ω, then the events A = "exactly two heads" and B = "the first coin is heads". List A ∩ B, A ∪ B, and Ā.
2. (Intermediate) A client retries a failing request until it succeeds, making at most 3 attempts. Describe Ω as sequences of S/F results, and write the event "the client gives up". Is {1, 2, 3, gave up} (the number of attempts that ended in success, or giving up) also a valid sample space? What information does it lose?
3. (Interview) You want to test a load balancer that sends each of k requests to one of n servers at random. Describe a sample space for one test run and give its size. Write the event "server 1 receives no requests" as a set. Hint: an outcome is a function from requests to servers (2.2). Count them with the product rule (5.1).

Probability Axioms

Why. People mean different things by "probability": a long-run frequency ("this coin lands heads half the time"), a degree of belief ("70% chance of rain tomorrow"), or a ratio of symmetric cases ("one of six faces"). All three interpretations follow the same rules. In 1933, Kolmogorov wrote these rules down as three axioms, which turned probability into an ordinary branch of mathematics. Every other rule in this chapter is a theorem derived from them.

What. A probability function P assigns a number to each event so that:

1. Nonnegativity: P(E) ≥ 0 for every event E.
2. Normalization: P(Ω) = 1.
3. Additivity: if E₁, E₂, … are pairwise disjoint, then P(E₁ ∪ E₂ ∪ ⋯) = P(E₁) + P(E₂) + ⋯. This must hold for countably many events, not only for two.

The pair (Ω, P) is a probability space. In the discrete case, P is completely determined by the probabilities of the individual outcomes. Choose any weights p(ω) ≥ 0 with Σ_{ω ∈ Ω} p(ω) = 1, and define

  P(E) = Σ_{ω ∈ E} p(ω)

Analogy. Picture one kilogram of sand spread over the outcomes. The probability of an event is the amount of sand lying on it. Sand cannot be negative (axiom 1), there is exactly one kilogram in total (axiom 2), and the sand on separate piles adds up (axiom 3).

Consequences. Each of these follows from the axioms in a line or two:

Rule	Reason
P(Ā) = 1 − P(A)	A and Ā are disjoint with union Ω, so P(A) + P(Ā) = P(Ω) = 1
P(∅) = 0	∅ is the complement of Ω
If A ⊆ B then P(A) ≤ P(B)	B is the disjoint union of A and B − A, and P(B − A) ≥ 0
0 ≤ P(A) ≤ 1	A ⊆ Ω and monotonicity
P(A − B) = P(A) − P(A ∩ B)	A is the disjoint union of A − B and A ∩ B

Example 1: a loaded die. Faces 1 to 5 each have probability 1/10 and face 6 has probability 1/2. The weights are nonnegative and sum to 1, so this is a valid probability space. P(even) = 1/10 + 1/10 + 1/2 = 7/10.

Example 2: an infinite space. Toss a fair coin until the first head, and let the outcome be the number of tosses k. Put p(k) = 1/2ᵏ. The geometric series (3.1) gives Σ_{k ≥ 1} 1/2ᵏ = 1, so this is valid too, and P(k is even) = 1/4 + 1/16 + ⋯ = 1/3.

from fractions import Fraction as F

def make_space(weights):
    """weights: dict outcome -> probability. Checks the axioms for a finite space."""
    assert all(p >= 0 for p in weights.values()), "nonnegativity violated"
    assert sum(weights.values()) == 1, "probabilities must sum to 1"
    return weights

def P(space, event):                                  # event: a predicate on outcomes
    return sum((p for w, p in space.items() if event(w)), F(0))

die = make_space({1: F(1, 10), 2: F(1, 10), 3: F(1, 10), 4: F(1, 10), 5: F(1, 10), 6: F(1, 2)})
print(P(die, lambda w: w % 2 == 0))                   # 7/10
print(P(die, lambda w: w % 2 == 1))                   # 3/10   = 1 − P(even)
print(P(die, lambda w: w > 6))                        # 0      the impossible event

# Countably infinite space: first head on toss k with probability 1/2^k. Partial sums approach 1.
print(float(sum(F(1, 2**k) for k in range(1, 31))))  # 0.9999999990686774
print(float(sum(F(1, 2**k) for k in range(2, 61, 2))))   # 0.3333333333333333   P(k is even)

Key takeaways

• Three axioms: nonnegativity, P(Ω) = 1, and additivity over disjoint events.
• In a discrete space, it is enough to give each outcome a weight. The weights must be nonnegative and sum to 1.
• The complement rule, monotonicity, and 0 ≤ P ≤ 1 are theorems, not extra assumptions.
• Frequentist and Bayesian interpretations disagree about what probability means, but both use these rules.

Practice

1. (Beginner) Using only the axioms and the complement rule, prove that P(A ∪ B) ≤ P(A) + P(B) for any two events. Hint: write A ∪ B as the disjoint union of A and B − A.
2. (Intermediate) For which constant c is p(k) = c/k², for k = 1, 2, 3, …, a valid assignment of probabilities? Is there any c that makes p(k) = c/k valid? Hint: Σ 1/k² = π²/6, and the harmonic series diverges (3.1).
3. (Interview) Why is it impossible to "pick a positive integer uniformly at random"? Hint: suppose every integer has the same probability q. Consider q = 0 and q > 0 separately, and apply countable additivity.

Equally Likely Outcomes

Why. Many experiments are symmetric: fair coins, fair dice, well-shuffled decks, a random key under a good hash function. When no outcome is favored, a probability is simply a ratio of two counts, and every technique from Chapter 5 becomes a technique for computing probabilities.

What. On a finite sample space, the uniform probability gives every outcome the same weight 1/|Ω|. Then, for every event E,

  P(E) = |E| / |Ω| = (number of favorable outcomes) / (number of possible outcomes)

This formula is only valid when the outcomes really are equally likely. In 1754, d'Alembert argued that the probability of at least one head in two tosses is 2/3. He listed the outcomes as "head on the first toss", "head on the second toss", and "no head", and treated them as equally likely:

Model	Sample space	Equally likely?	P(at least one head)
Ordered tosses	{HH, HT, TH, TT}	Yes	3/4 (correct)
d'Alembert	{H, TH, TT}	No: H covers both HH and HT	2/3 (wrong)

How. Count |Ω| and |E| in the same way: both ordered, or both unordered. The counting tools of 5.2 do the rest.

Example 1: poker. A hand is 5 cards out of 52, so there are C(52, 5) = 2 598 960 equally likely unordered hands. To count full houses (three cards of one rank and two of another), choose the rank of the triple (13 ways), its suits (C(4, 3)), the rank of the pair (12 ways), and its suits (C(4, 2)):

  P(full house) = 13 · 4 · 12 · 6 / 2 598 960 = 3744 / 2 598 960 ≈ 0.00144

Example 2: a lottery. If 6 numbers are drawn from 49, a single ticket wins with probability 1/C(49, 6) = 1/13 983 816.

Example 3: sampling without replacement (hypergeometric). A batch of N items contains m defective ones. You inspect r items chosen at random without replacement. Every r-subset is equally likely, so

  P(exactly k defective in the sample) = C(m, k) · C(N − m, r − k) / C(N, r)

The numerator chooses k of the defective items and r − k of the good ones. Vandermonde's identity (5.3) shows that these probabilities sum to 1 over all k. For example, take 20 chips of which 4 are defective, and inspect 3: P(no defective found) = C(16, 3)/C(20, 3) = 560/1140 = 28/57 ≈ 0.49. An inspection of 3 chips misses all four defective chips almost half the time.

from itertools import combinations
from math import comb
from fractions import Fraction as F

# Brute force: label the chips 0..19 so they are distinguishable; chips 0..3 are defective
samples = list(combinations(range(20), 3))
clean = sum(all(chip >= 4 for chip in s) for s in samples)
print(len(samples), F(clean, len(samples)))           # 1140 28/57
print(F(comb(16, 3), comb(20, 3)))                    # 28/57   the hypergeometric formula agrees

full_house = 13 * comb(4, 3) * 12 * comb(4, 2)
print(full_house, F(full_house, comb(52, 5)))         # 3744 6/4165
print(round(full_house / comb(52, 5), 5))             # 0.00144

Key takeaways

• When all outcomes are equally likely, P(E) = |E|/|Ω|, so a probability is a ratio of two counts.
• Check that the outcomes really are equally likely. Sums of dice and d'Alembert's three outcomes are not.
• Count the event and the sample space in the same way, both ordered or both unordered.
• Sampling without replacement gives the hypergeometric formula C(m, k)C(N − m, r − k)/C(N, r).

Practice

1. (Beginner) Two fair dice are rolled. Find P(sum = 8) and P(at least one 6).
2. (Intermediate) A password consists of 8 characters chosen uniformly and independently from 62 (letters and digits). Find the probability that it contains at least one digit, and the probability that all 8 characters are different.
3. (Interview) A deck is shuffled uniformly. What is the probability that the top two cards are both aces? Solve it twice, once by counting ordered outcomes and once by counting unordered ones, and confirm that the answers agree. Hint: ordered: 4 · 3 favorable out of 52 · 51. Unordered: C(4, 2) out of C(52, 2).

Complements and Unions

Why. Questions of the form "at least one ..." are everywhere: at least one server fails, at least one hash collision occurs, at least one test is flaky. Counting "at least one" directly means adding up many overlapping cases. Its complement, "none", is usually a single clean case. When the complement does not help, inclusion–exclusion (5.5) and the union bound handle unions of overlapping events.

What.

• Complement rule: P(Ā) = 1 − P(A), so P(at least one) = 1 − P(none).
• Union of two events: P(A ∪ B) = P(A) + P(B) − P(A ∩ B). The outcomes in A ∩ B are counted twice by P(A) + P(B), so they are subtracted once:

  ┌──────────── Ω ─────────────────────┐
  │     ┌───── A ─────┐                │
  │     │        ┌────┼──── B ─────┐   │
  │     │  A − B │A∩B │   B − A    │   │      P(A) + P(B) counts A ∩ B twice
  │     │        └────┼────────────┘   │
  │     └─────────────┘                │
  └────────────────────────────────────┘

• Inclusion–exclusion (5.5) extends this to any number of events. For three events: P(A ∪ B ∪ C) = P(A) + P(B) + P(C) − P(A∩B) − P(A∩C) − P(B∩C) + P(A∩B∩C).
• Union bound (Boole's inequality): P(A₁ ∪ ⋯ ∪ Aₙ) ≤ P(A₁) + ⋯ + P(Aₙ). It holds for any events, with no assumption of independence. It is accurate when the events are rare and rarely overlap, which is exactly the "something bad happens" situation in system and algorithm analysis.

Example 1: the Chevalier de Méré (1654). A gambler noticed that betting on "at least one 6 in 4 rolls of a die" won slightly more often than it lost, but that betting on "at least one double 6 in 24 rolls of two dice" lost slightly more often than it won. He found this puzzling, because 4/6 = 24/36. The complement explains it. There are 5⁴ sequences of 4 rolls with no 6, out of 6⁴, and 35²⁴ sequences of 24 double rolls with no double 6, out of 36²⁴:

  P(at least one 6 in 4 rolls)              = 1 − (5/6)⁴    ≈ 0.5177
  P(at least one double 6 in 24 rolls)      = 1 − (35/36)²⁴ ≈ 0.4914

The ratio argument "4/6 = 24/36" adds probabilities of overlapping events, which is the union bound used as if it were an equality.

Example 2: divisibility. Pick an integer uniformly from 1 to 100. P(divisible by 2 or 5) = 50/100 + 20/100 − 10/100 = 0.6.

Example 3: a union bound for reliability. A job runs on 10 machines, and each fails with probability 0.01. Whether or not failures are correlated, P(some machine fails) ≤ 10 · 0.01 = 0.1. If failures happen to be independent (see Independence below), the exact value is 1 − 0.99¹⁰ ≈ 0.0956, so the bound is close.

from fractions import Fraction as F

print(round(1 - (5/6)**4, 4), round(1 - (35/36)**24, 4))    # 0.5177 0.4914

N = range(1, 101)
A = {n for n in N if n % 2 == 0}
B = {n for n in N if n % 5 == 0}
print(F(len(A | B), 100), F(len(A) + len(B) - len(A & B), 100))   # 3/5 3/5

union_bound = 10 * 0.01
exact_if_independent = 1 - 0.99**10
print(union_bound, round(exact_if_independent, 4))                 # 0.1 0.0956

Key takeaways

• For "at least one", compute 1 − P(none).
• P(A ∪ B) = P(A) + P(B) − P(A ∩ B). Inclusion–exclusion extends this to more events.
• The union bound P(∪ Aᵢ) ≤ Σ P(Aᵢ) needs no independence. It is nearly exact when the events are rare.
• Adding probabilities of events that overlap overestimates the probability of their union.

Practice

1. (Beginner) A 4-digit PIN from 0000 to 9999 is chosen uniformly. What is the probability that some digit appears more than once?
2. (Intermediate) An integer is chosen uniformly from 1 to 1000. Use inclusion–exclusion to find the probability that it is divisible by 3, 5, or 7.
3. (Interview) A batch job runs on 200 machines. Each machine fails during the job with probability 10⁻⁴, and failures may be correlated (shared racks, shared power). Give an upper bound on the probability that the job sees any failure. Can you give a lower bound without further assumptions? Hint: the union bound needs no independence. For the lower bound, think about what happens if all failures occur together.

Conditional Probability

Why. New information changes probabilities. Once you know the first card dealt is an ace, the chance that the second is an ace drops. Once a test comes back positive, the chance that the patient is sick rises. Conditional probability is the rule for updating, and it also lets you compute the probability of a sequence of steps one step at a time.

What. For events A and B with P(B) > 0, the conditional probability of A given B is

  P(A | B) = P(A ∩ B) / P(B)

Conditioning on B restricts the sample space to B and rescales, so that B has total probability 1. Outcomes outside B are ruled out. Inside B, the outcomes keep their relative weights:

  before:  Ω = all 36 dice outcomes          after learning B = "first die is 3":
           every outcome has weight 1/36       Ω' = {(3,1), (3,2), …, (3,6)}, each with weight 1/6

With equally likely outcomes, this simplifies to P(A | B) = |A ∩ B| / |B|. For a fixed B, the function P(· | B) satisfies all three axioms, so the complement rule, inclusion–exclusion, and every other rule also hold for conditional probabilities.

Rearranging the definition gives the multiplication rule P(A ∩ B) = P(B) · P(A | B). Applying it repeatedly gives the chain rule:

  P(A₁ ∩ A₂ ∩ ⋯ ∩ Aₙ) = P(A₁) · P(A₂ | A₁) · P(A₃ | A₁ ∩ A₂) ⋯ P(Aₙ | A₁ ∩ ⋯ ∩ Aₙ₋₁)

Example 1: order matters. Roll two dice. P(sum = 8 | first die is 3) = 1/6, since only (3, 5) works out of 6 outcomes. On the other hand, P(first die is 3 | sum = 8) = 1/5, since (3, 5) is one of the 5 outcomes with sum 8. In general, P(A | B) and P(B | A) are different numbers. Confusing them is the "prosecutor's fallacy", and Bayes' theorem (7.2) is the correct way to convert one into the other.

Example 2: the chain rule for card draws. The probability that the first two cards dealt are both aces is P(1st ace) · P(2nd ace | 1st ace) = 4/52 · 3/51 = 1/221.

Example 3: two children. A family has two children, and each child is independently a boy or a girl with probability 1/2. Take Ω = {BB, BG, GB, GG}, listing the older child first.

| What you learn | Restricted sample space | P(both boys | information) | |---|---|---| | The older child is a boy | {BB, BG} | 1/2 | | At least one child is a boy | {BB, BG, GB} | 1/3 |

The two statements sound alike, but they rule out different outcomes. What matters is exactly which event you condition on.

from itertools import product
from fractions import Fraction as F

def cond(omega, A, B):                       # P(A | B) on an equally likely sample space
    in_B = [w for w in omega if B(w)]
    return F(sum(1 for w in in_B if A(w)), len(in_B))

dice = list(product(range(1, 7), repeat=2))
print(cond(dice, lambda w: sum(w) == 8, lambda w: w[0] == 3))    # 1/6
print(cond(dice, lambda w: w[0] == 3, lambda w: sum(w) == 8))    # 1/5

kids = list(product("BG", repeat=2))                              # (older, younger)
both = lambda w: w == ("B", "B")
print(cond(kids, both, lambda w: w[0] == "B"))                   # 1/2
print(cond(kids, both, lambda w: "B" in w))                      # 1/3

print(F(4, 52) * F(3, 51))                                       # 1/221   chain rule

Key takeaways

• P(A | B) = P(A ∩ B)/P(B): restrict to B, then rescale.
• P(· | B) is itself a probability function, so every rule still applies after conditioning.
• Multiplication rule and chain rule: P(A ∩ B) = P(B)P(A | B). Use them for sequential experiments.
• P(A | B) ≠ P(B | A) in general. Be precise about what is known.

Practice

1. (Beginner) Two cards are drawn without replacement. Find P(second card is a heart | first card is a heart) and P(both cards are hearts).
2. (Intermediate) Two dice are rolled. Compute P(at least one 6 | sum = 9) and P(sum = 9 | at least one 6), and explain why they differ.
3. (Interview) A family has two children, and you learn that at least one of them is a boy born on a Tuesday. What is the probability that both are boys? Assume sexes and weekdays are uniform and independent. Hint: the sample space has 14 × 14 = 196 equally likely outcomes. Count the outcomes with at least one Tuesday boy, taking care not to count the outcome "both are Tuesday boys" twice. The answer is not 1/3.

Independence

Why. Many systems are built from parts that do not affect each other: separate coin flips, separate dice, servers in different data centers. When parts are independent, probabilities multiply, and the reliability of a large system can be computed from the reliability of its parts. Assuming independence when it does not hold is also one of the most expensive mistakes in engineering. Disks from the same manufacturing batch fail together, and mortgage defaults rose together in 2008.

What. Events A and B are independent if

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

When P(B) > 0, this is equivalent to P(A | B) = P(A): learning that B happened does not change the probability of A. The definition is symmetric in A and B. If A and B are independent, then so are A and B̄, Ā and B, and Ā and B̄.

Events A₁, …, Aₙ are mutually independent if the product rule holds for every subcollection: P(A_{i₁} ∩ ⋯ ∩ A_{iₖ}) = P(A_{i₁}) ⋯ P(A_{iₖ}). Independence of every pair is not enough:

  Two fair coin tosses.   A = "first is H",   B = "second is H",   C = "the two tosses agree"

  P(A) = P(B) = P(C) = 1/2
  P(A∩B) = P(A∩C) = P(B∩C) = 1/4         every pair is independent
  P(A∩B∩C) = P(HH) = 1/4 ≠ 1/8           but A and B together determine C

Disjoint is not independent. If A ∩ B = ∅ and both events have positive probability, then P(A ∩ B) = 0 ≠ P(A)P(B). Disjoint events are strongly dependent: learning that one happened tells you for certain that the other did not.

Application: reliability. Suppose components fail independently, and component i works with probability pᵢ.

  series (all must work):          ──[p₁]──[p₂]──[p₃]──       P(works) = p₁ · p₂ · p₃

  parallel (one is enough):           ┌──[p₁]──┐
                                   ───┼──[p₂]──┼───          P(works) = 1 − (1 − p₁)(1 − p₂)(1 − p₃)
                                      └──[p₃]──┘

A request that passes through 5 services in series, each 99.9% available, succeeds with probability 0.999⁵ ≈ 0.995. Three independent replicas, each 99% available, give 1 − 0.01³ = 0.999999. Correlated failures, such as a shared network switch, break both calculations. This is why replicas are placed in separate failure domains.

from itertools import product
from fractions import Fraction as F
from math import prod

omega = list(product("HT", repeat=2))
def P(E): return F(sum(1 for w in omega if E(w)), len(omega))

A = lambda w: w[0] == "H"
B = lambda w: w[1] == "H"
C = lambda w: w[0] == w[1]
both = lambda E1, E2: (lambda w: E1(w) and E2(w))

print(all(P(both(X, Y)) == P(X) * P(Y) for X, Y in [(A, B), (A, C), (B, C)]))   # True: pairwise independent
print(P(lambda w: A(w) and B(w) and C(w)), P(A) * P(B) * P(C))               # 1/4 1/8: not mutually independent

series   = lambda ps: prod(ps)
parallel = lambda ps: 1 - prod(1 - p for p in ps)
print(round(series([0.999] * 5), 4), round(parallel([0.99] * 3), 6))         # 0.995 0.999999

Key takeaways

• A and B are independent iff P(A ∩ B) = P(A)P(B), which is equivalent to P(A | B) = P(A) when P(B) > 0.
• Mutual independence requires the product rule for every subcollection. Pairwise independence is weaker.
• Disjoint events with positive probability are dependent.
• Series systems multiply success probabilities, and parallel systems multiply failure probabilities. Both calculations assume independent failures.

Practice

1. (Beginner) Two dice are rolled. Are the events "the first die is even" and "the sum is 7" independent? What about "the first die is even" and "the sum is 8"?
2. (Intermediate) Prove that if A and B are independent, then Ā and B are independent. Hint: P(B) = P(A ∩ B) + P(Ā ∩ B).
3. (Interview) A checkout request calls 5 microservices in sequence, each with 99.9% availability. What is the end-to-end availability? How many independent replicas of each service would raise it above 99.99%, and why might the real improvement be smaller than your calculation predicts? Hint: each replicated service is a parallel subsystem inside a series system. Then think about which failures replicas share.

<a id="72-bayes-and-applications"></a>

7.2 Bayes and Applications

This subchapter puts conditional probabilities together. It splits a probability into cases, reverses a conditional with Bayes' theorem, and applies both tools to two famously counterintuitive problems.

Law of Total Probability

Why. Some probabilities are hard to compute directly but easy once you know which case you are in. What fraction of all shipped parts are defective? That is hard to say until you split the parts by the machine that made them, because each machine has a known defect rate. The law of total probability is the probabilistic version of the sum rule (5.1): split into cases, solve each case, and combine the results with weights.

What. Events B₁, …, Bₙ form a partition of Ω if they are pairwise disjoint, their union is Ω, and each has positive probability. Then, for any event A:

  P(A) = P(A | B₁)·P(B₁) + P(A | B₂)·P(B₂) + ⋯ + P(A | Bₙ)·P(Bₙ)

Proof. The sets A ∩ Bᵢ are disjoint and their union is A. By additivity, P(A) = Σ P(A ∩ Bᵢ), and by the multiplication rule, each term equals P(A | Bᵢ)P(Bᵢ). The most common case uses the partition {B, B̄}: P(A) = P(A | B)P(B) + P(A | B̄)P(B̄).

A tree diagram shows the computation. Each branch is labeled with a conditional probability, a leaf's probability is the product of the labels along its path, and P(A) is the sum of the leaves where A happens:

                      ┌── defective  0.01   → 0.50 × 0.01 = 0.005
          ┌─ M1 0.50 ─┤
          │           └── ok         0.99
          │           ┌── defective  0.02   → 0.30 × 0.02 = 0.006
  part ───┼─ M2 0.30 ─┤
          │           └── ok         0.98
          │           ┌── defective  0.05   → 0.20 × 0.05 = 0.010
          └─ M3 0.20 ─┤
                      └── ok         0.95
                                              P(defective) = 0.021

Example: conditioning on the first step. Total probability can also be applied to the first step of a process that repeats. In tennis, from "deuce" a player must win two points in a row to win the game. Suppose the player wins each point independently with probability p, and let w be the probability of winning the game from deuce. Split on the next two points:

  win, win    (prob p²)          → game won
  lose, lose  (prob (1 − p)²)    → game lost
  split       (prob 2p(1 − p))   → back to deuce, winning with probability w

  w = p² + 2p(1 − p)·w     ⇒     w = p² / (p² + (1 − p)²)

For p = 0.6, this gives w ≈ 0.692. A small advantage on each point grows into a larger advantage over the whole game. The same "condition on the first step, get an equation" method solves random walks, gambler's ruin, and expected waiting times (7.3).

import random

def total_probability(cases):            # cases: list of (P(B_i), P(A | B_i))
    assert abs(sum(pb for pb, _ in cases) - 1) < 1e-12, "the B_i must partition Ω"
    return sum(pb * pa for pb, pa in cases)

print(round(total_probability([(0.50, 0.01), (0.30, 0.02), (0.20, 0.05)]), 3))   # 0.021

def deuce_exact(p):
    return p * p / (p * p + (1 - p) ** 2)

def deuce_sim(p, games, rng):
    wins = 0
    for _ in range(games):
        lead = 0
        while abs(lead) < 2:                 # play points until someone is two points ahead
            lead += 1 if rng.random() < p else -1
        wins += lead == 2
    return wins / games

print(round(deuce_exact(0.6), 4))                            # 0.6923
print(round(deuce_sim(0.6, 200_000, random.Random(3)), 2))   # 0.69

Key takeaways

• If B₁, …, Bₙ partition Ω, then P(A) = Σ P(A | Bᵢ)P(Bᵢ), a weighted average of the conditional probabilities.
• Draw a tree: multiply along the branches, then add the leaves where A happens.
• Conditioning on the first step of a repeating process gives an equation, or a recurrence (6.1), for the unknown probability.

Practice

1. (Beginner) 60% of a site's traffic comes from mobile devices, with an error rate of 3%. The rest comes from desktops, with an error rate of 1%. What is the overall error rate?
2. (Intermediate) Urn 1 holds 3 red and 2 blue balls, and urn 2 holds 1 red and 4 blue balls. A random ball is moved from urn 1 to urn 2, and then a ball is drawn from urn 2. What is the probability that it is red?
3. (Interview) A gambler starts with k dollars and bets 1 dollar on fair coin flips until reaching either 0 or N dollars. What is the probability of reaching N? Hint: conditioning on the first flip gives pₖ = (pₖ₋₁ + pₖ₊₁)/2, with p₀ = 0 and p_N = 1. Solve this linear recurrence with the methods of 6.1.

Bayes' Theorem

Why. We often know how likely the evidence is under each explanation, but we want to know how likely each explanation is given the evidence. A test's manufacturer reports P(positive | sick), but a patient wants P(sick | positive). A spam filter can estimate P(word | spam) from training data, but it needs P(spam | words). Bayes' theorem reverses the direction of a conditional probability. It also guards against base-rate neglect, the common error of ignoring how rare the explanation was to begin with.

What. From P(A ∩ B) = P(A | B)P(B) = P(B | A)P(A):

                 P(A | B) · P(B)                          P(A | Bⱼ) · P(Bⱼ)
  P(B | A)  =  ─────────────────          P(Bⱼ | A) = ───────────────────────────
                     P(A)                                Σᵢ P(A | Bᵢ) · P(Bᵢ)

The second form uses the law of total probability for the denominator. Each piece has a name:

Term	Name	Meaning
P(B)	Prior	Belief in B before seeing the evidence
`P(A	B)`	Likelihood
P(A)	Evidence	Total probability of the evidence, over all explanations
`P(B	A)`	Posterior

Example: a medical test. A disease affects 1% of the population. A test detects it 99% of the time (sensitivity) and gives a false positive in 5% of healthy people. If someone tests positive:

                  0.99 × 0.01                    0.0099
  P(sick | +) = ─────────────────────────── = ──────── ≈ 0.167
                0.99 × 0.01 + 0.05 × 0.99      0.0594

Only about one positive in six is a real case. Natural frequencies make this easy to see. Imagine 10 000 people:

	Test positive	Test negative	Total
Sick	99	1	100
Healthy	495	9405	9900
Total	594	9406	10 000

The 495 false positives from the large healthy group outnumber the 99 true positives from the small sick group.

Odds form. Writing the theorem in terms of odds makes repeated updates simple:

  posterior odds = prior odds × likelihood ratio,      likelihood ratio = P(A | B) / P(A | B̄)

For the test, the prior odds are 1 : 99, and the likelihood ratio is 0.99/0.05 = 19.8. After one positive, the odds are 19.8 : 99 = 1 : 5, which is probability 1/6. A second, independent positive multiplies the odds by 19.8 again, giving 3.96 : 1, which is about 80%. Each posterior becomes the prior for the next piece of evidence.

Application: naive Bayes classification. To classify a message as spam or not ("ham"), compare P(spam) · P(words | spam) with P(ham) · P(words | ham). The denominator P(words) is the same for both classes, so it can be ignored. The naive assumption treats words as independent given the class, so P(words | class) becomes a product of per-word probabilities. Real systems add the logarithms of these probabilities, so that a product of many small numbers does not underflow. They also use Laplace smoothing, adding 1 to every word count, so that a word never seen in one class does not force that class's probability to 0.

from collections import Counter
from math import log, exp

def bayes(prior, sensitivity, false_positive_rate):
    evidence = sensitivity * prior + false_positive_rate * (1 - prior)   # law of total probability
    return sensitivity * prior / evidence

p = bayes(0.01, 0.99, 0.05)
print(round(p, 4), round(bayes(p, 0.99, 0.05), 4))   # 0.1667 0.7984   yesterday's posterior is today's prior

train = [("win money now", "spam"), ("cheap money offer", "spam"), ("win a cheap prize", "spam"),
         ("meeting moved to noon", "ham"), ("lunch money for friday", "ham"), ("project meeting notes", "ham")]
labels = Counter(label for _, label in train)
words = {c: Counter() for c in labels}
for text, label in train:
    words[label].update(text.split())
vocab = {w for c in words for w in words[c]}

def log_score(text, c):                  # log P(c) + Σ log P(word | c); the shared evidence term is omitted
    total = sum(words[c].values())
    score = log(labels[c] / len(train))
    for w in text.split():
        score += log((words[c][w] + 1) / (total + len(vocab)))   # Laplace smoothing
    return score

for msg in ["cheap money", "meeting at noon"]:
    s, h = log_score(msg, "spam"), log_score(msg, "ham")
    print(msg, "->", round(1 / (1 + exp(h - s)), 3))          # normalized P(spam | words)
# cheap money -> 0.829
# meeting at noon -> 0.157

Key takeaways

• P(B | A) = P(A | B)P(B) / P(A): posterior is proportional to likelihood times prior.
• Compute the denominator with the law of total probability, summing over every explanation.
• A positive result for a rare condition can still leave the condition unlikely. Always account for the base rate.
• In odds form, posterior odds = prior odds × likelihood ratio, which makes repeated independent updates easy.

Practice

1. (Beginner) 2% of transactions are fraudulent. A detector flags 95% of fraudulent transactions and 3% of legitimate ones. What fraction of flagged transactions are actually fraudulent?
2. (Intermediate) A bag contains three coins: a fair coin, a two-headed coin, and a coin that lands heads with probability 3/4. You pick a coin at random, flip it, and get heads. What is the probability that you picked the two-headed coin? What is it if a second flip is also heads?
3. (Interview) An alert fires for 90% of deploys that cause a real regression, and for 2% of healthy deploys. One deploy in 200 causes a regression. The on-call team wants to roll back automatically on every alert. What fraction of rollbacks would be unnecessary, and what would you change? Hint: compute the posterior first. Then consider requiring two independent signals and recompute with the odds form.

Birthday Problem

Why. In a group of only 23 people, the probability that two share a birthday is above 50%. That sounds wrong because people compare 23 with 365. What matters is the number of pairs, which grows quadratically: 23 people form 253 pairs. The same mathematics determines when hash values collide, how long cryptographic hashes must be (a "birthday attack" on an n-bit hash takes about 2^{n/2} attempts, not 2ⁿ), and how many random bits an ID needs.

What. Put n items into d equally likely slots, independently. The slots are days of the year, hash buckets, or random IDs. Count the outcomes with no collision as falling factorials (5.2):

                         d · (d − 1) ⋯ (d − n + 1)     n−1
  P(no collision)  =  ───────────────────────────  =   ∏  (1 − k/d)
                                   dⁿ                  k=0

Approximation. Using 1 − x ≤ e^{−x} (which is very accurate when x is small):

  P(no collision) ≈ e^{−n(n−1)/(2d)}          P(collision) ≈ 1 − e^{−n²/(2d)}

  50% threshold:  n ≈ √(2d · ln 2) ≈ 1.18 √d
  union bound:    P(collision) ≤ C(n, 2)/d      (useful when the probability is small)

So collisions become likely once n reaches about √d. For d = 365 the threshold is n ≈ 22.5, so 23 people suffice.

n people	10	20	23	30	50	70
P(shared birthday)	0.117	0.411	0.507	0.706	0.970	0.999

Where it matters.

Space of d values	≈ 50% collision after
365 birthdays	23
32-bit random IDs	about 77 000
64-bit hashes	about 5 × 10⁹
122 random bits (UUID v4)	about 2.7 × 10¹⁸

A different question. "Someone shares my birthday" compares one fixed day with the other n people, not all pairs with each other. That probability is 1 − (364/365)ⁿ, and it does not reach 50% until n = 253.

from math import prod, exp, log, sqrt, comb
import random

def p_collision(n, d):
    return 1 - prod(1 - k / d for k in range(n))

def p_approx(n, d):
    return 1 - exp(-n * (n - 1) / (2 * d))

print([round(p_collision(n, 365), 3) for n in (10, 20, 23, 30, 50, 70)])
# [0.117, 0.411, 0.507, 0.706, 0.97, 0.999]
print(round(p_approx(23, 365), 3), round(comb(23, 2) / 365, 3))   # 0.5 0.693   approximation, union bound
print(round(sqrt(2 * 365 * log(2)), 1), round(sqrt(2 * 2**32 * log(2))))   # 22.5 77163

rng = random.Random(0)
trials = 100_000
hits = sum(len({rng.randrange(365) for _ in range(23)}) < 23 for _ in range(trials))
print(round(hits / trials, 2))                                    # 0.51   simulation

n = 1
while 1 - (364 / 365) ** n < 0.5:                                  # someone shares MY birthday
    n += 1
print(n)                                                           # 253

Key takeaways

• With n items in d slots, P(no collision) = ∏_{k<n}(1 − k/d) ≈ e^{−n²/(2d)}.
• Collisions become likely at about √d items, because the number of pairs grows like n²/2.
• An n-bit hash withstands only about 2^{n/2} attempts to find a collision. Choose ID and hash sizes with this in mind.
• "Any two people share a birthday" and "someone shares my birthday" are very different events.

Practice

1. (Beginner) Assume birth months are uniform. What is the probability that at least two of 5 people share a birth month?
2. (Intermediate) A service tags each request with a random 32-bit ID. Roughly how many requests can it handle before a collision has a 1% probability? Hint: solve n²/(2d) ≈ 0.01.
3. (Interview) You generate random short-URL tokens from 62 characters (letters and digits). How many characters do you need so that 10⁹ links have a collision probability below 10⁻⁶? Hint: use the union bound C(n, 2)/d ≤ 10⁻⁶ to find the minimum d, then convert to a number of base-62 characters.

Monty Hall Problem

Why. This puzzle is famous because almost everyone's first answer is wrong. When Marilyn vos Savant published the correct answer in 1990, thousands of readers wrote in to say she was mistaken, including many with PhDs. Its lesson is general: what someone chooses to reveal is itself information. A probability model has to include the rules that produced the evidence, not just the evidence.

What. A car is behind one of three doors and goats are behind the other two. You pick a door. The host, who knows where the car is, then opens a different door that hides a goat. If both other doors hide goats, the host picks one of them at random. Then the host offers you the chance to switch to the remaining closed door. Switching wins with probability 2/3.

Argument 1: what switching does. If you stay, you win exactly when your first pick was right, which has probability 1/3. If you switch, you win exactly when your first pick was wrong, because the host has then removed the only other goat. That has probability 2/3.

Argument 2: Bayes' theorem. You pick door 1 and the host opens door 3. Let Cᵢ be the event that the car is behind door i.

  P(host opens 3 | C₁) = 1/2     (host chooses at random between doors 2 and 3)
  P(host opens 3 | C₂) = 1       (host must avoid your door and the car)
  P(host opens 3 | C₃) = 0       (host never reveals the car)

                          (1/2)(1/3)                     1/6
  P(C₁ | opens 3) = ──────────────────────────── = ───── = 1/3,     P(C₂ | opens 3) = (1/3)/(1/2) = 2/3
                      (1/2)(1/3) + 1(1/3) + 0           1/2

  car behind     host opens     probability    stay wins?   switch wins?
  door 1 ──┬──── door 2         1/3 · 1/2      yes          no
           └──── door 3         1/3 · 1/2      yes          no
  door 2 ─────── door 3         1/3 · 1        no           yes
  door 3 ─────── door 2         1/3 · 1        no           yes
                                               total 1/3    total 2/3

The host's rules decide the answer. Change the host's behavior and the same scene, "a goat behind door 3", means something else:

| Host's rules | P(switch wins | goat shown) | |---|---| | Always opens a goat door, choosing at random if there are two (standard) | 2/3 | | Opens one of the other two doors at random, and it happens to show a goat ("Monty Fall") | 1/2 | | Offers a switch only when your first pick is the car ("Monty Hell") | 0 | | n doors, and the host opens n − 2 goat doors | (n − 1)/n |

from fractions import Fraction as F
import random

def monty(host):
    """Exact P(stay wins), P(switch wins) given that a goat was revealed. The player picks door 0."""
    stay = switch = goat_shown = F(0)
    for car in range(3):
        options = [d for d in (1, 2) if d != car] if host == "knows" else [1, 2]
        for opened in options:
            p = F(1, 3) * F(1, len(options))
            if opened == car:
                continue                          # car revealed: excluded by conditioning on "goat shown"
            other = 3 - opened                    # the remaining closed door among {1, 2}
            goat_shown += p
            stay += p * (car == 0)
            switch += p * (car == other)
    return str(stay / goat_shown), str(switch / goat_shown)

print(monty("knows"))       # ('1/3', '2/3')   standard rules
print(monty("random"))      # ('1/2', '1/2')   Monty Fall

rng = random.Random(5)
wins = 0
for _ in range(100_000):
    car, pick = rng.randrange(3), rng.randrange(3)
    opened = rng.choice([d for d in range(3) if d not in (pick, car)])   # a goat door the player did not pick
    final = next(d for d in range(3) if d not in (pick, opened))         # switch
    wins += final == car
print(round(wins / 100_000, 2))                   # 0.67

Key takeaways

• Under the standard rules, switching wins with probability 2/3, because it wins exactly when the first pick was wrong.
• The host's choice is not random: it depends on where the car is, and that dependence is the information you gain.
• Bayes' theorem settles the question once the likelihoods P(host opens door j | car at i) are written down.
• The same observation can mean different things under different rules. Model how the evidence was produced.

Practice

1. (Beginner) There are four doors. After your pick, the host opens one goat door among the other three, and you switch to one of the two remaining closed doors at random. What is your probability of winning? Compare it with staying.
2. (Intermediate) Use Bayes' theorem to show that in "Monty Fall", where the host opens one of the two other doors at random and happens to reveal a goat, switching wins with probability 1/2.
3. (Interview) A colleague's Monty Hall simulation reports that switching wins 50% of the time. Name two bugs that could produce this result. Hint: check whether the simulation picks the host's door without looking at where the car is and then discards the runs where the car is revealed. Also check whether "switch" chooses at random between the two unopened doors.

<a id="73-random-variables"></a>

7.3 Random Variables

This subchapter moves from events to numbers. A random variable attaches a numerical value to every outcome, and its distribution, average, and spread are the quantities that algorithm analysis actually uses.

Discrete Random Variables

Why. Usually we care about a number rather than a raw outcome: the number of heads, the number of comparisons a sort makes, the length of a hash chain, the profit of a bet. Random variables let us do arithmetic with uncertain quantities: add them, average them, and bound them.

What. A random variable is a function X: Ω → ℝ (2.2). Despite the name, it is neither random nor a variable. The randomness lies in which outcome ω occurs, and X simply reports the number X(ω). X is discrete if it takes finitely or countably many values.

  Ω (3 tosses)        X = number of heads
  HHH  ─────────────►  3
  HHT, HTH, THH ─────►  2
  HTT, THT, TTH ─────►  1
  TTT  ─────────────►  0

• The notation {X = x} means the event {ω ∈ Ω : X(ω) = x}, the preimage of x. Likewise {X ≤ x}, {X ∈ A}, and so on.
• Arithmetic on random variables is done outcome by outcome: (X + Y)(ω) = X(ω) + Y(ω). Any function g(X), such as X², max(X, Y), or [X > 3], is again a random variable.
• Random variables X and Y are independent if P(X = x and Y = y) = P(X = x) · P(Y = y) for all values x and y. This means every event about X is independent of every event about Y.

Example: two dice. On Ω = {1, …, 6}², define X₁(a, b) = a, X₂(a, b) = b, and S = X₁ + X₂. The event {S = 4} is {(1, 3), (2, 2), (3, 1)}. X₁ and X₂ are independent. X₁ and S are not: P(X₁ = 1, S = 12) = 0, while P(X₁ = 1) · P(S = 12) = 1/6 · 1/36 > 0.

Quantity in an algorithm	Sample space	Random variable
Comparisons made by randomized quicksort	All sequences of pivot choices	The number of comparisons in that run
Longest chain in a hash table	All mⁿ assignments of keys to buckets	The largest bucket size
Retries until success	Sequences of attempt results	The index of the first success

from itertools import product
from fractions import Fraction as F

omega = list(product(range(1, 7), repeat=2))           # equally likely outcomes
X1 = lambda w: w[0]                                      # a random variable is a function on outcomes
X2 = lambda w: w[1]
S  = lambda w: X1(w) + X2(w)                             # the sum of two random variables

def P(event): return F(sum(1 for w in omega if event(w)), len(omega))

print(sorted(w for w in omega if S(w) == 4))             # [(1, 3), (2, 2), (3, 1)]   the event {S = 4}

def independent(X, Y):
    xs, ys = {X(w) for w in omega}, {Y(w) for w in omega}
    return all(P(lambda w: X(w) == x and Y(w) == y) == P(lambda w: X(w) == x) * P(lambda w: Y(w) == y)
               for x in xs for y in ys)

print(independent(X1, X2), independent(X1, S))           # True False

Key takeaways

• A random variable is a function from outcomes to numbers. The randomness comes from which outcome occurs.
• {X = x} is an event, namely the preimage of x. Sums and functions of random variables are random variables.
• Random variables are independent when the joint probabilities factor for every pair of values, not just for some.

Practice

1. (Beginner) Toss three coins. Let X be the number of heads and Y the length of the longest run of equal results (for HHT, Y = 2). Tabulate X and Y for all 8 outcomes.
2. (Intermediate) Two dice are rolled, with X₁ the first die and S the sum. Show that the events {X₁ = 1} and {S = 7} are independent, but the random variables X₁ and S are not.
3. (Interview) Define precisely, as a function on a sample space, the random variable "length of the longest chain" in a hash table with chaining that holds n keys in m buckets. What assumption about the hash function makes all outcomes equally likely? Hint: an outcome is an assignment of keys to buckets. The "simple uniform hashing" assumption is that each key goes to a uniformly random bucket, independently of the other keys.

Probability Mass Functions

Why. Once a random variable is defined, the full sample space is often more detail than you need. To answer any question about X alone, it is enough to know how much probability sits on each of its values. This summary is the probability mass function (PMF). It is also how distributions are stored in code and how random values are generated.

What.

• The probability mass function of X is p_X(x) = P(X = x). It satisfies p_X(x) ≥ 0 and Σ_x p_X(x) = 1, and P(X ∈ A) = Σ_{x ∈ A} p_X(x).
• The cumulative distribution function (CDF) is F_X(x) = P(X ≤ x). For a discrete X, it is a step function that rises from 0 to 1, jumping by p_X(x) at each value x.
• The joint PMF of two random variables is p_{X,Y}(x, y) = P(X = x, Y = y). Summing over one variable gives the marginal PMF of the other: p_X(x) = Σ_y p_{X,Y}(x, y).
• The PMF of a function of X is found by collecting values: p_{g(X)}(y) = Σ_{x : g(x) = y} p_X(x).

Example 1: the sum of two dice. p_S(s) = (6 − |s − 7|)/36 for s = 2, …, 12:

   s   p_S(s)
   2   1/36   █
   3   2/36   ██
   4   3/36   ███
   5   4/36   ████
   6   5/36   █████
   7   6/36   ██████
   8   5/36   █████
   9   4/36   ████
  10   3/36   ███
  11   2/36   ██
  12   1/36   █

Example 2: a joint PMF. Toss three coins. Let X be the number of heads, and let Y = 1 if the first toss is heads and Y = 0 otherwise:

	X = 0	X = 1	X = 2	X = 3	p_Y
Y = 0	1/8	2/8	1/8	0	1/2
Y = 1	0	1/8	2/8	1/8	1/2
p_X	1/8	3/8	3/8	1/8	1

The margins are the row and column sums. X and Y are not independent: the cell (X = 0, Y = 1) is 0, but the product of its margins is 1/16.

Sampling from a PMF (inverse transform). To draw a value with PMF p, list the values in order, compute the CDF, draw u uniformly from [0, 1), and return the first value whose CDF exceeds u. Value x is returned exactly when u falls in an interval of length p(x). With the CDF precomputed, a binary search makes each draw O(log n). This is how weighted load balancers and random.choices work.

from fractions import Fraction as F
from collections import defaultdict
from itertools import accumulate
import bisect, random

die = {k: F(1, 6) for k in range(1, 7)}

def pmf_sum(p, q):                        # PMF of X + Y for independent X ~ p, Y ~ q: a convolution (6.3)
    r = defaultdict(F)
    for x, px in p.items():
        for y, qy in q.items():
            r[x + y] += px * qy
    return dict(r)

S = pmf_sum(die, die)
print([str(S[s]) for s in range(2, 13)])
# ['1/36', '1/18', '1/12', '1/9', '5/36', '1/6', '5/36', '1/9', '1/12', '1/18', '1/36']
print(sum(S.values()))                    # 1

values = sorted(S)
cdf = list(accumulate(S[v] for v in values))      # cdf[i] = P(S ≤ values[i])
print(cdf[values.index(7)])                        # 7/12

def sample(rng):                          # inverse transform: first value whose CDF exceeds u
    return values[bisect.bisect_right(cdf, rng.random())]

rng = random.Random(1)
draws = [sample(rng) for _ in range(100_000)]
print(round(draws.count(7) / len(draws), 2))       # 0.17   close to 1/6

Key takeaways

• The PMF p_X(x) = P(X = x) describes everything about X alone. It is nonnegative and sums to 1.
• The CDF F_X(x) = P(X ≤ x) is a nondecreasing step function from 0 to 1.
• Marginal PMFs are the row and column sums of a joint PMF. X and Y are independent iff every cell equals the product of its margins.
• The PMF of the sum of independent random variables is a convolution. Inverse-transform sampling turns a CDF into a random generator.

Practice

1. (Beginner) Find the PMF and CDF of M, the maximum of two dice. Check that P(M = m) = (2m − 1)/36. Hint: P(M ≤ m) = (m/6)².
2. (Intermediate) In the joint table above, find the PMF of X − Y, the number of heads among the last two tosses, by grouping cells. Compare the result with the PMF of the number of heads in two tosses.
3. (Interview) Implement weighted random selection over n servers with integer weights, so that each pick takes O(log n) time. How would you support weight updates, and can you get O(1) per pick? Hint: prefix sums with binary search. For updates, use a Fenwick tree. For O(1) picks, look up the alias method.

Expected Value

Why. A whole distribution is a lot to carry around. Often one number is enough: the long-run average. Average-case running time, expected cost, the fair price of a bet, and average latency are all expected values. The expected value is the single most used summary of a random variable.

What. The expected value (mean) of a discrete random variable X is

  E[X] = Σ_x x · P(X = x) = Σ_{ω ∈ Ω} X(ω) · p(ω)

The first sum groups outcomes by value, and the second sums over outcomes. They are equal, and the second is often easier to use. Analogy: if each value x is a weight of mass p_X(x) placed on a number line, E[X] is the balance point.

• E[X] need not be a value that X can take: a die has E[X] = 3.5.
• Law of the unconscious statistician: E[g(X)] = Σ_x g(x) · p_X(x). You do not need the PMF of g(X).
• Conditional expectation and the law of total expectation: if B₁, …, Bₙ partition Ω, then E[X] = Σ E[X | Bᵢ] · P(Bᵢ). This is the law of total probability (7.2) for averages.
• An expected value can be infinite or undefined when Ω is infinite. See the St. Petersburg game below.

Example 1: a bet. You pay 1 to roll a die and receive 4 if it shows a 6. Your net gain is 3 with probability 1/6 and −1 with probability 5/6, so the expected gain is 3/6 − 5/6 = −1/3 per game.

Example 2: linear search. A key is equally likely to be at any of the n positions of an array. Finding it at position i takes i comparisons, so E[comparisons] = (1 + 2 + ⋯ + n)/n = (n + 1)/2.

Example 3: a cache (total expectation). Reads hit the cache 95% of the time, with latency 2 ms, and otherwise go to disk, with latency 100 ms. The expected latency is 0.95 · 2 + 0.05 · 100 = 6.9 ms. The rare misses dominate the average.

Example 4: waiting for a six (conditioning on the first step). Let W be the number of rolls until the first 6. The first roll is a 6 with probability 1/6. Otherwise the process starts over, with one roll already used:

  E[W] = 1 + (5/6)·E[W]      ⇒      E[W] = 6

The same argument with success probability p gives E[W] = 1/p. This waiting time has the geometric distribution (7.4).

Example 5: an infinite expectation (St. Petersburg). A fair coin is tossed until the first head. If that happens on toss k, the game pays 2ᵏ. Each k contributes 2ᵏ · 2⁻ᵏ = 1 to the sum, so E[payout] = 1 + 1 + 1 + ⋯ = ∞. Yet no reasonable person would pay a large amount to play. If the bank can pay out at most 2³⁰, the expected payout is only 31. The mean is a poor summary of distributions with very heavy tails.

from fractions import Fraction as F

def E(pmf, g=lambda x: x):                 # E[g(X)] = Σ g(x)·p(x)
    return sum(g(x) * p for x, p in pmf.items())

die = {k: F(1, 6) for k in range(1, 7)}
print(E(die), E(die, lambda x: x * x))     # 7/2 91/6

bet = {3: F(1, 6), -1: F(5, 6)}            # net gain: +3 on a six, −1 otherwise
print(E(bet))                              # -1/3

n = 100
linear_search = {i: F(1, n) for i in range(1, n + 1)}
print(E(linear_search))                    # 101/2

print(round(0.95 * 2 + 0.05 * 100, 1))     # 6.9   total expectation over hit / miss

p = F(1, 6)                                # E[W] = Σ k·(1 − p)^(k−1)·p, truncated after 400 terms
print(round(float(sum(k * (1 - p) ** (k - 1) * p for k in range(1, 400))), 6))   # 6.0

cap = 2 ** 30                              # St. Petersburg with a bank that pays at most 2^30
print(round(float(sum(min(2 ** k, cap) * F(1, 2 ** k) for k in range(1, 200))), 6))   # 31.0   30 terms of 1, plus a tail worth 1

Key takeaways

• E[X] = Σ x · P(X = x), the probability-weighted average of the values, or the balance point of the distribution.
• E[g(X)] = Σ g(x)p_X(x), so there is no need to find the PMF of g(X) first.
• Law of total expectation: E[X] = Σ E[X | Bᵢ]P(Bᵢ). Conditioning on the first step gives equations such as E = 1 + (1 − p)E.
• Waiting for an event of probability p takes 1/p trials on average.
• A mean can be infinite or misleading when the tail is heavy.

Practice

1. (Beginner) Compute E[X] and E[X²] for a fair die, and check that E[X²] ≠ (E[X])².
2. (Intermediate) A request succeeds on each attempt with probability 0.8, independently, and the client makes at most 3 attempts. Find the expected number of attempts and the probability that the request ultimately fails.
3. (Interview) A teammate reports a server's mean throughput as 1/E[latency]. Is E[1/X] = 1/E[X] in general? Hint: try X uniform on {1, 3}. For positive X, which side is always at least as large? (This is a special case of Jensen's inequality.)

Linearity of Expectation

Why. This is the most useful tool in the probabilistic analysis of algorithms. The expected value of a sum is the sum of the expected values, always, even when the terms are dependent. To find the expectation of a complicated count, split it into simple pieces, find the expectation of each piece, and add. You never need the distribution of the total.

What. For any random variables X and Y on the same space, and any constants a, b, c:

  E[aX + bY + c] = a·E[X] + b·E[Y] + c           and     E[X₁ + ⋯ + Xₙ] = E[X₁] + ⋯ + E[Xₙ]

Proof. E[aX + bY + c] = Σ_ω (aX(ω) + bY(ω) + c)p(ω) = aΣ_ω X(ω)p(ω) + bΣ_ω Y(ω)p(ω) + cΣ_ω p(ω). The proof uses only the rearrangement of a sum, so independence plays no role.

Statement	True in general?
E[X + Y] = E[X] + E[Y]	Yes, always
E[cX] = c·E[X]	Yes, always
E[XY] = E[X]·E[Y]	Only when X and Y are uncorrelated (for example, independent)
E[g(X)] = g(E[X])	Only when g is linear

Example 1: 100 dice. The distribution of the sum of 100 dice is complicated, but its mean is not: 100 · 3.5 = 350.

Example 2: fixed points (the hat-check problem). n hats are returned to n guests in a uniformly random order. Let Xᵢ = 1 if guest i gets their own hat back and Xᵢ = 0 otherwise. Then E[Xᵢ] = 1/n, so

  E[number of guests with their own hat] = E[X₁ + ⋯ + Xₙ] = n · (1/n) = 1       for every n

The Xᵢ are dependent (if n − 1 guests have their own hats, so does the last one), and linearity does not care. Compare this with the derangement count in 5.5, which needed inclusion–exclusion to find the full probability.

Example 3: the coupon collector. Each cereal box contains one of n coupon types, uniformly at random. How many boxes do you need to buy to collect all n types? Split the process into phases. Phase i starts when you hold i − 1 distinct coupons and ends when you get a new one. During phase i, each box is new with probability (n − i + 1)/n, so by Example 4 of Expected Value, the phase lasts n/(n − i + 1) boxes on average:

  E[boxes] = n/n + n/(n − 1) + ⋯ + n/1 = n · (1 + 1/2 + ⋯ + 1/n) = n·Hₙ ≈ n ln n

  n = 6 (faces of a die):  6 · H₆ = 14.7 rolls to see every face

The same calculation tells you how many random requests it takes to touch every server, how many random tests it takes to exercise every code path, and how many random samples it takes to see every bucket.

from itertools import permutations
from fractions import Fraction as F
import random

perms = list(permutations(range(5)))                               # all 120 permutations
print(F(sum(sum(p[i] == i for i in range(5)) for p in perms), len(perms)))   # 1

def coupon_expected(n):
    return n * sum(F(1, k) for k in range(1, n + 1))               # n·H_n

def coupon_sim(n, rng):
    seen, draws = set(), 0
    while len(seen) < n:
        seen.add(rng.randrange(n))
        draws += 1
    return draws

rng = random.Random(11)
print(float(coupon_expected(6)), round(float(coupon_expected(50)), 1))         # 14.7 225.0
print(round(sum(coupon_sim(6, rng) for _ in range(100_000)) / 100_000, 1))     # 14.7

Key takeaways

• E[Σ Xᵢ] = Σ E[Xᵢ] holds even when the Xᵢ are dependent.
• Strategy: write the quantity as a sum of simple pieces, find the expectation of each, and add.
• E[XY] = E[X]E[Y] requires independence (or at least zero correlation), and E[g(X)] ≠ g(E[X]) for nonlinear g.
• A random permutation has 1 fixed point on average. The coupon collector needs n·Hₙ ≈ n ln n draws on average.

Practice

1. (Beginner) A biased coin lands heads with probability 0.3. What is the expected number of heads in 10 tosses? What is the expected total of 3 dice plus twice a fourth die?
2. (Intermediate) A fair coin is tossed n times. Show that the expected number of positions where HH begins (overlapping occurrences count) is (n − 1)/4, even though occurrences overlap and are dependent.
3. (Interview) The loop best = -inf; for x in arr: if x > best: best = x runs over a uniformly random permutation of n distinct numbers. How many times, on average, does the assignment best = x execute? Hint: the assignment runs at position i exactly when arr[i] is the largest of the first i elements. By symmetry, that has probability 1/i. Add up these probabilities.

Indicator Random Variables

Why. Linearity of expectation is most powerful when combined with a simple trick: write the quantity you want as a count, and write the count as a sum of 0/1 variables. Each 0/1 variable has an expectation equal to a probability, which is usually easy to find. This recipe computes the average number of collisions, inversions, empty buckets, fixed points, and many other quantities, and it is the standard tool for analyzing randomized algorithms.

What. The indicator of an event A is the random variable

  I_A(ω) = 1 if ω ∈ A,  0 otherwise           and          E[I_A] = 1·P(A) + 0·P(Ā) = P(A)

Useful identities: I_{A ∩ B} = I_A · I_B, I_Ā = 1 − I_A, and I_A² = I_A.

How: the indicator recipe.

1. Write the quantity X as a count: X = the number of events among A₁, …, Aₘ that happen.
2. Write X = I_{A₁} + ⋯ + I_{Aₘ}.
3. By linearity, E[X] = P(A₁) + ⋯ + P(Aₘ). Find each probability, often by symmetry.

Example 1: hash tables. Hash n keys into m buckets uniformly and independently.

  empty buckets:     I_j = [bucket j is empty],   P = (1 − 1/m)ⁿ        E = m(1 − 1/m)ⁿ ≈ m·e^{−n/m}
  colliding pairs:   I_{ab} = [keys a, b collide], P = 1/m               E = C(n, 2)/m
  chain length:      key a's bucket contains a, plus each of the n − 1 other keys with prob 1/m
                                                                          E = 1 + (n − 1)/m

With n = m, about 1/e ≈ 37% of the buckets are empty, and the chain holding any given key contains 2 − 1/m keys on average. The expected number of colliding pairs, C(n, 2)/m, is the birthday problem's union bound (7.2), now read as an expectation.

Example 2: insertion sort. An inversion is a pair of positions i < j with a[i] > a[j]. Insertion sort performs exactly one swap per inversion. In a random permutation, each of the C(n, 2) pairs is inverted with probability 1/2, by symmetry. So

  E[inversions] = C(n, 2) · 1/2 = n(n − 1)/4

On a random input, insertion sort therefore averages about n²/4 swaps, half the worst case and still Θ(n²).

from itertools import permutations
from fractions import Fraction as F
from math import comb
import random

n = 5
perms = list(permutations(range(n)))
inv = lambda p: sum(p[i] > p[j] for i in range(n) for j in range(i + 1, n))
print(F(sum(inv(p) for p in perms), len(perms)), F(n * (n - 1), 4))   # 5 5

keys, m, trials = 100, 100, 2000
rng = random.Random(2)
empty = pairs = 0
for _ in range(trials):
    counts = [0] * m
    for _ in range(keys):
        counts[rng.randrange(m)] += 1
    empty += counts.count(0)
    pairs += sum(c * (c - 1) // 2 for c in counts)                    # colliding pairs inside each bucket
print(round(empty / trials, 1), round(m * (1 - 1 / m) ** keys, 1))    # 36.6 36.6
print(round(pairs / trials, 1), comb(keys, 2) / m)                    # 49.4 49.5

Key takeaways

• E[I_A] = P(A), which turns expected counts into sums of probabilities.
• Recipe: count events, write the count as a sum of indicators, and add the probabilities.
• Symmetry often gives the probabilities directly: 1/2 for a pair to be inverted, 1/m for two keys to collide, 1/n for a fixed point.
• Hashing n keys into m buckets: C(n, 2)/m expected colliding pairs and m(1 − 1/m)ⁿ expected empty buckets.

Practice

1. (Beginner) An ascent in a permutation is a position i with a[i] < a[i + 1]. Show that a random permutation of n elements has (n − 1)/2 ascents on average.
2. (Intermediate) An m-sided die is rolled n times. Find the expected number of distinct values seen. Evaluate it for m = n = 100, and compare the result with 100(1 − 1/e).
3. (Interview) In a random graph on n vertices where each possible edge is present independently with probability p, what is the expected number of triangles? For which p does this expectation stay bounded as n grows? Hint: use one indicator per 3-subset of vertices. A triangle needs all three of its edges.

Variance and Standard Deviation

Why. The expected value hides the spread. Two services can both average 100 ms, one always taking 100 ms and the other usually fast but sometimes taking half a second. Users notice the difference. Variance measures how far a random variable typically strays from its mean. It is the input to Chebyshev's inequality (7.4), and it explains why averaging many measurements gives a more precise result.

What. Let μ = E[X].

  Var(X) = E[(X − μ)²] = E[X²] − μ²            σ_X = √Var(X)   (standard deviation, in the units of X)

The second form follows from expanding the square and applying linearity: E[X² − 2μX + μ²] = E[X²] − 2μ² + μ². It is usually easier to compute.

Property	Note
Var(X) ≥ 0, and Var(X) = 0 iff X is constant with probability 1
Var(aX + b) = a²·Var(X)	Shifting does not change spread. Scaling by a scales σ by `
Var(X + Y) = Var(X) + Var(Y) + 2·Cov(X, Y)	Cov(X, Y) = E[XY] − E[X]E[Y]
Var(X₁ + ⋯ + Xₙ) = Σ Var(Xᵢ) if the Xᵢ are pairwise independent	Unlike the expectation, variance adds only without correlation
Var(I_A) = p(1 − p) where p = P(A)	Because I_A² = I_A
Mean of n independent copies: Var(X̄) = σ²/n	Averaging shrinks σ by a factor of √n

Example 1: a die. E[X] = 7/2 and E[X²] = 91/6, so Var(X) = 91/6 − 49/4 = 35/12 ≈ 2.92 and σ ≈ 1.71. For the sum of two independent dice, Var(S) = 35/6.

Example 2: two services with the same mean.

	Latency distribution	Mean	Standard deviation
Service A	Always 100 ms	100 ms	0 ms
Service B	50 ms with prob 0.9, 550 ms with prob 0.1	100 ms	150 ms

Example 3: covariance matters. The number of fixed points of a random permutation is X = Σ Iᵢ. The indicators are dependent: P(Iᵢ = Iⱼ = 1) = 1/(n(n − 1)), which is larger than 1/n². So each pair has a small positive covariance:

  Var(X) = Σᵢ Var(Iᵢ) + Σ_{i≠j} Cov(Iᵢ, Iⱼ)
         = n·(1/n)(1 − 1/n) + n(n − 1)·(1/(n(n − 1)) − 1/n²)
         = (1 − 1/n) + (1/n) = 1                                 for every n ≥ 2

Why averaging works. If one benchmark run has standard deviation σ, the mean of n independent runs has standard deviation σ/√n. This is the standard error. Four times as many runs halve the uncertainty.

from fractions import Fraction as F
from itertools import permutations
from statistics import pstdev, mean
import random

def E(pmf, g=lambda x: x): return sum(g(x) * p for x, p in pmf.items())
def Var(pmf): return E(pmf, lambda x: x * x) - E(pmf) ** 2

die = {k: F(1, 6) for k in range(1, 7)}
print(Var(die), round(float(Var(die)) ** 0.5, 2))       # 35/12 1.71

B = {50: F(9, 10), 550: F(1, 10)}
print(E(B), Var(B), float(Var(B)) ** 0.5)              # 100 22500 150.0

n = 5
perms = list(permutations(range(n)))
fixed = [sum(p[i] == i for i in range(n)) for p in perms]
print(mean(fixed), pstdev(fixed) ** 2)                   # 1 1.0   variance of fixed points is 1

rng = random.Random(4)
for runs in (1, 4, 16, 64):                              # standard error of the mean of `runs` dice
    means = [mean(rng.randint(1, 6) for _ in range(runs)) for _ in range(20_000)]
    print(runs, round(pstdev(means), 2), round(1.708 / runs ** 0.5, 2))
# 1 1.7 1.71
# 4 0.85 0.85
# 16 0.43 0.43
# 64 0.21 0.21

Key takeaways

• Var(X) = E[(X − μ)²] = E[X²] − (E[X])². The standard deviation σ = √Var(X) has the same units as X.
• Var(aX + b) = a²Var(X).
• Variances of independent (or merely uncorrelated) variables add. In general, the covariance terms must be included.
• The average of n independent measurements has standard deviation σ/√n.

Practice

1. (Beginner) Find the variance of an indicator variable with P(A) = p. For which p is it largest?
2. (Intermediate) Show that Var(X + X) = 4Var(X), while Var(X + Y) = 2Var(X) when Y is an independent copy of X. Explain the difference using covariance.
3. (Interview) You benchmark a function 25 times and get a mean of 40 ms with a standard deviation of 10 ms. How precise is the mean? How many runs would you need to halve that uncertainty, and why might the runs not be independent in practice? Hint: standard error σ/√n. Think about caches warming up, CPU frequency scaling, and noisy neighbors.

<a id="74-distributions-and-bounds"></a>

7.4 Distributions and Bounds

This subchapter names the distributions that keep reappearing in computing and bounds how far a random variable can stray from its mean when its full distribution is unknown. It ends by putting randomness to work in proofs and in algorithms.

Bernoulli and Binomial Distributions

Why. Many processes consist of independent yes/no trials that repeat: packets are lost or delivered, a flaky test passes or fails, a coin lands heads or tails, an item is defective or fine. The binomial distribution counts the successes. It answers questions such as "how likely is it that 3 or more of these 200 packets are lost?"

What.

• A Bernoulli(p) random variable takes the value 1 ("success") with probability p and 0 with probability 1 − p. It is the indicator of a single event, so E[X] = p and Var(X) = p(1 − p).
• A Binomial(n, p) random variable counts the successes in n independent Bernoulli(p) trials:

  P(X = k) = C(n, k) · pᵏ · (1 − p)ⁿ⁻ᵏ          k = 0, 1, …, n

Where the formula comes from. Any particular sequence with k successes and n − k failures has probability pᵏ(1 − p)ⁿ⁻ᵏ, by independence. There are C(n, k) such sequences (5.2), one for each choice of the positions of the successes. The probabilities sum to (p + (1 − p))ⁿ = 1 by the binomial theorem (5.3), which is how the distribution got its name.

Mean and variance. Write X = X₁ + ⋯ + Xₙ, where Xᵢ is the indicator of success on trial i. By linearity, E[X] = np. Because the trials are independent, the variances add: Var(X) = np(1 − p). No sums of binomial coefficients are needed.

  Binomial(10, 0.3)
  k   P(X = k)
  0   0.028  █
  1   0.121  ██████
  2   0.233  ████████████
  3   0.267  █████████████        mode near np = 3
  4   0.200  ██████████
  5   0.103  █████
  6   0.037  ██
  7   0.009
  8+  0.002

With versus without replacement. The binomial distribution assumes independent trials, as when sampling with replacement. Sampling without replacement gives the hypergeometric distribution (7.1):

	Binomial(n, p)	Hypergeometric
Model	n independent trials with success probability p	n draws without replacement from N items, m of them successes
PMF	C(n, k)pᵏ(1 − p)ⁿ⁻ᵏ	C(m, k)C(N − m, n − k)/C(N, n)
Mean	np	n·m/N (the same, with p = m/N)
Variance	np(1 − p)	Smaller, by a factor of (N − n)/(N − 1)
When to use it	Independent events, or a population much larger than the sample	Small populations, as in quality control and card games

Example 1: packet loss. Each of 200 packets is lost independently with probability 0.01. P(no loss) = 0.99²⁰⁰ ≈ 0.134, and P(at most 2 lost) ≈ 0.677. The expected number of losses is 2, with a standard deviation of about 1.41.

Example 2: majority voting. Five independent classifiers are each correct with probability 0.7. A majority vote is correct when at least 3 of them are correct: Σ_{k=3}^{5} C(5, k)·0.7ᵏ·0.3⁵⁻ᵏ ≈ 0.837. Combining weak but independent voters improves accuracy. That is the idea behind ensemble methods and quorum systems, and it fails when the voters make correlated errors.

from math import comb, sqrt

def binom_pmf(n, p):
    return [comb(n, k) * p**k * (1 - p)**(n - k) for k in range(n + 1)]

pmf = binom_pmf(10, 0.3)
mean = sum(k * q for k, q in enumerate(pmf))
var = sum(k * k * q for k, q in enumerate(pmf)) - mean**2
print(round(mean, 6), round(var, 6))                   # 3.0 2.1   np and np(1 − p)

loss = binom_pmf(200, 0.01)
print(round(loss[0], 3), round(sum(loss[:3]), 3), round(sqrt(200 * 0.01 * 0.99), 2))   # 0.134 0.677 1.41

vote = binom_pmf(5, 0.7)
print(round(sum(vote[3:]), 3))                         # 0.837   majority of 5 voters correct

# Without vs with replacement: 20 items sampled from 1000, of which 50 are defective
print(round(comb(950, 20) / comb(1000, 20), 4), round(0.95**20, 4))   # 0.3549 0.3585

Key takeaways

• Bernoulli(p): a single 0/1 trial, with mean p and variance p(1 − p).
• Binomial(n, p): the number of successes in n independent trials, P(X = k) = C(n, k)pᵏ(1 − p)ⁿ⁻ᵏ.
• E[X] = np and Var(X) = np(1 − p), both obtained by writing X as a sum of indicators.
• Use the hypergeometric distribution for sampling without replacement. It is close to the binomial when the population is much larger than the sample.

Practice

1. (Beginner) A fair coin is tossed 10 times. Find P(exactly 5 heads) and P(at least 8 heads).
2. (Intermediate) A flaky test fails independently in 2% of runs, and CI runs it 50 times a day. Find the probability of at least one failure in a day, the expected number of failures, and the standard deviation.
3. (Interview) Each replica of a database is up independently with probability 0.9. A write needs a majority of replicas to be up. How many replicas (an odd number) do you need for writes to succeed with probability at least 0.999? Why is "independently" the weakest part of this calculation? Hint: evaluate the binomial tail P(X ≥ (n + 1)/2) for n = 3, 5, 7, ….

Geometric Distribution

Why. Many algorithms repeat an attempt until it succeeds: retrying a request, rejection sampling, probing a hash table, restarting a randomized algorithm. The number of attempts is random, and the geometric distribution describes it. It answers both "how many attempts on average?" and "how likely is it that we need more than k?"

What. Repeat independent Bernoulli(p) trials, and let X be the number of the trial on which the first success occurs:

  P(X = k) = (1 − p)^{k−1} · p          k = 1, 2, 3, …        (k − 1 failures, then a success)
  P(X > k) = (1 − p)ᵏ                                         (the first k trials all fail)
  E[X] = 1/p                    Var(X) = (1 − p)/p²

The probabilities form a geometric series (3.1) that sums to 1, which gives the distribution its name. The mean also follows from the tail-sum formula for random variables that take values in {0, 1, 2, …}: E[X] = Σ_{k ≥ 0} P(X > k) = Σ_{k ≥ 0} (1 − p)ᵏ = 1/p. 7.3 derived the same result by conditioning on the first step.

Watch the convention. Some textbooks and libraries define the geometric distribution as the number of failures before the first success. That variable takes values 0, 1, 2, … and has mean (1 − p)/p. Check which definition a library function uses before relying on it.

Memorylessness. For all s, t ≥ 0:

  P(X > s + t | X > s) = (1 − p)^{s+t} / (1 − p)ˢ = (1 − p)ᵗ = P(X > t)

Having already failed s times does not change the distribution of the remaining wait. The geometric distribution is the only distribution on the positive integers with this property. For example, a roulette wheel that has shown red ten times in a row is exactly as likely to show black on the next spin as on any other spin. Expecting otherwise is the gambler's fallacy.

Example 1: rejection sampling. You need a fair die roll, but your source produces random bits. Draw 3 bits to get a number from 0 to 7. If it is below 6, return it plus 1; otherwise draw again. Each attempt succeeds with probability 6/8, so the number of attempts is Geometric(3/4). On average that is 4/3 attempts, or 4 bits per die roll, and the probability of needing more than 5 attempts is (1/4)⁵ < 0.1%.

Example 2: earlier results revisited. The derangement sampler of 5.5 retries a shuffle until it has no fixed point, so it makes about e ≈ 2.72 attempts on average. Each phase of the coupon collector (7.3) is geometric, with success probability (n − i + 1)/n.

import random

def roll_from_bits(rng):
    attempts = 0
    while True:
        attempts += 1
        v = rng.getrandbits(3)            # uniform on 0..7
        if v < 6:                         # accepted with probability 6/8
            return v + 1, attempts

rng = random.Random(42)
results = [roll_from_bits(rng) for _ in range(200_000)]
attempts = [a for _, a in results]
print(round(sum(attempts) / len(attempts), 3))                     # 1.333   E[X] = 1/p = 4/3
print(round(sum(a > 2 for a in attempts) / len(attempts), 4), (1/4) ** 2)   # 0.0625 0.0625   P(X > 2)

# Memorylessness, checked on the simulated data: P(X > 3 | X > 1) vs P(X > 2)
over1 = [a for a in attempts if a > 1]
print(round(sum(a > 3 for a in over1) / len(over1), 3))            # 0.062   ≈ P(X > 2) = 0.0625

Key takeaways

• Geometric(p) is the trial on which the first success occurs: P(X = k) = (1 − p)^{k−1}p, and P(X > k) = (1 − p)ᵏ.
• E[X] = 1/p and Var(X) = (1 − p)/p².
• It is memoryless: past failures do not change the distribution of the remaining wait.
• Rejection sampling with acceptance probability p takes 1/p attempts on average, and the probability of a long run of rejections decays exponentially.

Practice

1. (Beginner) You roll a die until you get a 6. What is the probability that this takes more than 6 rolls? What is the expected number of rolls?
2. (Intermediate) Prove memorylessness from P(X > k) = (1 − p)ᵏ. Then show that a random variable on the positive integers with P(X > s + t) = P(X > s)P(X > t) for all s and t must be geometric. Hint: let q = P(X > 1) and use induction on t.
3. (Interview) Implement rand7(), which returns a uniform integer from 1 to 7, using only rand5(). What is the expected number of calls to rand5()? Hint: two calls give 25 equally likely outcomes. Keep 21 of them and reject 4. Each attempt is a Bernoulli(21/25) trial.

Markov's Inequality

Why. Often you know only the mean of a quantity: an algorithm's expected running time, a server's average queue length, the average size of a file. How likely is the quantity to be much larger than its mean? Markov's inequality gives an answer that uses nothing but the mean and the fact that the quantity is never negative. It turns an expected running time into a probability guarantee.

What. If X ≥ 0 and a > 0, then

  P(X ≥ a) ≤ E[X] / a          equivalently      P(X ≥ k·E[X]) ≤ 1/k

Proof. Compare X with the random variable a · I_{X ≥ a}. When X ≥ a, the indicator is 1 and a ≤ X. Otherwise the indicator is 0 and 0 ≤ X, because X is nonnegative. So a · I_{X ≥ a} ≤ X for every outcome, and taking expectations gives a · P(X ≥ a) ≤ E[X].

Analogy. If the average household income is 50 000 and no income is negative, then at most 10% of households can earn 500 000 or more. If more did, they alone would raise the average above 50 000.

Properties.

• Nonnegativity is essential. If X can be negative, large negative values can balance large positive ones, and the bound fails. For example, X = ±100 with equal probability has E[X] = 0, yet P(X ≥ 100) = 1/2.
• The bound is tight. Let X = a with probability μ/a and X = 0 otherwise. Then P(X ≥ a) = μ/a exactly. Knowing only the mean, you cannot do better.
• It is usually loose. For a particular distribution, the true probability can be far smaller than the bound. Using more information (the variance, or independence) gives sharper bounds: Chebyshev's inequality and Chernoff bounds.

Application: from Las Vegas to Monte Carlo. A randomized algorithm always returns a correct answer, and its expected running time is T. Stop it after 2T steps. By Markov's inequality, it has not finished with probability at most 1/2. If you restart it with fresh random bits up to k times, the probability that every run times out is at most 2⁻ᵏ. A guarantee about the average becomes a guarantee about a deadline.

Example. Quicksort with random pivots makes about 2n ln n comparisons on average (Randomized Algorithms below). Markov's inequality alone says that the probability of making 10 times that many is at most 10%. The true probability is astronomically smaller, but the Markov bound needs no analysis beyond the mean.

from math import comb

n, p, a = 100, 0.5, 75
pmf = [comb(n, k) * p**k * (1 - p)**(n - k) for k in range(n + 1)]
mu = n * p
exact = sum(pmf[a:])
print(f"P(X >= {a}): exact {exact:.1e}, Markov bound {mu / a:.3f}")   # P(X >= 75): exact 2.8e-07, Markov bound 0.667

# Tightness: X = a with probability mu/a, else 0
mu, a = 2.0, 8.0
tight = {a: mu / a, 0.0: 1 - mu / a}
print(sum(x * q for x, q in tight.items()), tight[a], mu / a)       # 2.0 0.25 0.25   bound attained exactly

Key takeaways

• For X ≥ 0: P(X ≥ a) ≤ E[X]/a, so the probability of exceeding k times the mean is at most 1/k.
• The proof is one line: a · I_{X ≥ a} ≤ X, then take expectations.
• It needs only the mean and nonnegativity. It is tight in general but often loose for a specific distribution.
• Combined with independent restarts, it turns an expected running time into a deadline that is missed only with exponentially small probability.

Practice

1. (Beginner) A service's average response time is 200 ms. Give an upper bound on the fraction of requests that take at least 2 s. What extra assumption did you need?
2. (Intermediate) Apply Markov's inequality to X = number of heads in 100 fair tosses to bound P(X ≥ 60). Then apply it to Y = (X − 50)² to bound the same probability, and compare the two bounds. Hint: E[Y] = Var(X) = 25. The second bound is Chebyshev's inequality in disguise.
3. (Interview) A Las Vegas algorithm has an expected running time of 1 second. Design a procedure that always finishes within 20 seconds and fails with probability at most 10⁻³. Hint: choose a timeout per run and a number of independent runs, then apply Markov's inequality to each run.

Chebyshev's Inequality

Why. Markov's inequality uses only the mean, and its bounds are often weak. If you also know the variance, you can bound the probability of deviating from the mean in either direction, and the bound improves as the spread shrinks. Chebyshev's inequality is the reason that averages of many samples are reliable (the weak law of large numbers). It is also the simplest way to decide how many samples an estimate needs.

What. For any random variable X with mean μ and variance σ², and any a > 0:

  P(|X − μ| ≥ a) ≤ σ² / a²          equivalently      P(|X − μ| ≥ kσ) ≤ 1/k²

Proof. (X − μ)² is nonnegative and has mean σ². Apply Markov's inequality to it with threshold a²: P((X − μ)² ≥ a²) ≤ σ²/a², and (X − μ)² ≥ a² is the same event as |X − μ| ≥ a.

| k (in standard deviations) | Chebyshev: P(|X − μ| ≥ kσ) ≤ | Example: Binomial(100, 0.5), actual probability | |---|---|---| | 2 | 25% | about 5.7% | | 3 | 11.1% | about 0.35% | | 5 | 4% | about 6 × 10⁻⁷ |

The bound holds for every distribution, so it is conservative for well-behaved ones such as the binomial. Its value is that it needs nothing beyond the variance.

The weak law of large numbers. Let X̄ be the average of n independent samples, each with mean μ and variance σ². By 7.3, Var(X̄) = σ²/n, so

  P(|X̄ − μ| ≥ ε) ≤ σ² / (n ε²)  →  0     as n → ∞

The sample average concentrates around the true mean. This justifies estimating probabilities by simulation, as many code examples in this chapter do.

Application: sample size. To estimate a proportion p, such as the fraction of users who click a button, sample n users and use the observed fraction. Each sample is a Bernoulli(p) trial with variance p(1 − p) ≤ 1/4. To get within ±ε of p, except with probability δ, Chebyshev's inequality requires

  1 / (4 n ε²) ≤ δ      ⇒      n ≥ 1 / (4 ε² δ)          ε = 0.01, δ = 0.05:  n ≥ 50 000

Beyond Chebyshev. For sums of independent bounded variables, Chernoff and Hoeffding bounds give tails that decay exponentially. For example, Hoeffding's inequality gives P(|X̄ − μ| ≥ ε) ≤ 2e^{−2nε²}, and the same accuracy target then needs only about 18 445 samples. Use Markov's inequality when you know only the mean, Chebyshev's when you know the variance, and Chernoff when you have independence.

from math import comb, exp, log, ceil
import random

n, p = 100, 0.5
pmf = [comb(n, k) * p**k * (1 - p)**(n - k) for k in range(n + 1)]
mu, var, a = n * p, n * p * (1 - p), 75
print(f"exact {sum(pmf[a:]):.1e}")                                   # exact 2.8e-07
print(f"Markov {mu / a:.3f}  Chebyshev {var / (a - mu)**2:.3f}  Hoeffding {exp(-2 * (a - mu)**2 / n):.1e}")
# Markov 0.667  Chebyshev 0.040  Hoeffding 3.7e-06

for k in (2, 3, 5):                                                  # P(|X − 50| ≥ 5k) vs 1/k²
    actual = sum(q for x, q in enumerate(pmf) if abs(x - mu) >= 5 * k)
    print(k, round(1 / k**2, 3), f"{actual:.1e}")
# 2 0.25 5.7e-02
# 3 0.111 3.5e-03
# 5 0.04 5.6e-07

eps, delta = 0.01, 0.05
print(ceil(1 / (4 * eps**2 * delta)), ceil(log(2 / delta) / (2 * eps**2)))   # 50000 18445

rng = random.Random(8)                     # how often is a 50 000-sample estimate of p = 0.3 off by ≥ 0.01?
misses = sum(abs(rng.binomialvariate(50_000, 0.3) / 50_000 - 0.3) >= 0.01 for _ in range(2000))
print(misses)                              # 0   far fewer than Chebyshev's allowance of 5%

Key takeaways

• P(|X − μ| ≥ kσ) ≤ 1/k² for every distribution with finite variance. It is Markov's inequality applied to (X − μ)².
• It gives two-sided bounds that shrink as the variance shrinks.
• The average of n independent samples concentrates: P(|X̄ − μ| ≥ ε) ≤ σ²/(nε²).
• To estimate a proportion within ±ε except with probability δ, n = 1/(4ε²δ) samples suffice. Chernoff and Hoeffding bounds need far fewer, if the samples are independent.

Practice

1. (Beginner) X has mean 100 and standard deviation 10. Bound P(X ≤ 70 or X ≥ 130).
2. (Intermediate) A fair coin is tossed 10 000 times. Use Chebyshev's inequality to bound the probability that the number of heads differs from 5000 by 200 or more. Then find how many tosses guarantee that the observed fraction of heads is within 0.01 of 1/2 with probability at least 99%.
3. (Interview) How many users must you sample to estimate a conversion rate within ±2 percentage points with 95% confidence? Give Chebyshev's answer, and explain why the usual answer based on the normal approximation (about 2 400) is much smaller. Hint: Chebyshev uses only the variance bound p(1 − p) ≤ 1/4, while the normal approximation also assumes the shape of the distribution of the sum.

Probabilistic Method

Why. Sometimes the easiest way to prove that an object exists is to show that a random object has the desired property with positive probability. If a random choice succeeds with nonzero probability, then at least one successful choice must exist. Paul Erdős developed this into a powerful method. It proves the existence of graphs, codes, and colorings that nobody knows how to construct explicitly, and in computer science it also produces approximation guarantees for problems such as MAX-CUT and MAX-SAT.

What. The method has two basic forms:

1. Probability form. Pick an object at random. If P(the object is bad) < 1, then some good object exists. The union bound often provides the estimate: if there are bad events B₁, …, Bₘ and Σ P(Bᵢ) < 1, then with positive probability none of them happens.
2. Expectation form. If E[X] = μ, then some outcome has X(ω) ≥ μ, and some outcome has X(ω) ≤ μ. A random variable cannot be below its average everywhere.

Example 1: a lower bound for Ramsey numbers (5.4). R(k, k) is the smallest n such that every red/blue coloring of the edges of K_n contains a monochromatic K_k. Claim: R(k, k) > 2^{k/2} for k ≥ 3. In other words, K_n with n = ⌊2^{k/2}⌋ has a coloring without any monochromatic K_k.

  Color each edge of K_n red or blue independently, with probability 1/2 each.
  For a fixed set S of k vertices:   P(all C(k,2) edges of S have the same color) = 2 · 2^{−C(k,2)}

  Union bound over all C(n, k) sets S:
    P(some monochromatic K_k) ≤ C(n, k) · 2^{1 − C(k,2)}
                              < (nᵏ / k!) · 2^{1 − k(k−1)/2}
                              ≤ (2^{k²/2} / k!) · 2^{1 − k²/2 + k/2}      using n ≤ 2^{k/2}
                              = 2^{1 + k/2} / k!  <  1                   for k ≥ 3

  So some coloring has no monochromatic K_k, and R(k, k) > n.

This argument from 1947 shows that good colorings exist and that a random coloring is almost always one. Yet no one knows how to construct a coloring of comparable size explicitly.

Example 2: large cuts. A cut splits the vertices of a graph into two sides. Its size is the number of edges that go between the sides. Claim: every graph with m edges has a cut with at least m/2 edges. Put each vertex on a random side, independently with probability 1/2. Each edge crosses the cut with probability 1/2, so by linearity the expected cut size is m/2, and some cut is at least that large.

Example 3: MAX-3SAT. Every 3-CNF formula (1.1) with m clauses, each containing three distinct variables, has an assignment that satisfies at least 7m/8 clauses. Under a uniformly random assignment, each clause fails only when all three of its literals are false, which has probability 1/8. Consequently, every such formula with at most 7 clauses is satisfiable.

From existence to algorithms. The expectation form also yields algorithms. You can sample repeatedly until you find an object at least as good as the average (a Las Vegas algorithm, analyzed with the geometric distribution). Alternatively, the method of conditional expectations fixes the random choices one at a time, each time keeping the conditional expectation at or above the target. That makes the algorithm deterministic.

from math import comb
import random

# Ramsey: for each k, the largest n for which the union bound C(n, k)·2^(1 − C(k, 2)) stays below 1
for k in range(4, 11, 2):
    n = k
    while comb(n + 1, k) * 2.0 ** (1 - comb(k, 2)) < 1:
        n += 1
    print(k, n, int(2 ** (k / 2)))         # k, n with R(k, k) > n, and Erdős's simpler 2^(k/2)
# 4 6 4
# 6 17 8
# 8 42 16
# 10 100 32

rng = random.Random(7)
V = 12
edges = [(u, v) for u in range(V) for v in range(u + 1, V) if rng.random() < 0.4]

def cut_size(side):                        # side[v] is 0 or 1
    return sum(side[u] != side[v] for u, v in edges)

tries = 0
while True:                                # Las Vegas: a cut of size ≥ m/2 exists, so keep sampling until one appears
    tries += 1
    side = [rng.randrange(2) for _ in range(V)]
    if cut_size(side) >= len(edges) / 2:
        break
print(len(edges), cut_size(side), tries)   # 29 19 1

The output shows that the clean bound 2^{k/2} gives away a lot: the last step of the inequality chain discards the k! in the denominator, and the exact union bound certifies much larger n. A careful estimate of the same calculation gives R(k, k) > (1 + o(1)) · k · 2^{k/2} / (e√2). Even so, no lower bound with an exponential base larger than √2 is known, and no explicit construction comes close.

Key takeaways

• If a random object is good with positive probability, then a good object exists.
• Common tools: the union bound (Σ P(Bᵢ) < 1, so no bad event is possible) and the expectation argument (some outcome is at least the mean).
• Results: R(k, k) > 2^{k/2}, every graph has a cut with at least m/2 edges, and every 3-CNF formula has an assignment satisfying at least 7/8 of its clauses.
• Existence proofs of this kind often become algorithms by resampling or by the method of conditional expectations.

Practice

1. (Beginner) Show that every graph with m edges has a bipartite subgraph with at least m/2 edges. Hint: this is the cut argument stated differently.
2. (Intermediate) Prove the MAX-3SAT claim: a 3-CNF formula with m clauses, each with three distinct variables, has an assignment satisfying at least 7m/8 clauses. Give a formula with 8 clauses that is unsatisfiable.
3. (Interview) A tournament is a complete graph in which every edge is directed, so that every pair of players has a winner. Show that if C(n, k)(1 − 2⁻ᵏ)^{n−k} < 1, then some tournament on n players has the following property: for every set of k players, some other player beats all of them. Hint: direct each edge at random. For a fixed set of k players, each outside player beats all of them with probability 2⁻ᵏ, independently of the other outside players. Then apply the union bound over all sets.

Randomized Algorithms

Why. An algorithm that makes random choices cannot be forced into its worst case by a particular input, because the adversary who chooses the input does not know which coins will be flipped. Randomness also gives simple, fast solutions where deterministic ones are complex or slow: primality testing (Miller–Rabin, 4.4), hashing (4.4), verifying matrix products, min-cut, load balancing, and sampling from streams. The tools of this chapter (indicators, linearity, Markov's inequality, and independent repetition) are exactly what is needed to analyze them.

What. A randomized algorithm reads random bits in addition to its input. Its guarantees hold for every input, with the probability taken over its own random choices, not over a distribution of inputs. There are two kinds:

	Las Vegas	Monte Carlo
Correctness	Always correct	Correct with high probability
Running time	Random. Analyze its expectation.	Bounded in the worst case
Examples	Randomized quicksort, rejection sampling, random-pivot quickselect	Miller–Rabin, Freivalds' check, Karger's min-cut, Bloom filters
How to improve it	Use Markov's inequality, with restarts, to meet a deadline	Repeat independently: an error of ε per run becomes εᵏ after k runs, for one-sided error
Converting to the other kind	Stop at a time limit (gives Monte Carlo)	If answers can be checked, repeat until a check passes (gives Las Vegas)

Example 1: randomized quicksort (Las Vegas). Choose the pivot uniformly at random. Let z₁ < z₂ < ⋯ < zₙ be the elements in sorted order. The elements zᵢ and zⱼ (i < j) are compared at most once, and they are compared exactly when the first pivot chosen from {zᵢ, …, zⱼ} is zᵢ or zⱼ. If any element strictly between them is chosen first, they end up on different sides and are never compared. All j − i + 1 elements of that range are equally likely to be chosen first, so with one indicator per pair:

                                   2                n  1
  E[comparisons] = Σ_{i<j}  ─────────── = 2(n + 1)·Σ  ─  − 4n = 2(n + 1)Hₙ − 4n ≈ 2n ln n ≈ 1.39 n log₂ n
                             j − i + 1             k=1 k

This holds for every input, including the already sorted input that drives a fixed "first element" pivot to Θ(n²) comparisons.

import java.util.Random;

public class RandomizedQuicksort {
    static long comparisons = 0;

    static void sort(int[] a, int lo, int hi, Random rng) {        // sorts a[lo..hi]
        if (lo >= hi) return;
        int r = lo + rng.nextInt(hi - lo + 1);                      // uniformly random pivot position
        swap(a, r, hi);
        int pivot = a[hi], i = lo;
        for (int j = lo; j < hi; j++) {                             // Lomuto partition
            comparisons++;                                          // each element is compared with the pivot once
            if (a[j] < pivot) swap(a, i++, j);
        }
        swap(a, i, hi);                                             // the pivot lands in its final position
        sort(a, lo, i - 1, rng);
        sort(a, i + 1, hi, rng);
    }

    static void swap(int[] a, int i, int j) { int t = a[i]; a[i] = a[j]; a[j] = t; }

    public static void main(String[] args) {
        int n = 1000, trials = 2000;
        Random rng = new Random(42);
        double harmonic = 0;
        for (int k = 1; k <= n; k++) harmonic += 1.0 / k;
        long total = 0;
        for (int t = 0; t < trials; t++) {
            int[] a = new int[n];
            for (int k = 0; k < n; k++) a[k] = k;                   // sorted input: the worst case for a fixed pivot
            comparisons = 0;
            sort(a, 0, n - 1, rng);
            total += comparisons;
        }
        System.out.printf("average %.0f, formula %.0f%n", (double) total / trials, 2 * (n + 1) * harmonic - 4 * n);
        // average 11004, formula 10986   (the standard error of the average is about 15)
    }
}

Example 2: Freivalds' algorithm (Monte Carlo). Checking whether A·B = C for n × n matrices by multiplying them takes Θ(n³) time, or about n^{2.37} with the best known algorithms. Freivalds' algorithm instead picks a random vector r ∈ {0, 1}ⁿ and checks whether A(Br) = Cr, which takes three matrix–vector products and O(n²) time.

• If AB = C, the check always passes.
• If AB ≠ C, let D = AB − C, and suppose D has a nonzero entry d_{ij}. Fix every coordinate of r except rⱼ. At most one of the two values of rⱼ can make row i of Dr equal to 0. So the check fails with probability at least 1/2.

With k independent rounds, a wrong C passes with probability at most 2⁻ᵏ. The error is one-sided: a "no" answer is always correct.

import random

def matvec(M, v):
    return [sum(m * x for m, x in zip(row, v)) for row in M]

def freivalds(A, B, C, rounds, rng):
    """Monte Carlo check of A·B == C. A False answer is always right; True is wrong with probability ≤ 2^-rounds."""
    n = len(A)
    for _ in range(rounds):
        r = [rng.randrange(2) for _ in range(n)]            # random 0/1 vector
        if matvec(A, matvec(B, r)) != matvec(C, r):         # three O(n²) matrix–vector products
            return False                                     # found a witness: A·B ≠ C for certain
    return True

rng = random.Random(9)
n = 60
A = [[rng.randrange(10) for _ in range(n)] for _ in range(n)]
B = [[rng.randrange(10) for _ in range(n)] for _ in range(n)]
C = [[sum(A[i][k] * B[k][j] for k in range(n)) for j in range(n)] for i in range(n)]   # the true product
bad = [row[:] for row in C]
bad[17][42] += 1                                             # one wrong entry

print(freivalds(A, B, C, 20, rng))                           # True
caught = sum(not freivalds(A, B, bad, 1, rng) for _ in range(10_000))
print(round(caught / 10_000, 2))                             # 0.51   a single round catches it half the time
print(sum(freivalds(A, B, bad, 20, rng) for _ in range(1000)))   # 0    20 rounds: fooled with probability ≤ 2^-20

Other randomized algorithms covered in these notes.

Algorithm	Kind	Guarantee	Where
Miller–Rabin primality test	Monte Carlo, one-sided	A composite number passes k rounds with probability at most 4⁻ᵏ	4.4
Universal hashing	Randomized data structure	P(collision) ≤ 1/m per pair, so expected chain length ≤ 1 + n/m for every key set	4.4, Indicator Random Variables
Fisher–Yates shuffle	Exact sampler	Each of the n! permutations has probability 1/n!	2.2, 5.5
Random derangement by rejection	Las Vegas	About e attempts on average	5.5, Geometric Distribution
Random cut	Las Vegas	Finds a cut with at least m/2 edges after a few samples on average	Probabilistic Method

Key takeaways

• Randomized guarantees hold for every input. The probability is over the algorithm's own coin flips, so no particular input is bad.
• Las Vegas algorithms are always correct and have a random running time. Monte Carlo algorithms have a bounded running time and may be wrong with small probability.
• Repeating a Monte Carlo algorithm with one-sided error ε gives error εᵏ after k runs, so error probabilities shrink exponentially.
• Randomized quicksort makes 2(n + 1)Hₙ − 4n ≈ 1.39 n log₂ n comparisons on average, by indicators and linearity. Freivalds' algorithm checks a matrix product in O(n²) time per round.

Practice

1. (Beginner) A Monte Carlo algorithm with one-sided error is wrong with probability at most 1/3 per run. How many independent runs reduce the error probability below 10⁻⁶?
2. (Intermediate) Reservoir sampling. A stream of unknown length is read once. Keep the first item. When item i arrives (i ≥ 2), replace the kept item with item i with probability 1/i. Prove by induction that after n items, each item is the kept one with probability exactly 1/n.
3. (Interview) Why does randomized quicksort's O(n log n) expected bound hold for every input, while quicksort with the first element as pivot has inputs that force Θ(n²) comparisons? Can an adversary still make randomized quicksort slow, and what changes if the adversary can predict your random number generator? Hint: separate "expected over inputs" from "expected over coin flips". Use Markov's inequality for the probability of a slow run. Consider what a seeded, predictable PRNG gives an attacker (hash-flooding attacks are the analogous problem for hash tables).

<a id="8-graph-theory"></a>