Probability & statistics

Most of the probability worth knowing here serves one practical question: my new version scored better, should I believe it? This section covers what you need to answer that, distributions and expectation, conditional probability and Bayes, and how much a measurement drifts purely because of which cases happened to be sampled. It is the difference between reporting a result and reporting a fluctuation.

What you actually need#

The probability that matters in practice is not the theory. It is being able to answer one question: my new version scored better, is that real?

A random variable, in plain words#

A random variable is neither a variable nor random in the everyday sense. It is a rule that turns something you cannot predict into a number. Roll a die and read the top face: that is a random variable. Serve a web request and record how many milliseconds it took: that is a random variable. Send a prompt to a model, have a grader score the answer out of two, and the score is a random variable. The rule never changes. What changes is which outcome the world hands you this time.

A distribution is the whole picture of one such variable: every value it can take, and how often. Write the values in one column and their long-run shares in another, and you have written the distribution down. The shares add up to one, because some outcome always happens. Nothing else is required of it, and the word carries no more mystery than that.

Make it concrete. A grader awards 0, 1 or 2 points per answer, and across everything your system produces, 20 per cent score 0, 50 per cent score 1 and 30 per cent score 2. Those three numbers are the distribution of the score. Every other quantity in this section — the average score, how much it wobbles, how far a 50-case eval can drift from the truth — is arithmetic on those three numbers.

import numpy as np
rng = np.random.default_rng(0)

values = np.array([0, 1, 2])          # what the grader can award
shares = np.array([0.2, 0.5, 0.3])    # how often each score happens; they sum to 1

draws = rng.choice(values, size=200_000, p=shares)
print(np.bincount(draws) / len(draws))
# [0.200215 0.50064  0.299145]

Expectation and variance, by hand then by counting#

The expectation — also called the mean, or the expected value — is the long-run average: each value multiplied by its share, added up. Here that is 0 times 0.2, plus 1 times 0.5, plus 2 times 0.3, which is 0 plus 0.5 plus 0.6, so 1.1. Notice what 1.1 is not. It is not a score any single answer can receive. An expectation does not have to be an attainable value, which is why an average of anything counted needs reading with care.

The variance measures spread: it is the expectation of the squared distance from the mean. The distance is squared so that two below the mean counts the same as two above, and so that one large miss counts for more than several small ones. By hand: a 0 is 1.1 away, squared 1.21, weighted by its share 0.2, contributing 0.242. A 1 is 0.1 away, squared 0.01, weighted 0.5, contributing 0.005. A 2 is 0.9 away, squared 0.81, weighted 0.3, contributing 0.243. Adding the three gives 0.49.

Squaring left the units wrong — 0.49 is in points squared — so the number usually quoted is its square root, 0.7, the standard deviation. That one is in points and can be compared directly with the mean of 1.1. If you remember one thing from the pair, it is that the mean tells you where a quantity sits and the standard deviation tells you how much a single observation of it is worth trusting.

mean = (values * shares).sum()
variance = (shares * (values - mean) ** 2).sum()
print(mean, variance, variance ** 0.5)
# 1.1 0.49 0.7

print(draws.mean(), draws.var())
# 1.09893 0.4895728551000001

Distributions#

• Normal: heights, measurement error, averages of many things. Symmetric around a mean.
• Binomial: counts of successes in trials. Pass rates on an eval set are binomial.
• Long-tailed: latency, request sizes, income. The mean is misleading and the median plus a high percentile is honest.

Conditional probability#

The probability of something given that something else is true. The definition is one line: the probability of A given B is the probability that A and B both happen, divided by the probability of B. In words, and this is the version to actually hold in your head: throw away every case where B did not happen, then ask what share of the cases left over had A. Conditioning is a filter, not a formula.

The place it goes wrong is that two facts which sound the same are different filters. Flip two fair coins. Four outcomes are equally likely: heads-heads, heads-tails, tails-heads, tails-tails. Told that at least one coin came up heads, you throw away tails-tails and three outcomes survive, of which one is two heads, so the answer is one in three. Told instead that the first coin came up heads, you keep heads-heads and heads-tails, so the answer is one in two. Same coins, same question, different filter, and most people give one half to both.

flips = rng.integers(0, 2, size=(200_000, 2))     # 1 is heads
both = flips.sum(axis=1) == 2
at_least_one = flips.sum(axis=1) >= 1
first_is_head = flips[:, 0] == 1

print(both[at_least_one].mean())     # 0.3345792524485309   -> one third
print(both[first_is_head].mean())    # 0.5008677438659486   -> one half

Independence is the case where the filter changes nothing. Two events are independent when the probability of both is the product of the two separate probabilities, which is the same as saying that learning one of them leaves your belief about the other exactly where it was. The first coin and the second coin are independent. The first coin and 'at least one head' are not, because knowing the first is heads guarantees the second fact.

second_is_head = flips[:, 1] == 1

print(first_is_head.mean() * second_is_head.mean(),
      (first_is_head & second_is_head).mean())
# 0.2507678055 0.251085     -> the product matches, so independent

print(first_is_head.mean() * at_least_one.mean(),
      (first_is_head & at_least_one).mean())
# 0.376200585 0.5013        -> nowhere near, so not independent

One more, because interviewers reach for it and because it breaks the most common wrong mental model. Uncorrelated does not mean independent. Correlation only detects a straight-line relationship. Take x spread evenly between minus one and one, and let y be x squared. Knowing x tells you y exactly — there is no relationship more dependent than that — and the correlation is zero.

x = rng.uniform(-1, 1, 200_000)
y = x ** 2

print(np.corrcoef(x, y)[0, 1])       # 0.0022397797187538398  -> effectively zero
print((y > 0.25).mean())             # 0.499755   unconditionally, a coin flip
print((y[x > 0.5] > 0.25).mean())    # 1.0        given x > 0.5, a certainty

Bayes' theorem#

Bayes' theorem reverses the direction of a conditional probability. You usually know how the evidence behaves given the state of the world — how often the test fires when someone is ill — and you want the state of the world given the evidence, which is the opposite. The rule is: the probability of A given B equals the probability of B given A, times the probability of A, divided by the probability of B. Read as a procedure it is even shorter. Take what you believed before the evidence, multiply by how likely that evidence would be if the belief were true, do the same for every competing belief, and rescale so the results add to one.

