Last time, we introduced Stochastic Gradient Descent (SGD), which requires choosing samples at random. Recall from the linear regression note that random number generators in computers depend on the selection of a seed. If we were to start with different seeds, our random selection would be different, which would change how well the parameters get updated, ultimately influencing the final performance of the model.
Running SGD therefore may give different results depending on the selection of seed (we could just have gotten lucky). In this case, how can we evaluate the performance of a model impartially? Am I seeing the results of a lucky draw? Can I expect a similar performance if I run my experiment with a different seed?
To answer these questions, weâd need to run the same experiment multiple times with different seeds and somehow aggregate and summarize the results. In this note, we look at how this is done.
Estimating the Mean via the Law of Large Numbers
Suppose that for $i=0,\dots,n$, we have $X_i$ as a performance metric of our experiment (like the validation accuracy of a logistic regression model). Assume that the variables are identically and independently distributed (i.i.d.), with mean $\mu$ and variance $\sigma^2$.
Then by the Law of Large Numbers, the sample mean $\bar X_n:=\frac{X_1+\dotsb+X_n}{n}$ converges to $\mu$ as $n\to\infty$. This means that the average performance $\bar X_n$ out of $n$ independent experiments becomes a better estimate of the true expected performance $\mu$ as we increase the number $n$ of experiments. This makes intuitive sense: the more experiments we run, the less we rely on luck, and we can more accurately gauge the true performance of the model.
Clearly, thereâs a relationship between $n$ and our approximation accuracy. But how close can these approximations get as a function of $n$?
Getting Confidence Intervals
Expand for a refresher on statistics
The variance of a random variable $X$, written $\mathrm{Var}(X)$ or $\sigma^2$, measures how spread out $X$ tends to be around its mean $\mu$. Itâs defined as the average squared distance from the mean:
$$ \mathrm{Var}(X) := \mathbf E[(X-\mu)^2]. $$
We square the distance rather than just averaging $X-\mu$ directly, because $X-\mu$ is sometimes positive and sometimes negative and so averaging it directly would let those cancel out and wonât yield variance properly. Squaring makes every deviation positive, so they can accumulate.
A reasonable question is: why square, rather than just take the absolute value $\lvert X-\mu \rvert$, which also fixes the sign-cancellation problem and is arguably more intuitive? Both are valid ways to measure spread, but squaring turns out to be much better algebraically. Variance has a clean additivity property: for independent variables, $\mathrm{Var}(X+Y) = \mathrm{Var}(X) + \mathrm{Var}(Y)$, knowing the two variances separately is enough to pin down the variance of the sum exactly. Average absolute deviation doesnât have this property. Moreover, squaring also makes the function $x\mapsto x^2$ smooth (differentiable everywhere, unlike $\lvert x \rvert$ which has a sharp corner at $0$), which, if you recall, matters once youâre optimizing things like loss functions.
The catch is that squaring also changes the units: whatever units $X$ is measured in, $\mathrm{Var}(X)$ ends up in squared units. This is no longer directly comparable to $\mu$, which is still in the original units, so we take the square root:
$$ \sigma := \sqrt{\mathrm{Var}(X)} $$
and this is called the standard deviation. Itâs back in the same units as $X$ and $\mu$, and roughly tells you how far a typical draw of $X$ lands from the mean.
Recall that $\bar X_n$ is our estimate of the true mean $\mu$, and like any estimate from a finite sample, it wonât land exactly on $\mu$. The spread of $\bar X_n$ just means: across the many different samples we could have drawn, how much does $\bar X_n$ tend to bounce around? A small spread means $\bar X_n$ is a more reliable estimate of $\mu$; a large spread means we got somewhat lucky or unlucky with our particular sample, and $\bar X_n$ could easily have come out quite different. This is the quantity weâll need to build a confidence interval.
Because variance is additive over independent variables, i.e. $\mathrm{Var}(X+Y)=\mathrm{Var}(X)+\mathrm{Var}(Y)$ when $X,Y$ are independent, and that constants pull out squared, i.e. $\mathrm{Var}(cX) = c^2\mathrm{Var}(X)$, we can compute the spread of $\bar X_n$ directly:
$$ \begin{aligned} \mathrm{Var}(\bar X_n) &= \mathrm{Var}\left(\frac{X_1+\dotsb+X_n}{n}\right) \\ &= \frac{1}{n^2}\mathrm{Var}(X_1+\dotsb+X_n) \\ &= \frac{1}{n^2}\Big(\mathrm{Var}(X_1)+\dotsb+\mathrm{Var}(X_n)\Big) \\ &= \frac{1}{n^2}\sum_{i=1}^n \mathrm{Var}(X_i) \\ &= \frac{1}{n^2}\cdot n\sigma^2 \\ &= \frac{\sigma^2}{n}. \end{aligned} $$
Taking the square root of both sides gives the standard deviation of $\bar X_n$:
$$ \sqrt{\mathrm{Var}(\bar X_n)} = \frac{\sigma}{\sqrt n}. $$
This is usually called the standard error to distinguish it from $\sigma$ itself. Itâs the spread of our estimate $\bar X_n$.
Z-scores: suppose we want to know whether $\bar X_n$ landed close to $\mu$ or far from it. âCloseâ and âfarâ only mean something relative to the standard error (being $1$ unit off doesnât sound like much if the standard error is $10$, but sounds huge if the standard error is $0.01$). So instead of looking at the raw difference $\bar X_n - \mu$, we divide it by the standard error:
$$ \frac{\bar X_n - \mu}{\sigma/\sqrt n}. $$
This is called a Z-score. Itâs useful because itâs unit-free: a Z-score of $2$ means âtwo standard errors away from expected,â which is unusual regardless of whether weâre measuring accuracy, seconds, or dollars. This is what lets us attach a probability to how far off an estimate is.
Motivating the $t$-Distribution
If we happened to know the true standard deviation $\sigma$, life would be easy. Assuming (for now) that each $X_i \sim \mathscr N(\mu, \sigma)$, the quantity
$$ Z_n:=\frac{\bar X_n-\mu}{\sigma/\sqrt n} $$
would follow a standard normal distribution exactly, and we could read off confidence intervals directly from a normal table.
The issue is that we never actually know $\sigma$. All we have is the data itself, so we estimate $\sigma^2$ with the sample variance
$$ S_n^2=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar X_n)^2. $$
We actually canât just plug in $S_n$ for $\sigma$. $S_n$ is itself a noisy estimate computed from the same finite sample, so weâve now stacked two sources of randomness on top of each other (the âwobbleâ of $\bar X_n$ around $\mu$, and the âwobbleâ of $S_n$ around $\sigma$). The resulting statistic
$$ T_n:=\frac{\bar X_n-\mu}{\sqrt{S^2_n / n}} $$
is therefore not exactly normal; it has heavier tails, since every so often $S_n$ underestimates $\sigma$ and inflates $T_n$.
It turns out $T_n$ follows a Student $t$-distribution with $n-1$ degrees of freedom1.
As $n\to\infty$, $S_n$ becomes an increasingly accurate estimate of $\sigma$, so the extra noise from estimating $\sigma$ washes out and $T_n$âs distribution converges back to $\mathscr N(0,1)$.
Constructing the Interval
We get a relation between the precision of our estimate and our confidence with that precision:
$$ \mathbf P(\lvert\bar X_n-\mu\rvert<c\sqrt{S^2_n / n}) =\mathbf P(\lvert T_n\rvert<c)=\alpha. $$
Expand to see derivation
Recall from above that $T_n := \dfrac{\bar X_n - \mu}{\sqrt{S_n^2/n}}$. We can rearrange this to isolate $\bar X_n - \mu$: multiplying both sides by $\sqrt{S_n^2/n}$ gives $\bar X_n - \mu = T_n\sqrt{S_n^2/n}$, and taking absolute values,
$$ \lvert\bar X_n - \mu\rvert = \lvert T_n\rvert\sqrt{S_n^2/n}. $$
So the event $\lvert\bar X_n-\mu\rvert<c\sqrt{S_n^2/n}$ (our estimate landing within $c$ standard errors of $\mu$) is exactly the same event as $\lvert T_n\rvert<c$. So their probabilities are equal too:
$$ \mathbf P(\lvert\bar X_n-\mu\rvert<c\sqrt{S^2_n / n}) =\mathbf P(\lvert T_n\rvert<c). $$
The right-hand side is useful because we know $T_n$âs distribution exactly (Student $t$ with $n-1$ degrees of freedom, from above), so we can choose $c$ to make this probability equal to whatever confidence level $\alpha$ we want:
$$ \mathbf P(\lvert T_n\rvert<c)=\alpha. $$
Since the middle expression is a value of the CDF (cumulative distribution function, the function giving $\mathbf P(T_n < x)$ for each $x$) of a known distribution, $c$ can be looked up (or computed numerically) once we fix a $\alpha$ (confidence level) and $n$ (sample size)2.
In words: if we repeated the entire experiment (running $n$ trials, computing $\bar X_n$ and $S_n$) many times over, then in $\alpha$ proportion of those repetitions weâd have $\lvert\bar X_n-\mu\rvert<cS_n/\sqrt n$. The interval
$$ (\bar X_n-cS_n/\sqrt n,\ \bar X_n+cS_n/\sqrt n) $$
is called a confidence interval, and itâs customary to report the result as $\bar X_n\pm cS_n/\sqrt n$ (stating the confidence level $\alpha$ somewhere nearby).
Dropping the Normality Assumption
Everything above assumed each $X_i$ was individually normal. The LindenbergâLĂ©vy Central Limit Theorem lets us relax this: regardless of the distribution of the $X_i$ (as long as it has finite variance), the sequence $\sqrt n(\bar X_n-\mu)$ converges in distribution to $\mathscr N(0,\sigma^2)$ as $n\to\infty$. That is, for every $x\in\mathbf R$,
$$ \mathbf P(\sqrt n(\bar X_n-\mu) < x) \xrightarrow{n\to\infty} \mathbf P(\mathscr N(0, \sigma^2) < x). $$
So even without assuming normality up front, $\bar X_n$ behaves approximately normally for large $n$.
If you want a quick mental shortcut for confidence intervals under normality, the 68-95-99.7 rule is handy:
$$ \mathbf P(|\mathscr N(\mu, \sigma^2) - \mu| < \sigma) \approx 68.27\% $$
$$ \mathbf P(|\mathscr N(\mu, \sigma^2) - \mu| < 2\sigma) \approx 95.45\% $$
$$ \mathbf P(|\mathscr N(\mu, \sigma^2) - \mu| < 3\sigma) \approx 99.73\% $$
Writing an Ensemble Training Loop
Notice that throughout this note, weâve been assuming access to $n$ independently trained models (one per random seed). This collection is called an ensemble. Simply put, this is a set of models trained on the same task, differing only in some source of randomness (like the seed, or the data order). Ensembles show up for a few reasons:
- Performance metrics with confidence intervals, you need multiple independent trials to estimate spread at all, and an ensemble is that set of trials (what we just saw above).
- Aggregating predictions from an ensemble (e.g. averaging their outputs, or taking a majority vote) usually beats any individual model in the ensemble. The intuition is similar to why $\bar X_n$ is a better estimate of $\mu$ than any single $X_i$: each modelâs predictions carry some noise or bias tied to its particular random seed, and averaging across independently-seeded models tends to cancel out the noise thatâs uncorrelated across them, leaving something closer to the âtrueâ underlying signal.
- Comparing hyperparameters, training several models that differ in some hyperparameter (like learning rate) lets us see which setting performs best, using the same multiple-trials machinery from this note to know whether an apparent difference is real or just noise.
When we write the ensemble training loop, weâll focus on making it efficient: instead of training the ensemble models one by one in a Python loop, weâll vectorize the process. Weâll be training all $n$ models simultaneously by treating the ensemble as an extra batch dimension. A naive loop wastes the hardwareâs ability to do many identical operations at once: GPUs are built around massive parallelism, and even on CPU, SIMD vectorization lets a single instruction operate on multiple pieces of data simultaneously. Training $n$ models sequentially ignores all of this, so we get a substantial speedup by restructuring the computation to update every model in the ensemble in one shot.
In Practice
Suppose we train a small logistic regression model with $n=10$ different random seeds, and record the validation accuracy each time:
import numpy as np
from scipy import stats
# validation accuracies from 10 runs with different seeds (these are made up btw)
accuracies = np.array([0.842, 0.851, 0.839, 0.847, 0.855,
0.836, 0.849, 0.844, 0.858, 0.841])
n = len(accuracies)
x_bar = accuracies.mean() # sample mean, \bar X_n
s = accuracies.std(ddof=1) # sample std, S_n (ddof=1 for the n-1 in the denominator)
se = s / np.sqrt(n) # standard error, S_n / sqrt(n)
alpha = 0.95
# We want c such that P(|T_n| < c) = alpha, where T_n ~ Student-t with n-1 degrees of freedom.
# Since this is two-sided, the leftover probability (1 - alpha) is split evenly
# between the left and right tails, so we need the (1 - (1-alpha)/2) = (1+alpha)/2 quantile.
# e.g. alpha = 0.95 -> we want the 0.975 quantile, leaving 2.5% in each tail.
quantile_needed = (1 + alpha) / 2
c = stats.t.ppf(quantile_needed, df=n - 1) # inverse CDF: "value below which `quantile_needed` of the mass lies"
margin = c * se
print(f"xÌ = {x_bar:.4f}, S = {s:.4f}, SE = {se:.4f}")
print(f"quantile needed = {quantile_needed}, c = {c:.3f}")
print(f"95% CI: ({x_bar - margin:.4f}, {x_bar + margin:.4f})")
This gives:
| Quantity | Value |
|---|---|
| Sample mean, $\bar{X}_n$ | $0.8462$ |
| Sample standard deviation, $S_n$ | $0.0071$ |
| Standard error, $S_n/\sqrt{n}$ | $0.0022$ |
| Critical value ($\alpha=95%$, $df=9$), $c$ | $2.262$ |
| Margin of error, $cS_n/\sqrt{n}$ | $0.0050$ |
| 95% Confidence interval | $(0.8412,\ 0.8512)$ |
So weâd report the modelâs accuracy as $0.846 \pm 0.005$ (95% CI), meaning if we reran this whole process (10 fresh seeds, recompute $\bar X_n$ and $S_n$) many times, about 95% of those intervals would contain the true expected accuracy $\mu$.
A note on one-sided vs. two-sided:
The code above uses quantile_needed = (1 + alpha) / 2, not regular alpha. This is because the interval we built is two-sided. We asked for the probability that $\bar X_n$ falls within $c$ standard errors of $\mu$ in either direction, i.e. $\mathbf P(\lvert T_n \rvert < c) = \alpha$, which splits the leftover $1-\alpha$ probability evenly into both tails ($\frac{1-\alpha}{2}$ on each side). This is the right choice when you have no reason to expect the estimate to be biased in a particular direction; you just want to know how far off it could be either way, which is the usual case when reporting a validation accuracy.
A one-sided interval would instead ask a directional question, like âwhatâs a lower bound $\mu$ is at least $\alpha$ proportion of the time?â. This is useful if you only care about one direction (e.g. a safety-critical setting where you only care whether accuracy might be worse than expected, not better). This version puts the entire $1-\alpha$ probability in a single tail, so quantile_needed would just be alpha itself (stats.t.ppf(alpha, df=n-1)), and the resulting $c$ will be smaller than the two-sided one at the same $\alpha$.
stats.t.ppf(q, df) is the percent point function, which is the inverse of the CDF. Where the CDF answers âwhatâs the probability of landing below $x$?â, ppf answers the reverse question: âwhat value $x$ has probability $q$ of landing below it?â Thatâs exactly the $c$ we need from the Constructing the Interval section above, computed numerically instead of read off a printed $t$-table.
-
The $n-1$, rather than $n$, comes from the fact that the $n$ residuals $X_i-\bar X_n$ used to build $S_n^2$ arenât fully independent (theyâre constrained to sum to zero, since $\bar X_n$ was computed from the same data). So only $n-1$ of them carry independent information. Learn more about degrees of freedom here. ↩
-
For example, at $\alpha=95\%$ with $n=10$, we get $c\approx 2.262$. ↩