Statistical Inference for Machine Learning

You fit a model on the data you happen to have. You then use it on people, transactions or images you have never seen. Statistical inference is the set of tools for saying what the first thing licenses you to claim about the second.

Machine learning uses this field’s estimators constantly, while performance figures are often reported without any statement of uncertainty. A model reports 91.2% accuracy and the number gets written down as though it were measured with a ruler. It was not; it was estimated from one sample, and a second sample would have produced a different figure.

The number you compute is not the number you want

Suppose we want the average height of adults in a country. Measuring everyone is out of the question, so we measure twenty people and take their average. Run that study five times, with a different twenty people each time.

import numpy as np

rng = np.random.default_rng(0)
true_average = 170.0                    # the real average height, in cm

for study in range(5):
    sample = rng.normal(true_average, 8.0, size=20)
    print(f"  study {study + 1}: measured average {sample.mean():6.2f} cm")
# study 1: measured average 168.53 cm
# study 2: measured average 170.50 cm
# study 3: measured average 172.82 cm
# study 4: measured average 172.00 cm
# study 5: measured average 169.39 cm

Five honest studies of the same population, five different answers, spread over four centimeters. Nobody made a mistake. The truth was 170 the whole time and no study reported it.

The average of your sample is a fact about that sample, not about the population it came from; it is a roll of dice whose outcome depends on which twenty people walked in. Ask what its typical value is, and how far it usually lands from the truth, and you are doing statistical inference.

Three words get used for three different things, and they are worth separating before going further. The parameter is the number you want and cannot see — 170 here. The estimator is the recipe you apply to data — “average the sample”. The estimate is what the recipe returned this time — 168.53. Only the estimate is in front of you; only the parameter is what you care about.

Two different ways of being wrong

Since the estimate moves around, we should ask how it moves. There are two separate failures, and they need separate names because the fixes are different.

Compare three recipes on the same problem. The first averages the twenty measurements. The second does the same but with a tape measure that reads 3 cm too long. The third measures one person and goes home.

import numpy as np

rng = np.random.default_rng(1)
truth = 170.0

def sample_average(s): return s.mean()
def biased_ruler(s):   return s.mean() + 3.0      # every measurement 3 cm too tall
def first_person(s):   return s[0]                # measure one person, go home

for name, method in (("sample average", sample_average),
                     ("biased ruler  ", biased_ruler),
                     ("first person  ", first_person)):
    answers = np.array([method(rng.normal(truth, 8.0, size=20)) for _ in range(20_000)])
    print(f"  {name}  average answer {answers.mean():7.3f}"
          f"   spread {answers.std():5.3f}   typical miss {np.abs(answers - truth).mean():5.3f}")
# sample average  average answer 169.978   spread 1.787   typical miss 1.426
# biased ruler    average answer 173.007   spread 1.809   typical miss 3.080
# first person    average answer 170.048   spread 8.002   typical miss 6.372

The columns measure the two failure modes. Average answer is the mean of the 20,000 estimates; subtracting 170 gives an estimate of the method’s bias. Spread is their standard deviation — the square root of their variance — expressed in centimeters. Typical miss is their mean absolute error: the average absolute difference from 170 across the runs. Read the middle two columns against each other. The biased ruler is reliable — its spread of 1.809 is as tight as the good method’s 1.787 — and it is reliably wrong, averaging 173 when the truth is 170. Repetition alone does not remove the systematic offset — and it is not that every run overshoots: about 4.7% of these runs still land below 170. Without an external reference, a tightly clustered set of estimates can conceal the bias entirely.

The one-person method fails the opposite way. Its average answer is 170.048, essentially perfect — average a million of its estimates and that average will typically be close to the truth. But any single run is off by 6.4 cm on average, because it inherits the spread of one human being instead of averaging it away.

Those two failures have standard names. Bias is the gap between where the estimator lands on average and the truth, written \(\mathbb{E}[\hat{\theta}]-\theta\). Variance is how much it scatters around its own average. They combine into one summary of typical error, the mean squared error. It averages squared errors rather than absolute ones, so it is a different measure from the “typical miss” column above:

\[\text{MSE}=\text{bias}^{2}+\text{variance}\]