The example every interviewer uses is the medical test, so work it properly once. A test is 99 per cent accurate for a condition that affects 1 person in 1000. Spell out what 99 per cent accurate means, because the phrase hides two numbers: the test catches 99 per cent of the people who have the condition, and it correctly clears 99 per cent of the people who do not, so it wrongly flags 1 per cent of them. You test positive. What is the chance you have it?

A test is 99% accurate. The condition affects 1 in 1000.

Out of 100,000 people:
  100 have it        ->  99 test positive        (99% of 100)
  99,900 do not      -> 999 test positive        (1% of 99,900)

  1,098 positives in total, of which 99 are real.
  99 / 1,098 = 0.090

A positive result means about a 9% chance of actually having it.

The symbolic version gives the identical number. The probability of being ill given a positive is 0.99 times 0.001, divided by that same product plus 0.01 times 0.999. The numerator is 0.00099 and the denominator is 0.00099 plus 0.00999, which is 0.01098, so the answer is 0.0902. It is the counts argument with the 100,000 divided out.

def posterior(prior, sensitivity, specificity):
    hit = sensitivity * prior                       # ill, and flagged
    false_alarm = (1 - specificity) * (1 - prior)   # healthy, and flagged anyway
    return hit / (hit + false_alarm)

for prior in (0.001, 0.01, 0.10):
    print(prior, round(posterior(prior, 0.99, 0.99), 4))
# 0.001 0.0902
# 0.01 0.5
# 0.1 0.9167

Check it by simulation rather than believing the algebra, which takes one more block. Two million people, one in a thousand ill, everyone tested.

n = 2_000_000
has_it = rng.random(n) < 0.001
tests_positive = np.where(has_it, rng.random(n) < 0.99, rng.random(n) < 0.01)

print(tests_positive.sum(), has_it[tests_positive].mean())
# 21848 0.09044306114976199

There is a form of the rule worth memorising because it makes repeated evidence trivial. Work in odds rather than probabilities: prior odds times the likelihood ratio gives posterior odds, and the likelihood ratio for a positive is sensitivity divided by the false positive rate, here 0.99 over 0.01, which is 99. Prior odds of 1 to 999, multiplied by 99, gives 0.099 to 1, which is the 9 per cent again. Test a second time with an independent test and multiply by 99 once more: odds of 9.81 to 1, or 91 per cent.

odds = 0.001 / 0.999            # prior odds: 1 case per 999 healthy people
for test in (1, 2, 3):
    odds = odds * 99            # each independent positive multiplies by the ratio
    print(test, round(odds, 3), round(odds / (1 + odds), 4))
# 1 0.099 0.0902
# 2 9.811 0.9075
# 3 971.27 0.999

None of this is medical trivia. It is the arithmetic of every rare-event classifier you will ship. A fraud model at a 1 per cent base rate with a 1 per cent false positive rate sends an analyst one real case for every one false alarm, and that is the good version; at a 0.1 per cent base rate the analyst sees ten false alarms per real case and stops reading the queue. The quantity that decides whether an alerting system is usable is precision at the operating point, and precision is this calculation.

Sampling variation#

Measuring on 50 examples gives a number with real uncertainty in it. Two systems scoring 82 percent and 86 percent on 50 examples may well be identical, and the difference is noise.

import math

def margin(p, n):
    """Rough 95% margin of error for a proportion."""
    return 1.96 * math.sqrt(p * (1 - p) / n)

print(margin(0.85, 50), margin(0.85, 500), margin(0.85, 5000))
# 0.09897514839594836 0.031298690068435775 0.009897514839594836

That formula is worth understanding rather than memorising, and it comes straight out of the previous two sections. A pass on one eval case is a random variable worth 1 or 0; if the true pass rate is p its variance is p times one minus p. Average n of them and the variance divides by n, so the standard deviation of the average is the square root of p times one minus p over n. The 1.96 turns a standard deviation into a 95 per cent range, for reasons the central limit theorem supplies later in this topic.

Now put a number on the mistake in the headline. Take two systems that are genuinely identical, both passing 85 per cent of the time, and score each on its own 50 cases. How often does the pair differ by four points or more — enough for someone to write that one is better?

a = rng.binomial(50, 0.85, 200_000) / 50      # 200,000 evals of system A
b = rng.binomial(50, 0.85, 200_000) / 50      # and of an identical system B

print((np.abs(a - b) >= 0.04).mean())    # 0.624755
print((np.abs(a - b) >= 0.10).mean())    # 0.13899

This is the most useful thing on this page. Before claiming an improvement, check whether your eval set is large enough for the difference to mean anything.

Paired comparison is stronger#

The 62 per cent above is pessimistic on purpose: it assumed the two systems were scored on different samples. In practice you score both on the same examples, and that is a much better experiment, because the hard cases are hard for both and the easy cases are easy for both. All the shared difficulty cancels. So throw away the aggregate scores and look only at the cases where the two systems disagree.

Say A scores 41 out of 50 and B scores 43, the familiar 82 against 86. On 45 of the 50 cases they gave the same verdict. Of the five where they differed, B won three and A won two. If the two systems were equally good, each disagreement would go either way like a coin flip, so the question is whether 3 out of 5 is a surprising number of coin flips. It is not: the two-sided probability of a split at least that lopsided is 1.0, meaning every possible outcome is at least as extreme as this one.

from scipy import stats

print(stats.binomtest(3, 5, 0.5).pvalue)   # 1.0      the actual result: 3 wins of 5
print(stats.binomtest(5, 5, 0.5).pvalue)   # 0.0625   a clean sweep of all 5
print(stats.binomtest(6, 6, 0.5).pvalue)   # 0.03125  six disagreements, all one way

That reframes what to ask for. The number that decides whether a comparison can conclude anything is not the size of the eval set, it is the number of cases the two systems actually disagree on. Two prompts that differ on five of fifty cases cannot be separated no matter how the five fall. Either find more cases where they behave differently — which usually means harder cases — or accept that the change is too small for this eval to see.

Loss functions can look like arbitrary choices somebody made. They are not: nearly all of them come from one principle, which is to pick the parameters that make the data you actually observed most probable. This section works from there to why cross-entropy is the classification loss, then on to two practical tools: a way to put a range around any measurement without any distribution theory, and a way to check whether a model's stated confidence means anything. The second matters the moment you show a probability to a user.

Maximum likelihood#

