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.
An accuracy of 91.2% describes the test cases on which it was measured. How precisely does it estimate performance on new cases from the same population? Another test sample could give a different result. We will study this sampling uncertainty using simple examples before returning to model comparisons. The code uses NumPy and SciPy; the probability article supplies the background on means, variances, and conditional probability.
Parameters, estimators, and estimates
Suppose we want the average height of adults in a country. In this simplified simulation, heights are independent draws from a Normal distribution with mean 170 cm and standard deviation 8 cm. We know the population mean because we set it in the code; a real survey would have to estimate it. We simulate five studies, each measuring twenty people.
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 cmThe five sample means range from 168.53 to 172.82 cm, even though the population mean stays at 170 cm. The differences arise from which observations were sampled.
A sample mean describes the observed people and provides evidence about the population mean under the sampling assumptions. Statistical inference asks how such estimates vary across samples and how far they tend to be from the population quantity. A real survey also needs a sampling design that represents the target population; these simulations assume that problem has been resolved.
The parameter is the population quantity being estimated: the mean of 170 cm here. The estimator is the rule applied to data, such as averaging the measurements. The estimate is its result for a particular sample, such as 168.53 cm. Repeating the sampling changes the estimate while the parameter stays fixed.
Bias and variance across repeated samples
Repeated samples let us distinguish an estimator whose answers are centered away from the truth from one whose answers are widely scattered. Both can contribute to estimation error.
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.372The columns summarize the repeated estimates. 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. The biased ruler has about the same spread as the sample average, but its estimates are centered near 173 cm. In theory their standard deviations are identical: adding a constant changes the mean but not the variance. The small difference between 1.809 and 1.787 comes from simulating separate sets of samples. 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, close to 170 — 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\). Here \(\theta\) is the fixed parameter, \(\hat{\theta}\) is the estimator, and the expectation averages over repeated samples. Variance is how much it scatters around its own average. For a scalar estimator with finite second moment, they combine into the mean squared error, \(\mathbb{E}[(\hat{\theta}-\theta)^2]\). 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}\]
For the height example, the sample average has theoretical variance \(8^2/20=3.2\) cm² and zero bias, so its MSE is 3.2 cm². The biased ruler has the same variance and bias 3 cm, giving MSE \(3.2+3^2=12.2\) cm². A biased estimator can nevertheless have lower MSE than an unbiased alternative if its variance reduction outweighs the added squared bias.
An estimator is consistent if the probability of being farther than any fixed positive tolerance from the truth tends to zero as the sample size grows. This is about increasingly concentrated estimates, not every estimate getting closer at each step. Unbiasedness concerns the expectation at a given sample size. The one-person estimator is unbiased for every \(n\), but inconsistent because it ignores the added observations. Conversely, the variance estimator below that divides by \(n\) is biased at each finite sample size and still consistent.
Why the variance formula divides by n minus one
The sample variance gives a concrete example of bias whose size depends on the sample size. 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.9997The true variance is 1. At \(n=2\), the simulated mean of the estimates is 0.5032 when dividing by \(n\), and 1.0064 when dividing by \(n-1\). Their exact expectations are 0.5 and 1. The columns labeled E[MLE] and E[unbiased] estimate expectations by averaging 200,000 runs; expected gives the theoretical value \((n-1)/n\) for the first estimator.
The sample mean \(\bar X\) minimizes the sum of squared distances to the observed values. Centering on it therefore gives a sum no larger than centering on the unknown population mean \(\mu\). The identity \(\sum_i(X_i-\bar X)^2=\sum_i(X_i-\mu)^2-n(\bar X-\mu)^2\) makes the correction explicit. For \(n>1\) iid observations with finite variance \(\sigma^2\), the expected terms on the right are \(n\sigma^2\) and \(n\operatorname{Var}(\bar X)=\sigma^2\). The expected centered sum is therefore \((n-1)\sigma^2\), which explains the divisor.
Both estimators are consistent. Dividing by \(n\) has bias \(-\sigma^2/n\), or a relative bias of \(-1/n\) when \(\sigma^2>0\). The expected shortfall is 50% at \(n=2\), 5% at \(n=20\), and 1% at \(n=100\).
Maximum likelihood: which setting best explains what we saw
Maximum likelihood starts with a statistical model for the observations. For each allowed parameter setting, it evaluates the probability or density of the data that were observed. A setting that maximizes this likelihood is a maximum likelihood estimate. For independent observations under the model, we can write it as follows:
\[\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.
Independence lets the joint probability or density factor into a product. Taking logs turns that product into a sum without changing which parameter value maximizes it. Computing log probabilities directly also avoids underflow from multiplying many tiny probabilities; computing the product first and then taking its log would not solve that problem. The notation arg max means the parameter setting that maximizes the expression.
For example, let \(\theta\) be the probability of heads on independent coin flips, and observe H, H, H, T. The likelihood is \(\theta^3(1-\theta)\): it equals 0.0625 at \(\theta=0.5\) and about 0.1055 at \(\theta=0.75\). Differentiating its log gives \(3/\theta-1/(1-\theta)\), which is zero at \(\theta=3/4\), the maximum. The likelihood compares parameter settings for fixed data; it is not a probability distribution over those settings.
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. Several common losses are equivalent to negative log-likelihoods up to additive constants and positive scale factors. 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 squared differences on the original scale, 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 require the actual errors to be Normal — squared error remains a reasonable way to target a conditional mean when the relevant second moments are finite, even when the error variance changes with the inputs. A restricted model family may be unable to represent that mean exactly.
Maximum likelihood has an asymptotic guarantee attached: under the usual regularity conditions and a correctly specified model, it is asymptotically efficient among regular estimators: in large samples, its variance reaches the relevant asymptotic lower bound. This is a large-sample result, not a guarantee for every dataset. 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 for Normal data with an unknown mean and it is biased low, badly so at \(n=2\).
Adding what you already believed
Maximum likelihood uses the data within a chosen model. Maximum a posteriori (MAP) estimation also specifies a prior distribution over the parameters, then selects a value where the posterior distribution is largest. For example, a prior can favor moderate coefficients over extremely large ones.
\[\hat{\theta}=\arg\max_{\theta}\left[\log p(D \mid \theta)+\log p(\theta)\right]\]
Here \(D\) denotes the observed dataset. Bayes’ rule gives \(p(\theta\mid D)\propto p(D\mid\theta)p(\theta)\), so maximizing the posterior is equivalent to maximizing the sum of the log-likelihood and log-prior. The normalization constant does not depend on \(\theta\). For the priors in the table, negating the objective turns the coefficient-dependent part of the log-prior into a regularization penalty.
| Prior belief about the coefficients | Log-prior term, up to a constant | Corresponding penalty or objective |
|---|---|---|
| Independent zero-centered Gaussian coefficients, each with variance \(\tau^2\) | \(-\frac{\|\theta\|_2^2}{2\tau^2}\) | L2 penalty (ridge with squared-error loss) |
| Independent zero-centered Laplace coefficients, each with scale \(b\) | \(-\|\theta\|_1/b\) | L1 penalty (lasso with squared-error loss) |
| Constant density over the allowed parameter range | constant on that range | same maximizers as likelihood on that range |
In the table, \(\|\theta\|_2^2=\sum_j\theta_j^2\) and \(\|\theta\|_1=\sum_j|\theta_j|\). The Gaussian prior contributes \(\lambda\|\theta\|_2^2\) when minimizing the negative log-posterior, with \(\lambda=1/(2\tau^2)\) for an unscaled negative log-likelihood. Thus a more concentrated prior gives a stronger penalty. Averaging the data loss instead of summing it changes the matching coefficient. Cross-validation selects a penalty strength using held-out performance. A flat density is still a prior choice: over the entire real line it cannot integrate to 1 and is called an improper prior.
MAP reports a point estimate, possibly a vector of parameters. Reporting that point alone does not describe posterior uncertainty; doing so requires additional information, such as a posterior interval.
How sure are we?
Return to the five height studies. Alongside each point estimate, we can report an interval produced by a procedure designed to contain the population mean in a stated fraction of repeated studies. That is a confidence interval; the coverage claim depends on the procedure’s assumptions.
One term first. The standard error of an estimate is the standard deviation of that estimate across repeated samples — the quantity estimated by the spread column in the repeated-sampling experiment. 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.
For the Normal mean example, the interval is \(\bar x\pm c\,s/\sqrt n\), where \(c\) is the chosen multiplier. If \(\bar x=170\), \(s=8\), and \(n=20\), the estimated standard error is about 1.79 cm. A multiplier of 2.093 gives a margin of about 3.74 cm and an interval of roughly [166.26, 173.74] cm. The following simulation constructs intervals from 10,000 different samples and measures coverage. stats.t.ppf(0.975, df=19) selects the 97.5th percentile of the t distribution, leaving 2.5% in each tail; the degrees of freedom are \(20-1=19\).
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.948In this run, the 1.96 multiplier gives 93.6% coverage and the t multiplier 2.093 gives 94.8%. Under the iid Normal model, using the estimated standard deviation makes the standardized sample mean follow a t distribution with 19 degrees of freedom. Its 97.5th-percentile cutoff, printed as 2.093, gives exactly 95% theoretical coverage; the simulated 94.8% is an estimate of that coverage. Estimating a standard error alone does not make a t interval exact for arbitrary data distributions.
The 95% is a property of the procedure, not of the interval in front of you. Under the stated model, each repetition has a 95% probability of producing an interval that contains the truth; over many independent repetitions, the observed coverage approaches 95%. 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.
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.878Both procedures target 95% coverage, but their estimated coverages here are 87.1% and 87.8%. About one interval in eight misses the mean. The percentile interval can be asymmetric, so the failure cannot be explained just by using a symmetric interval. The log-normal population is strongly right-skewed, and samples of 30 can poorly represent the upper tail that contributes substantially to its mean.
This ordinary bootstrap resamples individual observations as if they were iid draws from the population. It does not specify a parametric distribution shape, but its accuracy still depends on the sampling design, the statistic, and how well the empirical distribution captures relevant population features. Resampling cannot supply large tail values absent from the original sample. Studentized or bias-corrected variants can improve coverage in suitable settings; neither guarantees accurate coverage for every small sample or heavy-tailed population.
p-values, and what they do not say
A hypothesis test specifies a null hypothesis, such as equal population means for two groups, together with a statistical model and a test statistic. The p-value is the probability, under that null model, of a test statistic at least as extreme as the observed one. What counts as extreme depends on the chosen test, including whether it is one-sided or two-sided. A null hypothesis need not mean “no effect”; it can specify a nonzero value too.
A p-value is not the probability that the null hypothesis is true or that the result arose “by chance.” It also does not measure effect size or practical importance. A difference of 0.02 percentage points can have a small p-value in a sufficiently large study, yet still have little operational value. To assess it, report the estimated difference and its uncertainty, then consider the application’s costs. A small p-value alone does not establish that the study measured the intended effect precisely or without bias.
The code below runs twenty thousand experiments in which, by construction, the null hypothesis is true: both groups contain independent observations from the same Normal distribution. The two-sided t-test compares their means.
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.9941For 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. The independent-test calculation therefore does not give their exact false-positive probability. Additional tests can increase the chance of at least one false positive, but the amount depends on their joint behavior; perfectly redundant tests add no new chance.
Two common corrections have different aims. For a fixed family of \(m\) valid tests, Bonferroni uses a threshold of \(\alpha/m\) for each test, controlling the probability of any false rejection at no more than \(\alpha\), without requiring independence. This can be conservative. 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 false-discovery proportion is defined as zero when nothing is declared significant. 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.
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
ddofdefault. numpy.var - Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall.
Discover more from Insightful Data Lab
Subscribe to get the latest posts sent to your email.