Both terms count, and that has a consequence people find surprising: a biased estimator can beat an unbiased one. If accepting a small systematic error buys a large reduction in scatter, the trade is worth making. The first exercise constructs exactly such a case, and it is the same trade ridge regression makes on purpose.

One more property is worth naming. An estimator is consistent if it closes in on the truth as the sample grows. Consistency concerns what happens as the sample grows; unbiasedness concerns the estimator’s expectation at a given sample size. Neither implies the other, and the one-person method is the counterexample in one direction: it is unbiased at every \(n\) and never gets closer, because it ignores the extra data.

Why the variance formula divides by n minus one

Here is a case where the bias is real, small, and visible. To estimate spread you take each measurement’s distance from the center, square it, and average. The question is what to divide by: the number of measurements, or one less.

import numpy as np

rng = np.random.default_rng(0)
for n in (2, 5, 20, 100):
    s = rng.normal(0, 1, size=(200_000, n))
    centered = s - s.mean(1, keepdims=True)
    mle = (centered ** 2).mean(1)                       # divide by n
    unbiased = (centered ** 2).sum(1) / (n - 1)         # divide by n-1
    print(f"n={n:4d}  E[MLE] {mle.mean():.4f}  expected {(n - 1) / n:.4f}"
          f"   E[unbiased] {unbiased.mean():.4f}")
# n=   2  E[MLE] 0.5032  expected 0.5000   E[unbiased] 1.0064
# n=   5  E[MLE] 0.7987  expected 0.8000   E[unbiased] 0.9984
# n=  20  E[MLE] 0.9495  expected 0.9500   E[unbiased] 0.9994
# n= 100  E[MLE] 0.9897  expected 0.9900   E[unbiased] 0.9997

The true variance is 1. Dividing by \(n\) gives 0.5032 at \(n=2\) — half the right answer — while dividing by \(n-1\) gives 1.0064. The shortfall is not random noise; it is exactly the factor \((n-1)/n\), which the third column predicts in advance.

The reason is worth understanding because the same mechanism appears whenever a model is scored on the data it was fitted to. We measured spread around the sample average, and the sample average is by construction the point that sits closest to this particular sample. Distances measured from it are therefore too small. Dividing by \(n-1\) compensates for having spent one measurement’s worth of information locating the center.

Both estimators are consistent. Dividing by \(n\) underestimates the variance by \(1/n\) in expectation: 50% at \(n=2\), 5% at \(n=20\), and 1% at \(n=100\).

Maximum likelihood: which setting best explains what we saw

So far we picked recipes by intuition. There is a general method for inventing them. Suppose the data came from some process with a knob on it, and ask: for which knob setting is the data we actually observed most probable? Turn the knob to that setting and you have the maximum likelihood estimate.

\[\hat{\theta}=\arg\max_{\theta}\sum_i \log p(x_i \mid \theta)\]

\(p(x_i \mid \theta)\) is how plausible one observation is under the knob setting \(\theta\): a probability when the outcome is discrete, a density when it is continuous. Optimizers minimize by convention, so in practice the same quantity is written with a minus sign in front and minimized, which is where the phrase negative log-likelihood comes from.

The sum of logs is there for convenience — probabilities of many independent observations multiply, and logs turn multiplication into addition, which is easier to differentiate and does not underflow to zero on a computer. The \(\arg\max\) is just “the setting that makes this biggest”.

The probability article already noted that common losses are negative log-likelihoods; from the estimation side the same fact reads as a statement about what your loss is estimating. Many common loss functions are a negative log-likelihood wearing a different name. Squared error is what maximum likelihood gives you if you assume the errors are Normal with constant spread. Cross-entropy is what it gives you for a yes/no or pick-one outcome. Poisson deviance is what it gives you for counts.

A likelihood interpretation connects a loss to a noise model, but a loss can also be chosen simply for the errors it penalizes. Squared error penalizes absolute errors, so the same percentage error receives a much larger penalty on a larger forecast: 10% of \$1M against 10% of \$1K. Whether that matches the task depends on whether absolute or relative errors are what matter. Deriving it from a Normal model does not mean the model must be Normal — squared error remains a reasonable way to estimate a conditional mean even when the error variance changes with the inputs.