Pick the parameters that make the observed data most probable. It sounds abstract and it is what nearly all training does. Minimising a loss function is usually maximising a likelihood written the other way round.

The word likelihood is doing precise work and is worth pinning down, because it is the same arithmetic as a probability read in the opposite direction. A probability answers: given this coin is 70 per cent heads, how often do I see seven heads in ten flips? A likelihood answers: given that I did see seven heads in ten flips, how well does 70 per cent explain it? The formula is identical. What changes is which part is fixed. The data is fixed and known; the parameter is the thing being varied.

So the procedure is: write the probability of the data you actually got as a function of the parameter, then find the parameter that makes it largest. Seven heads in ten flips has probability p to the seventh times one minus p cubed, up to a counting constant that does not depend on p and so cannot move the answer. Try p equals 0.5 and you get 0.000977. Try 0.7 and you get 0.002224, more than twice as good. Try 0.9 and it falls to 0.000478, because nine-tenths heads struggles to explain the three tails.

import numpy as np
rng = np.random.default_rng(0)

heads, flips = 7, 10

def likelihood(p):
    return p ** heads * (1 - p) ** (flips - heads)

for p in (0.3, 0.5, 0.6, 0.7, 0.8, 0.9):
    print(p, round(likelihood(p), 6))
# 0.3 7.5e-05
# 0.5 0.000977
# 0.6 0.001792
# 0.7 0.002224
# 0.8 0.001678
# 0.9 0.000478

grid = np.linspace(0.001, 0.999, 999)
print("peak at", round(grid[likelihood(grid).argmax()], 3))   # peak at 0.7

Now the result that connects it to everything you train. Suppose your data came from a normal distribution and you want its mean. The log-likelihood of a normal, with the terms that do not involve the mean dropped, is minus the sum of squared distances from the mean, divided by twice the variance. That is worth reading twice: maximising the likelihood of a normal is exactly minimising a sum of squared errors. Squared error is not a convention somebody chose for convenience. It is what maximum likelihood says to do if you believe your noise is normal.

Differentiate to find where it peaks. The derivative of the sum of squared distances with respect to the mean is minus twice the sum of the deviations, and setting that to zero says the deviations must cancel — which happens at the arithmetic average and nowhere else. So the sample mean is the maximum likelihood estimate of a normal's mean. Check it on five numbers: 4.1, 5.3, 4.8, 6.2 and 5.1 add to 25.5, so the mean is 5.1, and the log-likelihood should be highest there and symmetric either side of it.

x = np.array([4.1, 5.3, 4.8, 6.2, 5.1])

def log_likelihood(mu, sigma):
    return (-0.5 * np.log(2 * np.pi * sigma ** 2) - (x - mu) ** 2 / (2 * sigma ** 2)).sum()

for mu in (4.9, 5.0, 5.1, 5.2, 5.3):
    print(mu, round(log_likelihood(mu, 1.0), 4))
# 4.9 -5.8647
# 5.0 -5.7897
# 5.1 -5.7647     <- the peak
# 5.2 -5.7897
# 5.3 -5.8647

print(x.mean(), ((x - x.mean()) ** 2).sum())
# 5.1 2.3400000000000007

The same derivation gives the variance: divide that 2.34 by the five observations to get 0.468. And here is the wrinkle interviewers like, because it is the one place maximum likelihood is visibly imperfect. That estimate is biased low. It measures spread around the sample mean rather than the true mean, and the sample mean sits, by construction, in the middle of the very data being measured, so the distances come out too small. The correction is to divide by n minus one instead of n, which is where that famous and otherwise unmotivated minus one comes from.

mle, unbiased = [], []
for _ in range(20_000):
    sample = rng.normal(10.0, 2.0, 5)          # true variance is 4.0
    mle.append(((sample - sample.mean()) ** 2).mean())     # divide by n
    unbiased.append(sample.var(ddof=1))                    # divide by n - 1

print(round(np.mean(mle), 3), round(np.mean(unbiased), 3))
# 3.205 4.006

Where cross-entropy comes from#

For classification, the model outputs a probability for each class. The likelihood of the data is the product of the probabilities assigned to the correct classes. Taking the log turns the product into a sum, and negating it gives cross-entropy. It is not an arbitrary choice; it is the likelihood.

Correct class given probability 0.9   ->  loss -log(0.9)  = 0.11
Correct class given probability 0.5   ->  loss -log(0.5)  = 0.69
Correct class given probability 0.01  ->  loss -log(0.01) = 4.6

A confident wrong answer is punished very hard, which is
exactly the behaviour you want.

Trace it once through real numbers, because the chain from raw model output to loss is short and worth being able to draw from memory. The network emits one unbounded score per class, called a logit. The softmax exponentiates them and divides by the total, which forces them positive and makes them sum to one, so they can be read as probabilities. Cross-entropy then picks out the single probability assigned to the correct class and takes its negative logarithm. Three logits of 2.0, 1.0 and 0.1 become probabilities of 0.659, 0.242 and 0.099, so the loss is 0.417 if the answer was the first class and 2.317 if it was the third.

logits = np.array([2.0, 1.0, 0.1])          # raw model outputs, any sign, any size
e = np.exp(logits - logits.max())           # subtract the max: identical result, no overflow
probs = e / e.sum()
print(np.round(probs, 4))                   # [0.659  0.2424 0.0986]

print(round(-np.log(probs[0]), 4), round(-np.log(probs[2]), 4))
# 0.417 2.317

Once you see losses as likelihoods, the choice of loss stops being a menu and becomes a statement about the noise. Assume the target is the truth plus normal noise and you get squared error. Assume each label is a coin flip whose bias the model predicts and you get cross-entropy. Assume counts arriving at a rate the model predicts and you get the Poisson loss used for demand and traffic forecasting. Assume noise with heavier tails than normal and you get absolute error, which is why mean absolute error is the robust default on data with outliers. Picking a loss is picking a noise model, and if you cannot say what noise model your loss implies, you have not chosen it, you have inherited it.

Bootstrapping#

When you cannot derive a confidence interval, resample your data with replacement many times, recompute the statistic each time, and look at the spread. It is general, it needs no distributional assumptions, and it answers the question most people actually have about an eval score.

