Probability & statistics
AI systems are probabilistic: a classifier outputs P(class | input), an LLM samples the next token from a distribution, and every eval number is an estimate with error bars. This lesson builds the theory — probability axioms, Bayes, random variables and distributions, expectation/variance, the law of large numbers & CLT, maximum likelihood, and hypothesis testing — and shows why temperature sampling works and when an A/B result is real.
random, math) so every result is reproducible — each block seeds the RNG, so the printed numbers below are exactly what you will get. The A/B-test section uses numpy for convenience and is labelled numpy-required; the same math is done with the stdlib elsewhere. No fabricated numbers.Learning objectives
- State the probability axioms and compute with conditional probability.
- Apply Bayes’ theorem and connect it to a probabilistic classifier.
- Work with Bernoulli, Binomial and Normal random variables.
- Compute expectation and variance and interpret them.
- Explain why temperature/top-p sampling works in LLM decoding.
- Demonstrate the law of large numbers and the CLT by simulation.
- Derive maximum-likelihood estimation and link it to model training.
- Run a hypothesis test + confidence interval for an A/B experiment and judge significance.
1 · Axioms & conditional probability intermediate
Probability assigns each event a number in [0, 1]. The axioms (Kolmogorov): P(Ω) = 1 for the whole sample space, P(A) ≥ 0, and for disjoint events P(A ∪ B) = P(A) + P(B). Conditional probability P(A | B) = P(A ∩ B) / P(B) is the probability of A once you know B happened — the update that all inference is built on. Events are independent when P(A ∩ B) = P(A)P(B).
python# stdlib only
# A fair pair of dice; sample space = 36 equally likely outcomes.
outcomes = [(i, j) for i in range(1, 7) for j in range(1, 7)]
def prob(event): # event: a predicate on (i, j)
hits = [o for o in outcomes if event(o)]
return len(hits) / len(outcomes)
pA = prob(lambda o: o[0] + o[1] == 8) # sum is 8
pB = prob(lambda o: o[0] == 5) # first die is 5
pAB = prob(lambda o: o[0] + o[1] == 8 and o[0] == 5)
print("P(sum=8) ", round(pA, 4)) # 0.1389
print("P(first=5) ", round(pB, 4)) # 0.1667
print("P(sum=8|first=5)", round(pAB / pB, 4)) # 0.1667
P(sum=8) 0.1389
P(first=5) 0.1667
P(sum=8|first=5) 0.1667
P(sum=8) = 5/36 ≈ 0.1389. Given the first die is 5, the second must be 3 — exactly one of six equally likely values — so P(sum=8 | first=5) = 1/6 ≈ 0.1667. Conditioning shrank the sample space from 36 outcomes to the 6 where the first die is 5. That shrink-the-space move is the engine of the Bayes update next.
2 · Bayes’ theorem — the classifier’s update advanced
Bayes’ theorem inverts a conditional: P(H | E) = P(E | H)·P(H) / P(E), where P(E) = Σh P(E | h)P(h) (the law of total probability). It turns a prior P(H) plus likelihood P(E | H) into a posterior P(H | E). Every probabilistic classifier — naive Bayes, and conceptually any model that outputs P(class | input) — is doing this update. The classic trap: a very accurate test for a rare condition still yields mostly false positives.
python# stdlib only
p_D = 0.01 # prior: 1% of people have the condition
sens = 0.99 # P(test + | disease) (true-positive rate)
false_pos = 0.05 # P(test + | no disease) (1 - specificity)
# Law of total probability for the denominator P(test +):
p_pos = sens * p_D + false_pos * (1 - p_D)
# Bayes:
p_D_given_pos = sens * p_D / p_pos
print("P(test +) ", round(p_pos, 5)) # 0.0594
print("P(disease | +) ", round(p_D_given_pos, 5)) # 0.16667
P(test +) 0.0594
P(disease | +) 0.16667
Despite a 99%-sensitive test, a positive result means only a 16.7% chance of disease — because the condition is rare, the 5% false positives among the 99% healthy people swamp the true positives. This base-rate reasoning is why a classifier tuned only for high recall on a rare class floods you with false alarms, and why you report precision, not just accuracy (see Ch 5).
A spam filter flags 95% of spam (P(flag | spam)=0.95) and wrongly flags 2% of ham (P(flag | ham)=0.02). If 20% of mail is spam, what fraction of flagged mail is actually spam?
Show answer
3 · Random variables & distributions advanced
A random variable maps outcomes to numbers. Three you must know: Bernoulli(p) — a single 0/1 trial (a coin, a click, a token being ‘correct’); Binomial(n, p) — the count of successes in n independent Bernoulli trials, PMF C(n,k)pk(1−p)n−k; and the Normal(μ, σ²) — the bell curve that the CLT (§5) makes universal. Discrete variables have a PMF (probability per value); continuous ones have a PDF (density).
python# stdlib only
import math
from math import comb
def binom_pmf(k, n, p):
return comb(n, k) * p**k * (1 - p)**(n - k)
n, p = 10, 0.3
print("P(3 successes in 10):", round(binom_pmf(3, n, p), 5)) # 0.26683
print("mean = n*p :", n * p) # 3.0
print("var = n*p*(1-p) :", round(n * p * (1 - p), 4)) # 2.1
def normal_pdf(x, mu=0.0, sigma=1.0):
return math.exp(-0.5 * ((x - mu) / sigma)**2) / (sigma * math.sqrt(2*math.pi))
print("N(0,1) density at 0 :", round(normal_pdf(0.0), 5)) # 0.39894
P(3 successes in 10): 0.26683
mean = n*p : 3.0
var = n*p*(1-p) : 2.1
N(0,1) density at 0 : 0.39894
The most likely single count for Binomial(10, 0.3) is around its mean np = 3, with PMF ≈ 0.267. The standard Normal peaks at 1/√(2π) ≈ 0.399. In AI, per-token correctness is Bernoulli, the number of correct answers on an eval set of n items is Binomial, and by the CLT your accuracy estimate is approximately Normal — which is what lets you put error bars on it (§7).
| Distribution | Models | Mean | Variance | AI use |
|---|---|---|---|---|
| Bernoulli(p) | one 0/1 trial | p | p(1−p) | single correct/incorrect, one click |
| Binomial(n,p) | successes in n trials | np | np(1−p) | #correct on an eval set of n |
| Normal(μ,σ²) | sums/means (via CLT) | μ | σ² | metric estimates, error bars |
4 · Expectation, variance & sampling — why temperature works advanced
Expectation E[X] = Σ x·P(x) is the long-run average; variance Var(X) = E[(X−μ)²] is the spread. LLM decoding is sampling from a distribution: the model produces logits, softmax(logits / T) turns them into a probability vector, and the next token is drawn from it. Temperature T reshapes that distribution — T<1 sharpens it (greedier, less variance), T>1 flattens it (more random). Top-p truncates to the smallest set of tokens whose probability sums to p, then renormalizes. This is applied conditional probability.
python# stdlib only
import math, random
def softmax(logits, T=1.0):
z = [x / T for x in logits]
m = max(z) # stability (see A4)
e = [math.exp(x - m) for x in z]
s = sum(e)
return [x / s for x in e]
logits = [2.0, 1.0, 0.1] # raw next-token scores
for T in (0.5, 1.0, 2.0):
probs = softmax(logits, T)
print(f"T={T}: {[round(x, 3) for x in probs]}")
# Draw one token from a distribution (inverse-CDF sampling):
def sample(probs, r):
c = 0.0
for i, pr in enumerate(probs):
c += pr
if r <= c:
return i
return len(probs) - 1
random.seed(7)
probs = softmax(logits, 1.0)
draws = [sample(probs, random.random()) for _ in range(10000)]
freq0 = draws.count(0) / len(draws)
print("empirical P(token 0):", round(freq0, 4), "vs", round(probs[0], 4))
T=0.5: [0.864, 0.117, 0.019]
T=1.0: [0.659, 0.242, 0.099]
T=2.0: [0.502, 0.304, 0.194]
empirical P(token 0): 0.6584 vs 0.659
At T=0.5 the top token gets 86% of the mass (near-greedy); at T=2.0 the distribution flattens toward uniform (more ‘creative’/random). The sampling check confirms the theory: drawing 10,000 tokens gives token 0 about 66% of the time, matching its 0.659 probability — the law of large numbers (§5) in action. This is the temperature knob the API exposes; see A4 for the decoding loop.
5 · Law of large numbers & the CLT expert advanced
The law of large numbers (LLN): the sample mean converges to the true mean as n grows — why more eval examples give a more trustworthy score. The central limit theorem (CLT): the sample mean of any finite-variance distribution is approximately Normal for large n, with variance σ²/n. Together they are why we can estimate metrics by averaging and attach Normal error bars — the foundation of the confidence intervals in §7.
python# stdlib only
import random
# --- LLN: sample mean of a fair coin -> 0.5 as n grows ---
random.seed(42)
def coin_mean(N):
return sum(1 for _ in range(N) if random.random() < 0.5) / N
print("LLN n=100 :", round(coin_mean(100), 3)) # 0.5
print("LLN n=10_000 :", round(coin_mean(10_000), 4)) # 0.4986
print("LLN n=1e6 :", round(coin_mean(1_000_000), 5)) # 0.49982
# --- CLT: mean of 30 Uniform(0,1) draws is ~Normal ---
random.seed(0)
def sample_mean(n=30):
return sum(random.random() for _ in range(n)) / n
means = [sample_mean() for _ in range(5000)]
mu = sum(means) / len(means)
var = sum((m - mu)**2 for m in means) / len(means)
print("CLT mean ~0.5:", round(mu, 4)) # 0.4998
print("CLT var ~sigma^2/n:", round(var, 6),
"vs theory", round((1/12) / 30, 6)) # 0.002702 vs 0.002778
LLN n=100 : 0.5
LLN n=10_000 : 0.4986
LLN n=1e6 : 0.49982
CLT mean ~0.5: 0.4998
CLT var ~sigma^2/n: 0.002702 vs theory 0.002778
As n grows the coin’s sample mean tightens onto 0.5 (LLN). And the mean of 30 uniforms — a flat distribution — produces a bell-shaped spread centered at 0.5 with variance ≈ σ²/n (Uniform variance 1/12, divided by 30 = 0.002778; simulated 0.002702). That the average of a non-Normal variable becomes Normal is the CLT, and it is precisely why your eval accuracy — an average of 0/1 correctness — has a Normal-shaped uncertainty you can bound.
6 · Maximum likelihood estimation expert expert
Maximum likelihood estimation (MLE) picks the parameter that makes the observed data most probable: maximize the likelihood L(θ) = ∏ P(xi | θ), or equivalently the log-likelihood (turns the product into a sum, avoids underflow). For a Bernoulli, the MLE of p is just the sample mean — provable by setting the derivative of the log-likelihood to zero. This is what training is: minimizing negative log-likelihood (cross-entropy, see A12) is maximizing the likelihood of the training data under the model.
python# stdlib only
import math
data = [1, 1, 0, 1, 0, 0, 1, 1, 1, 0] # 6 ones out of 10
def log_likelihood(p, data):
if p <= 0 or p >= 1: # avoid log(0)
return float("-inf")
return sum(math.log(p) if x else math.log(1 - p) for x in data)
# Analytic MLE for Bernoulli is the sample mean:
p_hat = sum(data) / len(data)
print("analytic MLE p:", p_hat) # 0.6
# Numeric grid search should agree:
best = max((round(0.01*i, 2) for i in range(1, 100)),
key=lambda p: log_likelihood(p, data))
print("numeric MLE p:", best) # 0.6
print("log-likelihood:", round(log_likelihood(p_hat, data), 5)) # -6.73012
analytic MLE p: 0.6
numeric MLE p: 0.6
log-likelihood: -6.73012
Six ones out of ten → the MLE is p̂ = 0.6, and a numeric search over the log-likelihood lands on the same value — confirming the analytic result. Training a neural network is the same principle at scale: adjust millions of parameters to maximize the log-likelihood of the training tokens, i.e. minimize cross-entropy loss. MLE is the bridge from probability to optimization (A12).
7 · Hypothesis testing & confidence intervals — is the A/B win real? expert expert
You ship a new prompt and conversion rises from 10% to 13%. Real, or noise? A hypothesis test asks: under the null hypothesis (no difference), how surprising is a gap this big? The p-value is that probability; below a threshold (commonly 0.05) you call it significant. A confidence interval gives the plausible range for the true difference. This is exactly how you decide whether an eval or A/B improvement (Ch 5) is a genuine gain or sampling luck.
python# numpy-required
import numpy as np
from math import erf, sqrt
def norm_cdf(x): # standard-normal CDF via erf
return 0.5 * (1 + erf(x / sqrt(2)))
# Control: 100/1000 convert. Variant: 130/1000 convert.
n_c, x_c = 1000, 100
n_v, x_v = 1000, 130
p_c, p_v = x_c / n_c, x_v / n_v
# Pooled two-proportion z-test (H0: p_c == p_v):
p_pool = (x_c + x_v) / (n_c + n_v)
se = sqrt(p_pool * (1 - p_pool) * (1/n_c + 1/n_v))
z = (p_v - p_c) / se
p_value = 2 * (1 - norm_cdf(abs(z))) # two-sided
print("rates:", p_c, p_v)
print("z =", round(z, 4), " p =", round(p_value, 5)) # z=2.1027 p=0.03549
# 95% CI for the difference (unpooled SE):
se_d = sqrt(p_c*(1-p_c)/n_c + p_v*(1-p_v)/n_v)
lo, hi = (p_v - p_c) - 1.96*se_d, (p_v - p_c) + 1.96*se_d
print("95% CI on lift:", (round(lo, 5), round(hi, 5))) # (0.00207, 0.05793)
rates: 0.1 0.13
z = 2.1027 p = 0.03549
95% CI on lift: (0.00207, 0.05793)
The z-statistic is 2.10, giving a two-sided p-value of 0.035 — below 0.05, so the 3-point lift is statistically significant: unlikely to be noise. The 95% CI on the lift is roughly (0.2%, 5.8%) — it excludes 0, the same conclusion, and additionally tells you the true gain could be as small as a fraction of a point. Report the CI, not just ‘significant’: a significant but tiny effect may not be worth shipping.
Checkpoint expert
✓ Checkpoint — you can move on when you can…
- State the probability axioms and compute a conditional probability.
- Apply Bayes’ theorem and explain the base-rate trap for rare classes.
- Give the mean and variance of Bernoulli, Binomial and Normal variables.
- Explain how temperature reshapes the LLM sampling distribution.
- Demonstrate LLN and CLT and state Var(mean) = σ²/n.
- Derive the Bernoulli MLE and connect NLL minimization to training.
- Run a two-proportion test + CI and judge whether an A/B lift is real.
Your eval set has 200 items and your agent scores 75% accuracy. Roughly how wide is the 95% confidence interval, and what does that imply about comparing it to a 78% baseline?
Show answer
🪜 Practice ladder beginner → industry
Six graded exercises, easy to real-world. Try each before opening its solution.
Context: Conditioning — restricting attention to the cases where B happened — is the core move behind every model that outputs P(y | x). Simulating it makes the definition concrete.
Your task: By simulating two dice, estimate P(sum = 8 | first die = 5) and compare to the exact value 1/6.
Requirements:
- Use
randomwith a fixed seed - Keep only trials where the first die is 5, then measure the fraction with sum 8
- Print the estimate and the exact 1/6
- No numpy
💡 Hint: Filter to the conditioning event first, then compute the fraction of those with the target event.
Show solution
Restrict to the conditioning event, then measure the target within it.
import random
random.seed(1)
N = 200_000
d = [(random.randint(1,6), random.randint(1,6)) for _ in range(N)]
given = [o for o in d if o[0] == 5] # first die = 5
hits = [o for o in given if sum(o) == 8]
print(round(len(hits)/len(given), 4), round(1/6, 4))0.1689 0.1667The estimate lands within ~0.001 of 1/6 — and tightens further as N grows (LLN).
Context: Turning likelihoods and a prior into a posterior is the naive-Bayes classifier in one line; the base rate can make a ‘good’ detector untrustworthy.
Your task: Compute P(spam | flagged) given P(flag|spam)=0.95, P(flag|ham)=0.02, and prior P(spam)=0.20.
Requirements:
- Compute P(flag) via the law of total probability
- Apply Bayes to get the posterior
- Print the posterior rounded
- No numpy
💡 Hint: P(flag) = P(flag|spam)P(spam) + P(flag|ham)P(ham); then Bayes.
Show solution
Total probability for the denominator, then Bayes.
p_spam = 0.20
p_flag_spam, p_flag_ham = 0.95, 0.02
p_flag = p_flag_spam*p_spam + p_flag_ham*(1-p_spam)
post = p_flag_spam*p_spam / p_flag
print(round(p_flag, 4), round(post, 4)) # 0.206 0.92230.206 0.922392% — high, because spam is common here; the same test at 1% prevalence would be far weaker.
Context: LLM decoding samples tokens from softmax(logits/T); showing the empirical frequencies match the probabilities links sampling theory to real generation.
Your task: Given logits [2.0, 1.0, 0.1], build the T=1.0 distribution, draw 20,000 tokens, and confirm the empirical frequency of each token matches its probability.
Requirements:
- Implement a numerically stable softmax with temperature
- Sample via inverse-CDF from a uniform draw
- Seed the RNG for reproducibility
- Compare empirical frequencies to the probabilities
- No numpy
💡 Hint: Accumulate probabilities and pick the first index whose cumulative sum exceeds a uniform draw.
Show solution
Softmax gives the distribution; inverse-CDF sampling draws from it; LLN makes frequencies converge.
import math, random
def softmax(z, T=1.0):
z=[x/T for x in z]; m=max(z); e=[math.exp(x-m) for x in z]; s=sum(e)
return [x/s for x in e]
def sample(probs, r):
c=0.0
for i,pr in enumerate(probs):
c+=pr
if r <= c: return i
return len(probs)-1
random.seed(3)
probs = softmax([2.0,1.0,0.1], 1.0)
draws = [sample(probs, random.random()) for _ in range(20000)]
freq = [round(draws.count(i)/len(draws),3) for i in range(3)]
print('probs', [round(x,3) for x in probs]) # [0.659, 0.242, 0.099]
print('freq ', freq)probs [0.659, 0.242, 0.099]
freq [0.66, 0.243, 0.098]
Context: The CLT is why averaging a non-Normal quantity (like 0/1 correctness) yields a Normal estimate you can bound. Watching it emerge from a skewed source drives the point home.
Your task: Draw sample means of n=40 exponential-like variables (use max of two uniforms, which is skewed), 5000 times, and check the distribution of means is roughly Normal with variance σ²/n.
Requirements:
- Use a clearly non-Normal base distribution
- Average n=40 draws, repeated 5000 times
- Report the empirical mean and variance of the sample means
- Compare variance to σ²/n
- No numpy
💡 Hint: Var(mean) = Var(base)/n; estimate Var(base) from many single draws first.
Show solution
Even a skewed base yields Normal-shaped sample means with variance shrinking as σ²/n.
import random
random.seed(5)
def base(): # skewed: max of two uniforms
return max(random.random(), random.random())
singles = [base() for _ in range(200000)]
mu_b = sum(singles)/len(singles)
var_b = sum((x-mu_b)**2 for x in singles)/len(singles)
n = 40
means = [sum(base() for _ in range(n))/n for _ in range(5000)]
mu = sum(means)/len(means)
var = sum((m-mu)**2 for m in means)/len(means)
print(round(mu,4), round(var,6), round(var_b/n,6))0.6669 0.001393 0.001395The sample-mean variance (0.001393) matches σ²/n (0.001395), and the means cluster in a bell around 2/3 — the CLT, from a skewed source.
Context: Training minimizes negative log-likelihood; implementing MLE for a simple model shows what that objective computes and why the log turns a fragile product into a stable sum.
Your task: Given Bernoulli data, compute the analytic MLE (sample mean) and confirm it maximizes the log-likelihood via a grid search.
Requirements:
- Implement the log-likelihood as a sum of logs
- Compute the analytic MLE (sample mean)
- Grid-search p and confirm the maximizer matches
- Print the maximized log-likelihood
- No numpy
💡 Hint: log L(p) = (#ones) log p + (#zeros) log(1-p); its maximizer is #ones / n.
Show solution
The Bernoulli MLE is the sample mean; a grid search over the log-likelihood confirms it.
import math
data = [1,0,1,1,0,1,0,1,1,1,0,1] # 8 ones / 12
def ll(p):
if p <= 0 or p >= 1: return float('-inf')
return sum(math.log(p) if x else math.log(1-p) for x in data)
p_hat = sum(data)/len(data)
grid = [round(0.01*i,2) for i in range(1,100)]
best = max(grid, key=ll)
print(round(p_hat,4), best, round(ll(p_hat),5))0.6667 0.67 -7.63817Analytic 0.6667 and the grid’s 0.67 agree (grid step is 0.01) — minimizing NLL is maximizing this.
Context: This is the daily decision behind shipping a prompt/model change: is the observed lift real, and how big could it be? A two-proportion test plus CI is the standard tool.
Your task: Given control 100/1000 and variant 130/1000 conversions, run a two-proportion z-test, report the p-value, and give the 95% CI on the lift — then state the ship/hold decision.
Requirements:
- Compute the pooled-SE z-statistic and two-sided p-value
- Compute the 95% CI on the difference (unpooled SE)
- State whether p < 0.05 and whether the CI excludes 0
- Give a one-line decision that mentions effect size, not just significance
- numpy allowed and labelled
💡 Hint: Pooled SE for the test; unpooled SE for the CI; z=1.96 for 95%.
Show solution
Pooled-SE z-test for significance, unpooled-SE interval for the effect size.
# numpy-required (numpy optional here; math suffices)
from math import erf, sqrt
def ncdf(x): return 0.5*(1+erf(x/sqrt(2)))
n_c,x_c,n_v,x_v = 1000,100,1000,130
p_c,p_v = x_c/n_c, x_v/n_v
pp = (x_c+x_v)/(n_c+n_v)
se = sqrt(pp*(1-pp)*(1/n_c+1/n_v))
z = (p_v-p_c)/se; pval = 2*(1-ncdf(abs(z)))
se_d = sqrt(p_c*(1-p_c)/n_c + p_v*(1-p_v)/n_v)
lo,hi = (p_v-p_c)-1.96*se_d, (p_v-p_c)+1.96*se_d
print(round(z,4), round(pval,5)) # 2.1027 0.03549
print(round(lo,5), round(hi,5)) # 0.00207 0.05793
print('ship' if pval < 0.05 and lo > 0 else 'hold')2.1027 0.03549
0.00207 0.05793
shipSignificant (p=0.035) and the CI excludes 0, so ship — but the lift could be as small as ~0.2 points, so weigh it against rollout cost.