14. Probability Distributions

Level: Intermediate You'll use this for: modeling request arrivals, error rates, time between failures, and measurement noise.

The key distributions, and when each shows up

Binomial

P(k successes in n trials) = C(n,k) · p^k · (1-p)^(n-k)

from scipy.stats import binom
# 1000 requests, each independently has 0.1% error rate — P(exactly 3 errors)?
binom.pmf(3, n=1000, p=0.001)

A small binomial computed fully by hand — 5 coin flips (p=0.5), what's P(exactly 2 heads)?

C(5,2) = 5! / (2!·3!) = 120 / (2×6) = 10     <- how many ways to choose which 2 of 5 flips are heads (Ch. 15)
p^k = 0.5² = 0.25                             <- probability those 2 specific flips are heads
(1-p)^(n-k) = 0.5³ = 0.125                    <- probability the other 3 specific flips are tails

P(exactly 2 heads) = 10 × 0.25 × 0.125 = 0.3125

Every binomial probability is this same three-part product: "how many arrangements achieve this outcome" × "probability of the success pattern" × "probability of the failure pattern."

Poisson — the request-arrival workhorse

Models "how many events happen in a fixed window" when events happen independently at a constant average rate λ.

P(k events) = (λ^k · e^-λ) / k!

from scipy.stats import poisson
# API averages 100 requests/sec. What's P(more than 130 arrive in one second)?
1 - poisson.cdf(130, mu=100)

A small Poisson computed by hand — a server averages λ=2 crashes per week, what's P(exactly 3 crashes this week)?

λ^k = 2³ = 8
e^-λ = e^-2 ≈ 0.1353         <- a fixed constant once λ is fixed, from a table or calculator
k! = 3! = 3×2×1 = 6

P(3 crashes) = (8 × 0.1353) / 6 = 1.0824 / 6 ≈ 0.1804

Same shape every time: raise λ to the k, multiply by the constant e^-λ, divide by k! to account for the fact that the k events could have landed in any order within the window.

Why this matters directly for capacity planning: if you provision for exactly the average load, you will be under capacity roughly half the time, because real traffic follows a Poisson-like distribution with variance, not a flat line. This is the mathematical justification for "provision for p99 load, not average load" (Ch. 25 builds directly on this with queueing theory).

Normal distribution and the Central Limit Theorem

The Central Limit Theorem says: the sum (or average) of many independent random variables tends toward a normal distribution, regardless of the shape of the original distribution. This is why so many real-world measurements (and so many statistical tests, Ch. 22) assume normality — it's not an arbitrary choice, it's a consequence of aggregating many small independent effects.

Exponential distribution — time between failures

If failures follow a Poisson process (constant hazard rate), the time between failures follows an exponential distribution. This is the basis of MTBF (Mean Time Between Failures) reasoning in reliability engineering.

from scipy.stats import expon
# Mean time between failures = 720 hours. P(next failure within 100 hours)?
expon.cdf(100, scale=720)

PMF, PDF, CDF: three ways to describe a distribution

PMF (discrete):    p(x) = P(X = x)                       Σₓ p(x) = 1
PDF (continuous):  f(x),   P(a ≤ X ≤ b) = ∫ₐᵇ f(x) dx     ∫ f(x) dx = 1
CDF (both):        F(x) = P(X ≤ x)                       F'(x) = f(x)

The CDF is the one engineers already use: a percentile is the inverse CDF. "p99 latency = 250 ms" means F(250 ms) = 0.99.

The normal density, every symbol explained

f(x) = 1/(σ·√(2π)) · e^( -(x - μ)² / (2σ²) )

Example: API latency ~ N(μ = 100 ms, σ = 20 ms). What fraction of requests exceed 140 ms?

z = (140 - 100) / 20 = 2       ->   P(Z > 2) = 1 - Φ(2) ≈ 0.0228   (2.3%)

Φ (capital phi) is the standard normal CDF. You look it up in a table or call scipy.stats.norm.cdf(2). The e^(-x²) shape decides how quickly tails die. At 4σ the probability is 0.003%, and at 6σ it's one in a billion. Real latency is not normal (see heavy tails below), and that difference is where outages come from.

Every distribution's mean and variance at a glance