The idea behind it is a sleight of hand that works. You wanted many samples from the population and you only have one sample of fifty. So treat the fifty you have as if they were the population, draw fifty from them with replacement, and score that. Repeat ten thousand times. The spread of those ten thousand scores mimics the spread you would have seen across ten thousand genuine eval sets. With replacement is the part people get wrong and the part that makes it work: without it, every resample is the original set reordered and the spread is zero.

pass_fail = rng.binomial(1, 0.85, 50)      # 50 graded eval cases, 1 = passed
print(pass_fail.sum(), pass_fail.mean())   # 44 0.88

def bootstrap_ci(results, n=10_000):
    results = np.asarray(results)
    samples = rng.choice(results, (n, len(results)), replace=True)
    means = samples.mean(axis=1)
    return np.percentile(means, [2.5, 97.5])

print(bootstrap_ci(pass_fail))             # [0.78 0.96]

Calibration#

A calibrated model's 80 percent predictions are right about 80 percent of the time. Most modern networks are overconfident. Calibration matters whenever a probability drives a decision, such as a threshold or a confidence display, and it can be measured by grouping predictions by confidence and checking accuracy in each group.

Measuring it takes about ten lines and no theory. Sort every prediction into buckets by the probability it stated, then within each bucket compare the average stated probability against the share that actually came true. A calibrated model produces two columns of matching numbers. The average gap, weighted by how many predictions fall in each bucket, is the expected calibration error, and it is the number usually reported.

n = 20_000
truth_rate = rng.uniform(0.05, 0.95, n)              # the real chance for each case
happened = (rng.random(n) < truth_rate).astype(int)  # what actually occurred

logit = np.log(truth_rate / (1 - truth_rate))
stated = 1 / (1 + np.exp(-2.0 * logit))              # a model that doubles its logits

edges = np.linspace(0, 1, 6)
bucket = np.clip(np.digitize(stated, edges) - 1, 0, 4)
ece = 0.0
for b in range(5):
    m = bucket == b
    says, actually = stated[m].mean(), happened[m].mean()
    ece += m.mean() * abs(says - actually)
    print(f"{edges[b]:.1f}-{edges[b+1]:.1f}  n={m.sum():5d}  says {says:.2f}  happens {actually:.2f}")
print("ECE", round(ece, 3))

# 0.0-0.2  n= 6303  says 0.07  happens 0.19
# 0.2-0.4  n= 2641  says 0.29  happens 0.39
# 0.4-0.6  n= 2224  says 0.50  happens 0.50
# 0.6-0.8  n= 2575  says 0.70  happens 0.60
# 0.8-1.0  n= 6257  says 0.93  happens 0.81
# ECE 0.102

Four numbers turn up constantly in model reporting and are routinely confused with one another. This section separates them, starting from the single quantity all four are averages of, and gives each a plain reading: what the floor is, what your model costs above that floor, what the gap between the two is used for, and what a perplexity figure actually tells you. It also covers when two of these numbers cannot be compared at all, which is the mistake most often made with them in public.

Surprise#

The surprise of an outcome is the negative logarithm of the probability you gave it. Something you said was certain surprises you not at all; something you gave one chance in a thousand surprises you a great deal. Everything on this page is an average of that one quantity, taken over a different set of assumptions.

surprise of an outcome     -log p
entropy        H(p)        average surprise, using the true distribution for both
cross-entropy  H(p, q)     average surprise of your q, when reality is p
KL             KL(p || q)  = H(p, q) - H(p)    the excess you pay for using q

log base 2 gives bits, natural log gives nats. Keep one base throughout.

Entropy is how uncertain the truth itself is: a fair coin has one bit, a loaded one has less, a certain outcome has none. It is the floor. No model, however good, can average less surprise than the entropy of what it is predicting, which is why a training loss that stops falling is not necessarily a model that stopped learning.

Compute it once by hand so the units mean something. A fair coin: each side has probability a half, whose negative log base two is one, and the average of one and one is one bit. A coin that lands heads nine times in ten: heads carries a surprise of 0.152 bits and tails carries 3.322, and the average weighted by how often each happens is 0.9 times 0.152 plus 0.1 times 3.322, which is 0.469 bits. Less than half the fair coin, because you can usually guess. Push it to 99 per cent heads and it falls to 0.081 bits, and at certainty it is zero.

import math

def entropy_bits(p):
    if p in (0.0, 1.0):
        return 0.0
    return -(p * math.log2(p) + (1 - p) * math.log2(1 - p))

for p in (0.5, 0.75, 0.9, 0.99, 1.0):
    print(p, round(entropy_bits(p), 4))
# 0.5 1.0
# 0.75 0.8113
# 0.9 0.469
# 0.99 0.0808
# 1.0 0.0

Cross-entropy and KL#

Cross-entropy is what you actually pay: the average surprise of your predicted distribution against outcomes drawn from the real one. It is the training loss for classification and for next-token prediction. KL divergence is the difference between what you paid and the floor, so it is the part of the loss that is your model's fault rather than the world's.

That split has a direct consequence for how you read a loss number, and it is the reason two training losses are so often compared when they should not be. A loss of 1.6 nats is excellent on text whose own entropy is 1.5 and dreadful on text whose entropy is 0.5, and nothing in the number itself distinguishes the two cases. So a loss curve is only comparable against another run on the same data. A loss that flattens at 1.6 may be a model that has run out of capacity, or a model that has learned everything there is and is now sitting on the floor, and the way to tell them apart is to change the model rather than to stare at the curve.

import math

def kl(p, q):
    return sum(pi * math.log(pi / qi) for pi, qi in zip(p, q))

p = [0.7, 0.2, 0.1]
q = [0.5, 0.3, 0.2]

kl(p, q)     # 0.085 nats
kl(q, p)     # 0.092 nats  -- a different number from the same pair

The asymmetry is not a technicality. KL(p || q) punishes q for putting almost no probability where p has some, so minimising it spreads q out to cover everything. KL(q || p) punishes q for putting probability where p has none, so minimising it makes q concentrate on one mode and ignore the rest. Which direction a method uses decides whether its approximation hedges or commits.

• Preference training uses KL from the starting model as a leash, penalising the tuned model for drifting too far from what it began as.
• Distillation trains the student by minimising KL against the teacher's full output distribution, which carries more information than the label alone.
• Monitoring uses it to detect drift, comparing the distribution of today's inputs against the distribution the model was validated on.

Perplexity#