Maximum likelihood has an asymptotic guarantee attached: under the usual regularity conditions and a correctly specified model, it is asymptotically efficient among regular estimators. These asymptotic guarantees do not imply unbiasedness or minimum MSE at a finite sample size. The variance estimator above shows the first part: dividing by \(n\) is the maximum likelihood answer and it is biased low, badly so at \(n=2\). It is not evidence of worse MSE — the first exercise finds that the MLE has the lower mean squared error of the two.

Adding what you already believed

Maximum likelihood listens only to the data. Often you know something before the data arrives — that a coefficient is probably not enormous, that a rate is probably not 0.99. Maximum a posteriori lets you write that down and adds it to the objective.

\[\hat{\theta}=\arg\max_{\theta}\left[\log p(D \mid \theta)+\log p(\theta)\right]\]

The first term is the likelihood from before. The second, the prior, scores how plausible a parameter value was in advance. And that second term, once you write it out for the usual choices, is regularization — the correspondence is exact, not an analogy.

This gives the tuning knob a modeling interpretation. For a fixed likelihood and scaling convention, a larger \(\lambda\) corresponds to a prior more concentrated near zero. Cross-validation does not express a belief; it selects that strength using held-out performance.

It also marks a limit. MAP returns a single number, and a single number is not a distribution. Whatever uncertainty the Bayesian framing was carrying is discarded the moment you report \(\hat{\theta}\) alone.

How sure are we?

Return to the five height studies. Each produced a single number and none produced 170. A more honest report would be a range wide enough to contain the truth in a stated fraction of repeated studies. That is a confidence interval.

One term first. The standard error of an estimate is the standard deviation of that estimate across repeated samples — exactly the spread column from the first table. For the mean of \(n\) independent, identically distributed measurements with finite variance, it is \(\sigma/\sqrt{n}\), estimated in practice by \(s/\sqrt{n}\); \(\sigma\) is the population standard deviation and \(s\) the sample standard deviation. It describes how much the average moves when you redraw the sample, not how much the individual measurements vary.

The recipe is: take the sample average, and reach out a certain number of standard errors on each side. The code below builds that interval ten thousand times and simply counts how often it caught the truth.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
truth, R = 170.0, 10_000
t_crit = stats.t.ppf(0.975, df=19)
hit_z = hit_t = 0

for _ in range(R):
    s = rng.normal(truth, 8.0, size=20)
    se = s.std(ddof=1) / np.sqrt(20)
    hit_z += abs(s.mean() - truth) <= 1.96 * se
    hit_t += abs(s.mean() - truth) <= t_crit * se

print(f"  multiplier 1.96 (normal)   coverage {hit_z / R:.3f}")
print(f"  multiplier {t_crit:.3f} (t, df=19)  coverage {hit_t / R:.3f}")
# multiplier 1.96 (normal)   coverage 0.936
# multiplier 2.093 (t, df=19)  coverage 0.948

Aiming for 95% and reaching 1.96 standard errors gets you 93.6%. Reaching 2.093 gets you 94.8%. The larger multiplier is correct here because the standard error was itself estimated from the same twenty measurements, and pretending an estimated quantity is known makes the interval too narrow. This is what the \(t\) distribution is for.

Now the part that trips people up, including people who use these intervals daily. The 95% is a property of the procedure, not of the interval in front of you. It says that if you ran the study over and over, 95% of the intervals produced would contain the truth. It does not say that this particular interval has a 95% chance of containing it — the truth is a fixed number and your interval either contains it or does not.

That sounds like hair-splitting until the guarantee fails, which it does quietly. Coverage holds only under the assumptions the recipe was built on. The exact \(t\) interval above relies on independent draws from a Normal distribution. Below, we test two intervals on right-skewed data instead: a normal-approximation interval built with the multiplier 1.96, and a percentile bootstrap interval.