Notice the Poisson's mean equals its variance. That gives you a free diagnostic. If your "requests per second" metric has a variance much larger than its mean (overdispersion), arrivals aren't independent. That points to bursty clients, retry storms, or cron jobs firing on the minute.

The Poisson limit theorem: where Poisson comes from

Theorem. If n → ∞ and p → 0 with n·p = λ held fixed, then Binomial(n, p) → Poisson(λ).

Many independent chances, each tiny, give Poisson. That's why it models requests (millions of users, each rarely clicking), disk failures (many disks, each rarely failing) and radioactive decay alike.

Example: 1000 requests at a 0.1% error rate, so λ = 1000·0.001 = 1.

Binomial: P(exactly 3) = C(1000,3)·0.001³·0.999⁹⁹⁷ = 0.06128
Poisson:  P(exactly 3) = 1³·e⁻¹/3!                 = 0.06131

The two agree to three decimal places, and the Poisson version can be worked out by hand.

Memorylessness: the exponential distribution's strange property

P(T > s + t | T > s) = P(T > t)

Proof in one line, using P(T > t) = e^(-λt) from the CDF:

P(T > s+t | T > s) = P(T > s+t) / P(T > s) = e^(-λ(s+t)) / e^(-λs) = e^(-λt)

If failures are exponential, a server that has run for 1,000 hours is exactly as likely to fail in the next hour as a brand-new one. That's why time-based replacement of exponentially-failing parts is pointless. It's also why "it's been quiet for an hour, so a request must be due" is the gambler's fallacy again. Real hardware isn't fully memoryless. Disks follow a "bathtub curve" with high infant mortality, a flat middle and a wear-out phase, and Backblaze's data shows exactly that shape.

The Law of Large Numbers and the Central Limit Theorem, stated properly

Let X₁, …, Xₙ be independent draws from any distribution with mean μ and standard deviation σ, and let X̄ₙ = (X₁ + … + Xₙ)/n be their average.

Law of Large Numbers. X̄ₙ → μ as n → ∞. Averages settle on the true mean.

Central Limit Theorem. For large n,

(X̄ₙ - μ) / (σ / √n)   ≈   N(0, 1)

Averages are approximately normal, centred on μ, with spread σ/√n, whatever the shape of the original distribution.

Example: individual request latencies are exponential with mean 50 ms, so σ = 50 ms as well. That's strongly skewed and nothing like a bell. Now average them in batches of 100:

standard error = 50 / √100 = 5 ms
95% of batch averages fall in  50 ± 1.96·5  =  [40.2, 59.8] ms

A simulation of 100,000 such batches gives a mean of 49.99, a standard deviation of 5.01, and 95.0% inside that interval. The skew has disappeared, as the theorem promises.

import numpy as np
rng = np.random.default_rng(0)
batch_means = rng.exponential(scale=50, size=(100_000, 100)).mean(axis=1)
print(batch_means.mean(), batch_means.std())          # ~50.0, ~5.0
print(np.mean(np.abs(batch_means - 50) < 1.96 * 5))   # ~0.95

The CLT is why A/B tests (Ch. 22) can use normal-based z-tests on conversion rates, and why mini-batch gradients in ML are noisy estimates of the true gradient whose noise shrinks like 1/√batch_size. Doubling the batch cuts gradient noise by only about 30%, which is why "just use bigger batches" has diminishing returns.

Concentration inequalities: guarantees without assuming normality

The CLT says "approximately, for large n". Sometimes you need a guarantee that holds for any distribution. These three inequalities give one:

Markov:      P(X ≥ a)          ≤ E[X] / a                 (X ≥ 0)
Chebyshev:   P(|X - μ| ≥ kσ)   ≤ 1 / k²
Hoeffding:   P(|X̄ₙ - μ| ≥ ε)   ≤ 2·e^(-2nε²)               (each Xᵢ in [0, 1])

Markov example: average queue length is 4. Without knowing anything else, P(queue ≥ 40) ≤ 4/40 = 10%. It's a crude bound, but it holds for any non-negative quantity.

Chebyshev example: for any distribution, at most 1/9 ≈ 11% of values lie beyond 3σ. The normal distribution puts only 0.3% out there, so a "3-sigma alert" on non-normal data can fire about 40× more often than you designed for.