Perplexity is the exponential of cross-entropy, taken in the same base. It exists because a loss of 1.9 nats per token means nothing to anyone, while its exponential does: a perplexity of 8 says the model is as uncertain, on average, as if it were choosing uniformly between 8 equally likely tokens at each step. That is the average branching factor, and it is the only intuition the number carries.

# Natural-log probabilities the model gave to each correct next token.
logprobs = [-0.10, -2.30, -0.69, -1.61, -0.36]

cross_entropy = -sum(logprobs) / len(logprobs)
perplexity = math.exp(cross_entropy)
print(round(cross_entropy, 3), round(perplexity, 3))     # 1.012 2.751

# Comparable across tokenizers, because bytes are not a modelling choice:
n_bytes = 23
print(round(-sum(logprobs) / math.log(2) / n_bytes, 3))  # 0.317

The branching-factor reading gives entropy a job at inference time as well as in training. Apply the same average-surprise calculation to the model's next-token distribution and its exponential tells you, in effect, how many tokens the model is really choosing between at that step. That number is what a sampling temperature moves. Dividing the logits by a temperature below one sharpens the distribution and lowers its entropy; dividing by a temperature above one flattens it and raises the entropy.

import numpy as np

logits = np.array([4.0, 3.0, 1.0, 0.5, 0.0])
for T in (0.5, 1.0, 2.0):
    z = logits / T
    e = np.exp(z - z.max())
    probs = e / e.sum()
    H = -(probs * np.log(probs)).sum()                 # entropy in nats
    print(T, np.round(probs, 3), round(float(H), 3), round(float(np.exp(H)), 2))

# 0.5 [0.878 0.119 0.002 0.001 0.   ] 0.389 1.48
# 1.0 [0.682 0.251 0.034 0.021 0.012] 0.857 2.36
# 2.0 [0.468 0.284 0.104 0.081 0.063] 1.327 3.77

It is also not a task metric. Perplexity measures how well the model predicts held-out text of the same kind as its training data, which is a different question from whether it does your job. It is a good early warning, a good check on whether quantization damaged something, and a poor argument for shipping.

Probability questions are easy to get wrong by reasoning about them and hard to get wrong by counting. This section takes that seriously: for each idea you simulate the situation tens of thousands of times in a few lines of Python and count what happened, then do the maths and check the two agree. The habit is worth more than any single result here, because it means you can settle an argument about a probability on your own instead of looking for someone to confirm it.

Count what happens#

Every simulation in this section has the same three steps, and they are worth naming because the third is the one people leave out. Generate the situation many times. Reduce each repetition to a true or false. Take the mean of those booleans, which in NumPy is the share that were true. That is the probability. Two dice summing to seven is the smallest example: there are 36 equally likely pairs and six of them sum to seven, so the answer is one sixth, and the simulation should land on 0.167.

import numpy as np
rng = np.random.default_rng(0)

rolls = rng.integers(1, 7, size=(100_000, 2))       # two dice, 100,000 times
print((rolls.sum(axis=1) == 7).mean())              # 0.16775, against 6/36 = 0.1667

Any probability is the long-run share of times something happens. If you can simulate the situation, you can estimate it by counting, and check an exact answer you worked out on paper.

How wrong can a simulation be?#

A simulated answer is itself a measurement, so it has its own error bar, and knowing it is what separates using this technique from waving at it. The estimate is a proportion from n independent trials, so the margin is the same one as for any proportion: 1.96 times the square root of p times one minus p over n. At 100,000 trials and a true answer near one sixth, that is about two parts in a thousand, which is why 0.16775 and 0.1667 agreeing to two decimals is a pass and agreeing to four would be luck.

import math

def mc_margin(p_hat, trials):
    return 1.96 * math.sqrt(p_hat * (1 - p_hat) / trials)

for trials in (1_000, 100_000, 10_000_000):
    print(trials, round(mc_margin(0.1667, trials), 5))
# 1000 0.0231
# 100000 0.00231
# 10000000 0.00023

More samples, closer answer#

That square-root law is the law of large numbers with a rate attached. The average of repeated draws converges on the expectation, but slowly and not monotonically: it can be further away after more draws than it was before. Watch it wander in on a fair coin.

for n in [10, 100, 1_000, 100_000]:
    print(n, rng.integers(0, 2, n).mean())   # coin flips: drifts towards 0.5
# 10 0.8
# 100 0.49
# 1000 0.509
# 100000 0.49946

The birthday problem#

Twenty-three people in a room, and the chance that two share a birthday is better than even. Nearly everyone guesses far lower, because they answer a different question — the chance that someone shares their birthday — when the real one counts all 253 pairs. It is the standard demonstration that intuition about compound events is not to be trusted, and it has an exact answer to check against: multiply 365/365 by 364/365 by 363/365 and so on, which is the probability that everyone is different, then subtract from one.

def shared_birthday(people, trials=50_000):
    days = rng.integers(0, 365, size=(trials, people))
    s = np.sort(days, axis=1)
    return (np.diff(s, axis=1) == 0).any(axis=1).mean()

def exact(people):
    none_shared = 1.0
    for i in range(people):
        none_shared *= (365 - i) / 365
    return 1 - none_shared

for k in (10, 23, 40, 60):
    print(k, round(shared_birthday(k), 3), round(exact(k), 3))
# 10 0.116 0.117
# 23 0.509 0.507
# 40 0.893 0.891
# 60 0.994 0.994

The engineering version of this is hash collisions, and the numbers are the same shape. With a 64-bit hash you expect a collision after about four billion keys rather than after eighteen quintillion, because it is pairs that matter, not items. Any time you reason about the chance that some pair of things collides, the birthday calculation is the one to reach for.

Monty Hall#

Three doors, a car behind one, you pick a door, the host — who knows where the car is — opens a different door showing a goat and offers you the switch. Almost everyone says it makes no difference. It doubles your chances, and the simulation is shorter than the argument. The key is that your first pick is right one time in three, and the host's behaviour cannot change that, so the other unopened door holds the car the remaining two times in three. Switching wins exactly when your first pick was wrong.

trials = 200_000
car = rng.integers(0, 3, trials)
first_pick = rng.integers(0, 3, trials)

print("stay  ", (first_pick == car).mean())     # 0.332165
print("switch", (first_pick != car).mean())     # 0.667835

Base rates#

A test catches 99% of people with a condition and wrongly flags 5% of people without it. The condition affects 1% of people. You test positive. How likely is it you have it?