To build the percentile bootstrap interval, draw 30 values with replacement from the 30 observed values, so an observation can appear several times or not at all. Record the resampled mean and repeat this 2,000 times. The 2.5th and 97.5th percentiles of those means become the interval endpoints. In the code, the outer loop generates a new sample from the population to test coverage. Within each iteration, rng.choice(x, size=(2000, 30)) creates the 2,000 bootstrap samples at once; mean(1) computes one mean per resampled row.

import numpy as np

rng = np.random.default_rng(1)
truth = np.exp(0.5)                            # mean of lognormal(0, 1)
normal_hits = bootstrap_hits = 0
R = 2000

for _ in range(R):
    x = rng.lognormal(0, 1, size=30)
    se = x.std(ddof=1) / np.sqrt(30)
    if abs(x.mean() - truth) <= 1.96 * se:
        normal_hits += 1
    resample = rng.choice(x, size=(2000, 30)).mean(1)
    lo, hi = np.percentile(resample, [2.5, 97.5])
    bootstrap_hits += lo <= truth <= hi

print(f"true mean {truth:.4f}")
print(f"normal CI coverage    {normal_hits / R:.3f}")
print(f"bootstrap CI coverage {bootstrap_hits / R:.3f}")
# true mean 1.6487
# normal CI coverage    0.871
# bootstrap CI coverage 0.878

Both intervals claim 95% and deliver about 87%. One in eight is wrong when one in twenty was promised, and nothing in the output warns you. Note that the percentile bootstrap interval is not symmetric, so symmetry alone does not explain the shortfall: at \(n=30\) the distribution of the sample mean of log-normal data is still strongly skewed, and a sample of 30 can fail to capture how much the upper tail contributes to the population mean.

The bootstrap deserves a note, because it is often reached for as the assumption-free option. It makes no assumption about the distribution’s shape, but it does assume the sample is representative enough that resampling it stands in for resampling the population. On heavy-tailed data at moderate \(n\) a sample can fail to represent the tail well, and this percentile interval then inherits the problem it appeared to solve. Other bootstrap variants — studentized or bias-corrected — address part of it.

p-values, and what they do not say

A hypothesis test starts by assuming there is nothing going on — no difference between the two models, no effect of the treatment — and asks how surprising the observed data would be under that assumption. The p-value is the answer: the probability of seeing a result at least this extreme if nothing is going on.

Three things it is not, all of them common in practice: it is not the probability that there is nothing going on, it is not the probability that the result was a fluke, and it says nothing about how large the effect is. A small p-value does not establish practical importance: a p-value of 0.001 on a difference of 0.02 percentage points means that difference was measured precisely, and whether it matters depends on the scale and cost of the application.

The code below runs twenty thousand experiments in which, by construction, nothing is going on — both groups are drawn from the identical distribution.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
p = np.array([stats.ttest_ind(rng.normal(size=30),
                              rng.normal(size=30)).pvalue
              for _ in range(20_000)])

print(f"P(p < 0.05) under the null {np.mean(p < 0.05):.4f}")
print(f"uniformity KS p-value      {stats.kstest(p, 'uniform').pvalue:.3f}")
for m in (1, 10, 20, 100):
    print(f"  {m:3d} tests: P(at least one p < 0.05) = {1 - 0.95 ** m:.4f}")
# P(p < 0.05) under the null 0.0488
# uniformity KS p-value      0.280
#    1 tests: P(at least one p < 0.05) = 0.0500
#   10 tests: P(at least one p < 0.05) = 0.4013
#   20 tests: P(at least one p < 0.05) = 0.6415
#  100 tests: P(at least one p < 0.05) = 0.9941

For an exact continuous test like this one, the p-value under a true null is uniform on \([0,1]\). The second line reports a KS test that does not reject uniformity, which is consistent with it rather than proof of it. Uniformity is why a 0.05 threshold produces a 5% false-positive rate on a single test, and the first line lands on 0.0488.

The last four lines compute \(1-0.95^{m}\), which assumes every null is true and the \(m\) tests are independent. Under those conditions a hundred tests give a 99.4% chance of at least one spurious “significant” result. These assumptions are not guaranteed when comparing model configurations or screening features on the same dataset: the comparisons share rows, and features can be correlated with one another. Dependence among the tests changes the probability of at least one false positive, so the figure shows the direction rather than a number to quote. The screening case is usually not even labeled as testing.