Hoeffding, in AI: how big should an eval set be? You want the model's measured accuracy within ε = 2 points of the truth, with 95% confidence (failure probability δ = 0.05). Set 2e^(-2nε²) = δ and solve:

n ≥ ln(2/δ) / (2ε²) = ln(40) / (2 · 0.02²) = 3.689 / 0.0008 ≈ 4,612 examples

Run it the other way and a 500-question benchmark only pins accuracy to ε = √(ln 40 / 1000) ≈ ±6 points. A leaderboard gap of 3 points on such a benchmark is noise. Hoeffding is also the backbone of PAC learning theory, the formal answer to "how much data do I need to learn?"

Heavy tails and the "tail at scale"

Many systems quantities follow a Pareto (power-law) distribution, where the tail decays polynomially instead of exponentially:

P(X > x) = (x_m / x)^α        for x ≥ x_m

With α ≈ 1.16 you get the "80/20 rule": 20% of the flows carry 80% of the bytes. Internet flow sizes, file sizes and request costs all look like this. In networking these are the "elephant" flows (few, huge) and "mice" flows (many, tiny), and datacenter switches schedule them differently.

Tail latency compounds under fan-out. A search query fans out to 100 leaf servers and must wait for the slowest. If each leaf is slow (above its p99) independently 1% of the time:

P(query is slow) = 1 - P(all 100 fast) = 1 - 0.99¹⁰⁰ ≈ 63%

One leaf's p99 becomes the whole service's median. This is the central result of Google's "The Tail at Scale" (Dean & Barroso, 2013). It's why hedged requests (send a backup after the p95 has passed) and tied requests exist, and why p99.9 matters more than p50 in microservice architectures.

AI in the wild: sampling text is sampling a distribution

At every step an LLM outputs logits z, turns them into a categorical distribution with softmax (Ch. 13), and draws a random sample from it. Temperature T reshapes that distribution before sampling:

pᵢ = e^(zᵢ / T) / Σⱼ e^(zⱼ / T)

Logits z = (2.0, 1.0, 0.1) for three candidate tokens:

As T → 0 this becomes argmax (greedy decoding). As T → ∞ it becomes uniform. Top-p (nucleus) sampling truncates to the smallest set of tokens whose probabilities sum to p (say 0.9), which is just the CDF from earlier in this chapter applied to the sorted probabilities.

Real-life engineering ties

• SWE/SRE: error-budget math, capacity planning against percentile — not average — load (Ch. 9, 25).
• AI: assumed noise models in regression (Gaussian noise), Poisson regression for count data (e.g., predicting clicks/purchases).
• Network: packet arrival modeling (Poisson), queueing theory inputs (Ch. 25 uses this directly — M/M/1 queues assume Poisson arrivals and exponential service times).
• Security: modeling time-between-attacks/scans for anomaly baselining; Poisson-based alerting ("we normally see ~5 failed logins/hour — 200 in the last hour is statistically extreme," tying back to Ch. 9's 3-sigma idea).

More real-life examples:

• Telephone exchanges — and modern call-center staffing tools — use the Erlang C formula, built on the same Poisson-arrival assumption as this chapter, to compute how many agents are needed to hit a target wait-time SLA; the unit "Erlang" is named for the mathematician who solved this exact problem in 1917.
• Backblaze's publicly published hard-drive failure statistics report Annualized Failure Rate, used with an exponential (constant hazard rate) model to estimate the probability a specific drive fails within the next year — the same math as the MTBF example above, with real, public data.
• Insurance actuaries model claim frequency with Poisson distributions and claim severity with heavy-tailed distributions (log-normal or Pareto) — nearly identical in structure to the risk model in Ch. 26.
• A/B testing platforms model each user's conversion as a Bernoulli trial and the aggregate conversion count as Binomial — the statistical foundation underneath every "statistically significant" banner those tools show.

Common mistake

Applying Normal-distribution intuition (symmetric, thin tails) to data that's actually heavy-tailed (network latency, financial returns, security incident severity) — heavy-tailed data has "black swan" events far more often than a Normal model predicts, badly underestimating tail risk.