n = 1_000_000
sick = rng.random(n) < 0.01
positive = np.where(sick, rng.random(n) < 0.99, rng.random(n) < 0.05)
print(sick[positive].mean())    # 0.16511298003850441: most positives are false alarms

Because healthy people vastly outnumber sick ones, 5% of them outnumber 99% of the sick. That is Bayes' theorem done in counts rather than symbols, and the same arithmetic explains why a fraud model with a low false positive rate still produces mostly false alerts.

Getting there by simulation rather than by formula is not a lesser route. The three lines above encode the story directly — who is sick, what the test does to each group — so a misreading of the question shows up as code you cannot write, whereas a misreading of the formula shows up as a confident wrong number. When a conditional probability argument matters, write both and make them agree.

Averages become bell curves#

Waiting times are about as un-bell-shaped as a quantity gets: most waits are short, a few are very long, and the histogram is a cliff falling to the right. Take thirty of them and average, though, and something strange happens to the average. Its distribution is close to symmetric, close to normal, and much narrower than the waits themselves.

waits = rng.exponential(scale=2.0, size=(100_000, 30))   # skewed, not bell-shaped
means = waits.mean(axis=1)                                 # the average of 30 waits
print(means.mean(), means.std())
# 2.00171649553826 0.3669219319476103
# the mean of the averages is the mean of the waits; 2.0 / sqrt(30) = 0.365

The central limit theorem is the claim that this happens for almost anything you average, whatever shape it started as, and that the limit is specifically the normal distribution. Claims like that deserve measurement rather than assent, so measure it. Standardise the average — subtract its mean, multiply by the square root of n — and a normal result must have its middle 95 per cent between minus 1.96 and plus 1.96, and zero skewness. Then watch those three numbers as n grows.

for n in (1, 2, 5, 30, 200):
    m = rng.exponential(1.0, size=(200_000, n)).mean(axis=1)
    z = (m - 1.0) * math.sqrt(n)          # standardised: normal means 2.5% at -1.96
    print(f"n={n:4d}  2.5% {np.percentile(z, 2.5):6.2f}   97.5% {np.percentile(z, 97.5):5.2f}"
          f"   skew {float(((z - z.mean()) ** 3).mean()):5.2f}")

# n=   1  2.5%  -0.97   97.5%  2.68   skew  2.00
# n=   2  2.5%  -1.24   97.5%  2.51   skew  1.36
# n=   5  2.5%  -1.51   97.5%  2.35   skew  0.90
# n=  30  2.5%  -1.78   97.5%  2.12   skew  0.37
# n= 200  2.5%  -1.89   97.5%  2.03   skew  0.14
# normal reference:  -1.96, 1.96, skew 0

That last point is the practical one, and it is why this is worth measuring rather than quoting. The average of many independent values is approximately normally distributed whatever the shape of the original values, with spread shrinking as the square root of the sample size, which is why error bars on averages work. But approximately is doing real work in that sentence. The rule of thumb that thirty is enough comes from mildly skewed data; on revenue per user, or latency, or anything where one observation in a thousand is a hundred times the others, thirty is nowhere near enough and the error bar you computed is too narrow. Simulate your own data before trusting a normal interval on it.

A distribution describes how a quantity tends to come out when it varies: request latencies, support tickets per hour, measurement error. Knowing which one fits a quantity tells you what counts as normal, what is genuinely unusual, and which summary numbers will mislead you. This section is a field guide to the few you actually meet, with the code to sample and fit each, and a warning about the heavy-tailed ones, where the average is not a useful number at all.

Two numbers summarise any of them#

Whatever the shape, two numbers describe it. The mean is the centre. The variance is the spread: take each value's distance from the mean, square it so that being below counts the same as being above, and average those. Squaring leaves the units wrong — a variance of latencies is in milliseconds squared — so the number usually quoted is its square root, the standard deviation, which is back in the original units and can be compared with the mean.

Both are expectations, which is the word for a probability-weighted average: the mean is the expectation of the value, and the variance is the expectation of the squared distance from the mean. When you have a distribution written down, that weighting is by each outcome's probability. When you have a pile of data instead, every observation is equally likely, so the weighting collapses into the ordinary average and the two definitions become the same arithmetic.

import numpy as np

latencies = np.array([90, 95, 100, 105, 110])
mean = latencies.mean()
variance = ((latencies - mean) ** 2).mean()
print(mean, variance, np.sqrt(variance))
# 100.0 50.0 7.0710678118654755

# Same mean, far more spread. The mean alone tells you nothing about this.
spiky = np.array([10, 20, 100, 180, 190])
print(spiky.mean(), spiky.std())
# 100.0 76.15773105863909

A field guide#

Quantity                                   Distribution to try first
Did the user convert? (yes/no)             Bernoulli; conversions out of n: binomial
Support tickets per hour                   Poisson
Time until the next request                exponential
Measurement error, heights, averages       normal
Request latency, file sizes, incomes       log-normal or other heavy-tailed
Which of k classes                         categorical

The rest of this section takes that table one row at a time. For each distribution: the signal in the problem statement that should make you reach for it, the two numbers that describe it, and a check in code. Recognition is the skill being trained here — the maths of each is a lookup, but noticing that a quantity is a count per period rather than a proportion is not.

Yes or no: Bernoulli and binomial#

The recognition signal is a single trial with two outcomes. Did the user convert, did the request succeed, did the eval case pass. That is a Bernoulli variable: it takes the value 1 with probability p and 0 otherwise. Its mean is p, because the average of a pile of ones and zeros is the share of ones. Its variance is p times one minus p, which is largest at p equals a half and shrinks to nothing at either end — a system that succeeds 99 per cent of the time varies far less from sample to sample than one that succeeds half the time.

The binomial is the same trial repeated n times, counting successes. Its mean is n times p and its variance is n times p times one minus p, which is just n independent Bernoulli variances added. Any pass rate on an eval set of fixed size is binomial, which is why the margin-of-error formula for a proportion works and where it comes from.

rng = np.random.default_rng(0)

one = rng.binomial(1, 0.3, 200_000)              # Bernoulli(0.3)
print(one.mean(), one.var())                     # p, then p(1 - p) = 0.21
# 0.299145 0.209657268975

k = rng.binomial(20, 0.3, 200_000)               # Binomial(20, 0.3)
print(k.mean(), k.var())                         # np = 6, then np(1 - p) = 4.2
# 6.00427 4.2180517671

Counts per period: Poisson#

The recognition signal is a count of events in a window, where there is no natural maximum. Support tickets per hour, errors per minute, arrivals per second. The binomial needs an n; the Poisson does not, which is exactly what makes it the right shape for counts. It has one parameter, the average rate, and its defining oddity is that its variance equals its mean. If you measure a count whose variance is much larger than its mean, the events are not arriving independently — they are bursty, and a Poisson model will badly understate your worst hours.

The Poisson and the exponential are the same story told twice. If events arrive independently at a steady rate, the count per window is Poisson and the gaps between consecutive events are exponential. You can derive one from the other in code rather than on paper: generate exponential gaps, add them up, and count how many land in the first hour.

lam = 12.0                                             # 12 tickets an hour on average
gaps = rng.exponential(1 / lam, size=(50_000, 40))     # exponential time between arrivals
arrival = np.cumsum(gaps, axis=1)
per_hour = (arrival <= 1.0).sum(axis=1)

print(per_hour.mean(), per_hour.var())     # a Poisson has variance equal to its mean
# 11.9993 11.954779509999998

from scipy import stats
print((per_hour == 20).mean(), stats.poisson(mu=12).pmf(20))
# 0.00976 0.00968203216821916

Waiting times: exponential#

The recognition signal is time until the next thing happens, when the thing has no memory of how long you have already waited. Its mean is one over the rate, so at twelve arrivals an hour the average gap is five minutes, and its standard deviation equals its mean, which is why waiting times feel so unpredictable.

Memorylessness is the property worth checking, because it is the one that violates intuition. Having already waited three minutes tells you nothing about how much longer you will wait: the distribution of remaining time is identical to the distribution you faced at the start. This is why a retry after a timeout on an exponential service does not improve your position, and it is exactly the assumption that fails for real queues, where a request that has already been slow is usually stuck on something and will keep being slow.

w = rng.exponential(2.0, 500_000)        # mean wait of 2

print((w > 2).mean())                    # 0.367862   chance of waiting past 2
print((w[w > 3] > 5).mean())             # 0.36882405745062835
# given you have already waited 3, the chance of another 2 is the same number

Measurements: normal#

The recognition signal is a quantity built from many small independent contributions: heights, measurement error, and above all averages of anything at all, by the central limit theorem. It takes two parameters, the mean and the standard deviation, and those two describe it completely — there is no third number to learn. The figures worth memorising are that 68 per cent of the mass lies within one standard deviation of the mean, 95 per cent within two, and 99.7 per cent within three, which is where the 1.96 in every confidence interval comes from.

z = rng.normal(0, 1, 500_000)
for k_ in (1, 2, 3):
    print(k_, (np.abs(z) < k_).mean())
# 1 0.682944
# 2 0.954696
# 3 0.997334

Which of k: categorical, and the softmax#

The recognition signal is a choice among k options with no ordering, and it is the distribution you will meet most often in this field, because it is what a language model emits at every step. A categorical distribution is just k probabilities that sum to one. Models do not produce those directly; they produce k unbounded scores, and the softmax turns scores into a categorical distribution by exponentiating and dividing by the total.

Sampling from it is then one call. The temperature knob divides the scores before the exponential, so a temperature below one exaggerates the gaps and concentrates the mass on the leader, while a temperature above one flattens everything towards uniform. Watching the sampled frequencies match the computed probabilities is the check that the sampler is doing what you think.

logits = np.array([4.0, 3.0, 1.0, 0.5, 0.0])       # raw model scores over 5 tokens

def softmax(z, T=1.0):
    e = np.exp(z / T - (z / T).max())
    return e / e.sum()

for T in (0.5, 1.0, 2.0):
    p = softmax(logits, T)
    draws = rng.choice(len(logits), size=100_000, p=p)
    print(T, np.round(p, 3), np.round(np.bincount(draws, minlength=5) / 100_000, 3))

# 0.5 [0.878 0.119 0.002 0.001 0.   ] [0.878 0.118 0.002 0.001 0.   ]
# 1.0 [0.682 0.251 0.034 0.021 0.012] [0.682 0.251 0.034 0.021 0.013]
# 2.0 [0.468 0.284 0.104 0.081 0.063] [0.467 0.285 0.104 0.082 0.063]

With SciPy#

SciPy gives every one of these the same four methods, which is worth knowing because it means learning one distribution object teaches you all of them. The pmf or pdf is the height at a value, the cdf is the mass at or below a value, ppf is the inverse of the cdf and is how you get a percentile, and rvs draws samples. Nearly every question you have about a distribution is one of those four.

print(stats.binom(n=1000, p=0.03).interval(0.95))   # (20.0, 41.0) conversions to expect
print(stats.poisson(mu=12).pmf(20))                 # 0.00968203216821916
print(stats.norm(loc=0, scale=1).cdf(1.96))         # 0.9750021048517795

Heavy tails break averages#

Latency is not normally distributed. Most requests are fast and a few are extremely slow. The mean is dragged up by the slow ones and describes almost nobody's experience. That is why services report the median and high percentiles such as p95 and p99.

lat = rng.lognormal(mean=5, sigma=1, size=100_000)       # a realistic latency shape
print(f"mean {lat.mean():.0f}  median {np.median(lat):.0f}  p99 {np.percentile(lat, 99):.0f}")
# mean 245  median 147  p99 1527

# What a normal model would have predicted for the p99:
print(round(float(lat.mean() + 2.326 * lat.std())))
# 996

Fitting the right family instead takes one line and gives back the parameters you generated with. stats.lognorm.fit(lat, floc=0) returns a shape of 1.003 and a log scale of 4.999, against the 1.0 and 5.0 used to make the data. The reason to fit rather than to eyeball is that the fitted object then answers questions the sample cannot: the chance of exceeding a timeout you have never yet hit, or the percentile a new SLA would land on.

You have measured a difference between two things — two models, two prompts, two versions of a page. This section is about the question that follows: could a gap that size have turned up by chance alone? It covers what a p-value does and does not say, the several ways it is routinely misread, how to run the standard tests in SciPy, and why a confidence interval is usually the more honest thing to put in the report.

The definition, carefully#

A p-value is the probability of seeing data at least as extreme as yours if there were really no effect. A small p-value means your data would be surprising in a world with no effect. That is all it means.