Two standard corrections exist. Bonferroni divides the threshold by the number of tests, which is simple and severe. Benjamini-Hochberg controls the false discovery rate — the expected proportion of false discoveries among the things you declare significant — under independence or certain positive dependence conditions. The two control different quantities, so which one fits depends on whether you need to bound the chance of any false positive at all or the expected share of them.

The third exercise measures what happens when you check a result repeatedly as data arrives, which is one way a well-intentioned analysis produces a false positive.

Exercises

1. A biased estimator that wins. Compare the sample variance divided by \(n-1\), by \(n\), and by \(n+1\) at \(n = 5\), reporting bias, variance, and mean squared error for each. Say which one you would use.

You should get: the unbiased estimator with the worst mean squared error of the three.

import numpy as np

rng = np.random.default_rng(0)
n, R = 5, 200_000
s = rng.normal(0, 1, size=(R, n))
ss = ((s - s.mean(1, keepdims=True)) ** 2).sum(1)

for name, div in (("unbiased 1/(n-1)", n - 1),
                  ("MLE 1/n", n),
                  ("MSE-optimal 1/(n+1)", n + 1)):
    e = ss / div
    print(f"{name:22s} bias {e.mean() - 1:+.4f}"
          f"  var {e.var():.4f}  MSE {((e - 1) ** 2).mean():.4f}")
# unbiased 1/(n-1)       bias +0.0022  var 0.5052  MSE 0.5052
# MLE 1/n                bias -0.1983  var 0.3233  MSE 0.3626
# MSE-optimal 1/(n+1)    bias -0.3319  var 0.2245  MSE 0.3347

2. What the 95% refers to. Simulate 20,000 confidence intervals for a known mean and report the coverage. Then report the interval width at \(n = 25, 100, 400, 1600\) and say what it costs to halve it.

You should get: a coverage close to 0.95 and a width that needs four times the data to halve.

import numpy as np

rng = np.random.default_rng(1)
R, n, sigma = 20_000, 25, 2.0
x = rng.normal(5, sigma, size=(R, n))
half = 1.96 * sigma / np.sqrt(n)

print(f"coverage {np.mean(np.abs(x.mean(1) - 5) <= half):.4f}"
      f"  width {2 * half:.4f}")
for m in (25, 100, 400, 1600):
    print(f"  n={m:5d}  width {2 * 1.96 * sigma / np.sqrt(m):.4f}")
# coverage 0.9486  width 1.5680
#   n=   25  width 1.5680
#   n=  100  width 0.7840
#   n=  400  width 0.3920
#   n= 1600  width 0.1960

3. Optional stopping. Run a two-sample test on pure noise while peeking at the data 1, 2, 5, 10, and 20 times as it accumulates, stopping at the first significant result. Report the false-positive rate for each.

You should get: a nominal 5% test that fires four times as often once you are allowed to look repeatedly.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
R, N = 20_000, 200
A = rng.normal(size=(R, N))                    # no true difference
B = rng.normal(size=(R, N))

for peeks in (1, 2, 5, 10, 20):
    sizes = np.unique(np.linspace(20, N, peeks).astype(int))
    significant = np.zeros(R, bool)
    for k in sizes:
        significant |= stats.ttest_ind(A[:, :k], B[:, :k], axis=1).pvalue < 0.05
    print(f"  {len(sizes):3d} looks: false positive rate {significant.mean():.4f}")
#     1 looks: false positive rate 0.0461
#     2 looks: false positive rate 0.0885
#     5 looks: false positive rate 0.1474
#    10 looks: false positive rate 0.1883
#    20 looks: false positive rate 0.2235

References

• Wasserstein, R. L., & Lazar, N. A. (2016). The ASA Statement on p-Values: Context, Process, and Purpose. The American Statistician. amstat.org
• Benjamini, Y., & Hochberg, Y. (1995). Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society B. Implementation: statsmodels.stats.multitest.fdrcorrection
• Scholz, F. W. Maximum Likelihood Estimation. University of Washington. faculty.washington.edu
• NumPy reference for the variance divisor and the ddof default. numpy.var
• Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall.