Estimation, likelihood and uncertainty

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