• It is not the probability that there is no effect.
• It is not the probability your result is a fluke.
• It says nothing about how big or important the effect is. A tiny effect in a huge sample can have a tiny p-value.

The definition has a consequence that makes the whole thing click, and it is easier to see by simulation than to argue. If there is genuinely no effect, the p-value is a uniformly distributed random number between 0 and 1. Not usually large — uniform. So 5 per cent of the time it lands below 0.05, half the time below 0.5, and there is nothing special about any particular value it takes. Run two identical arms against each other twenty thousand times and look.

import math
import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

def two_proportion_test(x1, n1, x2, n2):
    """Successes and trials for each arm. Returns the two-sided p-value."""
    pooled = (x1 + x2) / (n1 + n2)
    se = math.sqrt(pooled * (1 - pooled) * (1 / n1 + 1 / n2))
    z = (x2 / n2 - x1 / n1) / se
    return 2 * (1 - stats.norm.cdf(abs(z)))

null_p = np.array([
    two_proportion_test(rng.binomial(2_000, 0.03), 2_000,
                        rng.binomial(2_000, 0.03), 2_000)
    for _ in range(20_000)
])
print((null_p < 0.05).mean(), (null_p < 0.10).mean(), (null_p < 0.50).mean())
# 0.049 0.09855 0.49675

A worked A/B test#

Take a real-shaped question. The control page converts at 3.0 per cent and the variant at 3.4 per cent, a relative lift of 13 per cent, on 10,000 visitors each. Someone wants to ship it on Friday. The two-proportion test gives p equal to 0.108, so by the usual convention it does not clear the bar — and the more useful output is the interval on the difference, which runs from minus 0.09 to plus 0.89 percentage points. The data is consistent with the variant being slightly worse and with it being a 30 per cent relative improvement. Run it to 40,000 per arm and the same rates give p equal to 0.0013 and an interval of 0.16 to 0.64 points.

def ab(x1, n1, x2, n2):
    p1, p2 = x1 / n1, x2 / n2
    se = math.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    lift = p2 - p1
    return p1, p2, two_proportion_test(x1, n1, x2, n2), lift - 1.96 * se, lift + 1.96 * se

for n in (10_000, 40_000):
    p1, p2, pv, lo, hi = ab(int(0.030 * n), n, int(0.034 * n), n)
    print(n, p1, p2, round(pv, 4), round(lo * 100, 3), round(hi * 100, 3))

# 10000 0.03 0.034 0.108  -0.088  0.888
# 40000 0.03 0.034 0.0013  0.156  0.644

Two common tests#

Proportions are the common case in this field, but two others turn up constantly. Use a t-test when you are comparing averages of a continuous quantity, and pass equal_var=False, because the equal-variance version assumes something you rarely have any reason to believe. Use a chi-square test on a contingency table when you are asking whether two categorical variables are related.

group_a = rng.lognormal(mean=3.0, sigma=0.6, size=400)     # session length, minutes
group_b = rng.lognormal(mean=3.1, sigma=0.6, size=400)
t, p = stats.ttest_ind(group_a, group_b, equal_var=False)
print(round(group_a.mean(), 1), round(group_b.mean(), 1), round(float(t), 2), round(float(p), 4))
# 23.9 26.8 -2.39 0.0169

table = [[120, 880],    # plan A: churned, stayed
         [ 90, 910]]    # plan B
chi2, p2, dof, expected = stats.chi2_contingency(table)
print(round(float(chi2), 2), round(float(p2), 4), dof)
print(np.round(expected, 1).tolist())
# 4.47 0.0344 1
# [[105.0, 895.0], [105.0, 895.0]]

Prefer an interval#

A confidence interval gives a range of effect sizes consistent with the data. 'Plan B churn is 3 points lower, 95% interval 0.2 to 5.8' tells you the direction, the size and the uncertainty at once, which a p-value alone does not.

It also makes the two honest failure modes distinguishable, and a p-value does not. An interval from minus 0.1 to plus 0.9 points says you have not learned much and should keep the experiment running. An interval from minus 0.05 to plus 0.06 points says you have learned something definite: whatever the effect is, it is too small to care about. Both of those come back as not significant, and they call for opposite decisions.

Power#

Power is the chance your study detects an effect of a given size if it exists. Small studies have low power, so they miss real effects, and the effects they do report tend to be exaggerated. Estimate the sample you need before collecting data.

Power is simulable, which means you never have to trust a formula you half remember. Decide the effect you would want to catch, generate data with that effect in it, run your actual test, and count how often it fires. Here the effect is the 3.0 against 3.4 per cent from above, and the question is how many visitors per arm it takes to see it reliably.

for n in (5_000, 20_000, 80_000):
    wins = sum(
        two_proportion_test(rng.binomial(n, 0.030), n, rng.binomial(n, 0.034), n) < 0.05
        for _ in range(4_000)
    )
    print(n, wins / 4_000)
# 5000 0.2125
# 20000 0.61925
# 80000 0.99425

The closed form agrees and is quicker for planning. For two proportions, the sample per arm is 2.8 squared — that constant being 1.96 for the false positive rate plus 0.84 for 80 per cent power — times the sum of the two variances, divided by the squared difference. For these numbers it says 30,353 per arm, and simulating at exactly that n gives 80.2 per cent power.

p1, p2 = 0.030, 0.034
n_needed = (1.96 + 0.84) ** 2 * (p1 * (1 - p1) + p2 * (1 - p2)) / (p2 - p1) ** 2
print(round(n_needed))       # 30353

wins = sum(
    two_proportion_test(rng.binomial(30_353, p1), 30_353, rng.binomial(30_353, p2), 30_353) < 0.05
    for _ in range(4_000)
)
print(wins / 4_000)          # 0.80225

Twenty metrics, one winner#

The most common way a team fools itself is not a bad test, it is many tests. Every dashboard has twenty metrics, and each one gets a 5 per cent chance of a false alarm, so the chance that at least one fires when nothing at all has changed is one minus 0.95 to the twentieth, which is 64 per cent. On average one metric will be significant every time you look.

print(round(1 - 0.95 ** 20, 4))          # 0.6415

flags = rng.random((50_000, 20)) < 0.05   # 20 metrics, nothing happening
print(flags.any(axis=1).mean(), flags.sum(axis=1).mean())
# 0.64172 0.99986