Statistical Inference, from Scratch | Oh My Kode

Statistical Inference, from Scratch

18 Feb 2024

68 minutes read

You are cooking a large pot of soup. You want to know whether it needs salt. You do not drink the pot — you stir it, take one spoonful, and taste. That spoonful is a sample, the pot is the population, and the leap you just made without thinking about it is the whole of statistics : using a small piece to say something honest about a whole you will never see. This post builds that leap carefully, in three passes. The first needs nothing but arithmetic. The second adds the two ideas everybody misuses — confidence intervals and p-values — and takes the time to say what they actually mean. The third goes underneath, to where the formulas come from and where they stop working. Read the first part and you will understand what a statistician is doing ; read all three and you will be able to check their work.

PART I — THE FOUNDATIONS
nothing here needs more than arithmetic

1. The one problem statistics solves

Everything that follows is an answer to a single question :

There is a number I care about. I cannot measure it, because measuring it would mean examining everyone or everything. I can only afford to look at a few. What can I honestly say ?

Two words carry all the weight, so let us fix them properly.

The population is the entire collection you actually care about — every customer, every patient, every possible roll of the die. A parameter is a number describing it : the true average $\mu$, the true proportion $p$. Parameters are fixed and unknown.

A sample is the handful you actually observed, $x_1, \dots, x_n$. A statistic is any number you compute from it, such as the sample average $\bar{x}$. Statistics are known and random — random because a different sample would have given a different one.

In plain words : the pot has one true saltiness, and it is not going to change while you think about it. Your spoonful, on the other hand, could have come out saltier or blander by luck. The truth is fixed ; your view of it wobbles.

That single sentence is the hinge of the entire subject, and almost every classic misunderstanding comes from getting it backwards. Say it once more, because we will need it in section 4 : the parameter does not move. The data moves.

the population — you never get to see this four samples of 12 — this is all you ever get sample 1 sample 2 sample 3 sample 4 μ — fixed, unknown, and it never moves
Figure 1 - The whole setup in one picture. The dashed line is the truth : one number, fixed, invisible. Each sample is a handful of draws, and each produces its own estimate (the orange triangle). The estimates scatter around the line — that scatter is the only thing statistics is ever really talking about.

One more thing the soup teaches, and it is not a detail. You stirred before tasting. If you had skimmed a spoonful off the top of an unstirred pot you would have learned about the top of the pot, not the pot. Stirring is what makes the spoonful representative, and its statistical name is random sampling : every member of the population must have a known, non-zero chance of ending up in your sample. Everything in this post assumes you stirred. Section 12 is about what happens when you did not.

2. Estimation : making the guess

The first job is to produce a number. This is estimation, and it comes with a distinction worth keeping straight :

  • an estimator is the recipe — “add them up and divide by how many” ;
  • an estimate is the number that recipe produced this time — $\bar{x} = 4.2$.

The recipe is fixed ; the number it produces is random, because the data is. So we judge the recipe, not the number. And there are exactly two questions to ask of a recipe.

Question one : does it aim at the right place ? Run the recipe on every possible sample and average the answers. If that average lands on the truth, the recipe is unbiased. For the sample mean this is a one-line calculation :

\[\begin{equation} \mathbb{E}[\bar{X}] = \mathbb{E}\!\left[\frac{1}{n}\sum_{i=1}^{n} X_i\right] = \frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[X_i] = \frac{1}{n}\cdot n\mu = \mu. \label{eq:unbiased} \end{equation}\]

In plain words : the spoonful is sometimes saltier than the pot and sometimes blander, but there is no systematic pull in either direction. Over many spoonfuls the errors cancel.

Question two : how much does it wobble ? Unbiased is not enough — a recipe that returns “the truth plus or minus a thousand” is unbiased and useless. We want the spread of the estimates around the truth. Using the fact that variances of independent things add :

\[\begin{equation} \operatorname{Var}(\bar{X}) = \operatorname{Var}\!\left(\frac{1}{n}\sum_{i=1}^{n} X_i\right) = \frac{1}{n^{2}}\sum_{i=1}^{n}\operatorname{Var}(X_i) = \frac{n\sigma^{2}}{n^{2}} = \frac{\sigma^{2}}{n}. \label{eq:var} \end{equation}\]

Take the square root to get back to the original units and you have the single most useful quantity in applied statistics :

The standard error of the sample mean is the typical distance between your estimate and the truth : $$ \operatorname{SE}(\bar{X}) \;=\; \frac{\sigma}{\sqrt{n}}. $$ It is not the spread of the data ($\sigma$ is that). It is the spread of the estimate.

Read $\eqref{eq:var}$ again, because it contains a law that governs the cost of every study, survey and A/B test ever run. The error shrinks like $\sqrt{n}$, not like $n$.

To halve the standard error you must quadruple the sample. To divide it by ten you need a hundred times the data.
0.50 n = 4 0.25 n = 16 0.125 n = 64 each marked step multiplies n by 4 — and only halves the error sample size n standard error, in units of σ
Figure 2 - The square-root law, and why data gets expensive. Going from 4 observations to 16 halves the standard error ; halving it again costs another 48. The curve flattens fast — which is exactly why a poll of 1000 people and a poll of 1500 give you almost the same precision.

In plain words : the first hundred observations teach you a great deal ; the next hundred teach you much less ; going from ten thousand to eleven thousand is nearly pointless. This is why polls of a thousand people report a margin of about three points, and why polling ten thousand — ten times the cost — only takes it to one point.

You have met this law before if you read the HyperLogLog post : that algorithm's relative error is $1.04/\sqrt{m}$ with $m$ registers. Same $\sqrt{\cdot}$, same bargain. It is not a coincidence — HyperLogLog is estimating a parameter from a random sketch, and it pays the same price everyone pays.

Bias and wobble combine into one honest measure of a recipe’s quality, the mean squared error :

\[\begin{equation} \operatorname{MSE}(\hat\theta) = \mathbb{E}\big[(\hat\theta - \theta)^{2}\big] = \underbrace{\big(\mathbb{E}[\hat\theta] - \theta\big)^{2}}_{\text{bias}^{2}} + \underbrace{\operatorname{Var}(\hat\theta)}_{\text{wobble}}. \label{eq:mse} \end{equation}\]

Keep this decomposition in mind : collecting more data shrinks the second term and does absolutely nothing to the first. We will come back to that with some force in section 12.

3. Why the average behaves so well

Two theorems explain why the sample mean is the workhorse of the whole subject. Neither needs a proof here, but both deserve to be understood rather than memorised.

The law of large numbers says the wobble eventually goes away : as $n$ grows, $\bar{X}_n$ settles onto $\mu$ and stays there. That is the formal version of “a bigger spoonful is more reliable”, and $\eqref{eq:var}$ already told us why — the variance $\sigma^2/n$ goes to zero.

The second is genuinely surprising, and it is the reason the same handful of formulas works everywhere.

Let $X_1, \dots, X_n$ be independent draws from any distribution with mean $\mu$ and finite variance $\sigma^{2}$. Then as $n$ grows, $$ \frac{\bar{X}_n - \mu}{\sigma/\sqrt{n}} \;\xrightarrow{\;d\;}\; \mathcal{N}(0, 1), $$ that is, the sample mean behaves like a bell curve centred at $\mu$ with standard deviation $\sigma/\sqrt{n}$ — whatever the shape of the original population.

In plain words : the population can be wildly lopsided, spiky, or shaped like a camel. Average enough draws from it and the averages you get will pile up into a symmetric bell. The lopsidedness averages away.

n = 1 the raw population n = 5 already less lopsided n = 30 a bell, near enough the vertical line is the true mean ; it is in the same place every time
Figure 3 - The central limit theorem at work. The population on the left is as lopsided as it gets — an exponential, all its mass bunched near zero with a long tail. Average just 5 draws and the lopsidedness is already fading ; average 30 and the distribution of the average is almost exactly the bell curve drawn in dashes. Each panel is scaled to its own height.

This is what makes the subject practical. You almost never know the shape of the population — but you do not need to. As long as you are working with an average of enough independent observations, you know the shape of its wobble. Every formula in Part II is built on that one gift.

"Enough" is not a fixed number. For a roughly symmetric population, $n \approx 15$ is plenty. For something as skewed as an exponential, $n \approx 30$ is the traditional rule of thumb. For very heavy tails — income, insurance claims, network traffic — even $n = 1000$ can be too few, and section 12 explains why.
PART II — SAYING HOW SURE YOU ARE
the two tools everybody uses and most people misquote

4. Confidence intervals : a range instead of a number

An estimate alone is a bad answer, because it hides how much it wobbles. “The average is 4.2” and “the average is 4.2, give or take 0.05” and “the average is 4.2, give or take 3” are three completely different states of knowledge. A confidence interval, an idea Neyman set out precisely in 1937 (Neyman, 1937), reports the give-or-take.

Building one is now just algebra on the central limit theorem. The CLT says $Z = (\bar{X} - \mu)/(\sigma/\sqrt{n})$ is a standard bell, and a standard bell sits between $-1.96$ and $+1.96$ about 95 % of the time :

\[\mathbb{P}\!\left(-1.96 \;\le\; \frac{\bar{X} - \mu}{\sigma/\sqrt{n}} \;\le\; 1.96 \right) = 0.95 .\]

Multiply through by $\sigma/\sqrt{n}$, subtract $\bar{X}$, multiply by $-1$ (which flips the inequalities), and the same statement reads :

\[\begin{equation} \mathbb{P}\Big(\bar{X} - 1.96\,\tfrac{\sigma}{\sqrt{n}} \;\le\; \mu \;\le\; \bar{X} + 1.96\,\tfrac{\sigma}{\sqrt{n}}\Big) = 0.95 . \label{eq:ci} \end{equation}\]

So the recipe is estimate $\pm$ 1.96 standard errors. Nothing more.

Now the part that matters more than the formula. Look hard at $\eqref{eq:ci}$ and ask what is random in it. Not $\mu$ — $\mu$ is a fixed number that has never moved. The random things are $\bar{X}$, and therefore the two endpoints. The interval is what jumps around ; the truth stands still.

The number 95 % is a property of the procedure, not of the interval in front of you. It says : if you repeated the whole study many times, about 95 % of the intervals you built this way would contain $\mu$. Your particular interval either contains $\mu$ or it does not — there is no probability left in it once the data is in.
μ — the truth, fixed this one misses 19 of these 20 intervals cover μ — but from the inside, none of them can tell
Figure 4 - Twenty studies of the same fixed quantity. Every interval was built by the identical, correct procedure — yet one of them does not contain μ. That is not a mistake, it is the 5 % working as advertised. The confidence is in the method, not in any one interval, and nothing about the unlucky one looks wrong from the inside.

Figure 4 is worth more than any sentence about this. Twenty studies, twenty intervals, one fixed truth. Most of the intervals cover it ; one does not. Nothing about the unlucky interval looks different from the inside — that is the entire point, and it is why “95 % confidence” cannot be a statement about the interval you happen to be holding.

The width of the interval is $2 \times 1.96 \times \sigma/\sqrt{n}$, so the square-root law from section 2 governs precision too : to make your interval half as wide, get four times the data.

5. When you do not know $\sigma$ : Student’s t

Formula $\eqref{eq:ci}$ has a cheat in it. It uses $\sigma$, the true spread of the population — which you do not know, for the same reason you do not know $\mu$. In practice you estimate it from the sample :

\[\begin{equation} s^{2} = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^{2}. \label{eq:s2} \end{equation}\]

Why $n-1$ and not $n$ ? Because the deviations are measured from $\bar{x}$, not from the true $\mu$, and $\bar{x}$ sits by construction as close to your data as possible. That makes the squared deviations a little too small. Dividing by $n-1$ corrects exactly for it. The intuitive bookkeeping : you had $n$ independent numbers, you spent one of them working out where the centre was, so $n-1$ pieces of information remain — the degrees of freedom.

Substituting an estimate for $\sigma$ adds a second source of wobble, so the multiplier $1.96$ is no longer right — it is too optimistic. Gosset worked out the correct distribution in 1908, publishing under the pen name Student because his employer, Guinness, would not let him publish under his own (Student, 1908). The t-distribution with $n-1$ degrees of freedom looks like a bell with fatter tails, and the fatter the fewer the observations :

sample size $n$ multiplier for 95 %
5 2.78
10 2.26
30 2.05
100 1.98
$\infty$ 1.96

In plain words : with five observations you must be considerably more cautious than the textbook $1.96$, because you are unsure about the spread as well as the centre. By $n = 30$ the correction is small, and by $n = 100$ it hardly matters. Use $t$ by default ; it is right for small samples and costs nothing for large ones.

import math, statistics

def mean_ci(xs, mult=1.96):
    """95 % interval for a mean. Use a t table for small n."""
    n = len(xs)
    m = statistics.fmean(xs)
    se = statistics.stdev(xs) / math.sqrt(n)
    return m, (m - mult * se, m + mult * se)

6. Hypothesis tests : the same picture, read backwards

A confidence interval answers “what values are plausible ?”. A hypothesis test answers a narrower question : “is this specific value still plausible ?” — usually the boring value, the one meaning nothing is happening.

The logic is proof by contradiction, softened with a probability :

  1. Assume nothing is happening. This is the null hypothesis $H_0$ — the coin is fair, the drug does nothing, the two versions convert equally.
  2. Ask how surprising your data would be if that were true.
  3. If it would be very surprising, stop believing $H_0$. If it would not, you have learned nothing and you say so.

In plain words : a courtroom. The defendant is presumed innocent ; the evidence is weighed against that presumption ; only if the evidence would be very unlikely for an innocent person do we convict. And crucially — an acquittal is not a proof of innocence, it is a failure to prove guilt.

The measure of surprise is the p-value, and it deserves a formal definition because nearly everything said about it in casual conversation is wrong.

Let $T$ be a test statistic — a number summarising how far the data sits from what $H_0$ predicts. The p-value is $$ p \;=\; \mathbb{P}\big(\,T \text{ at least as extreme as observed} \;\mid\; H_0 \text{ is true}\,\big). $$ It is the probability of the data given the hypothesis. It is not the probability of the hypothesis given the data.

That last sentence is the one to tattoo somewhere. Conditional probabilities do not flip. Almost every dog has four legs, so $\mathbb{P}(\text{four legs} \mid \text{dog}) \approx 1$ ; but plenty of four-legged things are not dogs, so $\mathbb{P}(\text{dog} \mid \text{four legs})$ is small. A p-value is the first kind of statement. Reading it as the second — “p = 0.03, so there is a 3 % chance the drug does nothing” — is the single most common error in the applied literature (Wasserstein & Lazar, 2016).

if nothing were going on, the statistic would land like this observed t = 2.3 p = 0.021 — the shaded area in both tails observed t = 1.6 p = 0.110 — the shaded area in both tails observed t = 1.0 p = 0.317 — the shaded area in both tails the further out the observation, the smaller the shaded area, the more surprising the data
Figure 5 - What a p-value is. The curve is where the test statistic would fall if the null hypothesis were true. The shaded tails are everything at least as extreme as what you actually observed, and their area is $p$. Note what the picture does not contain : any information about whether $H_0$ is true. It only says how unusual your data would be if it were.

The mechanics are the same algebra as before. For a mean :

\[\begin{equation} t = \frac{\bar{x} - \mu_0}{s/\sqrt{n}} \quad\text{— how many standard errors is the data from the null value ?} \label{eq:tstat} \end{equation}\]

Then $p$ is the area in the tails of the null distribution beyond $\pm t$. Compare $p$ to a threshold $\alpha$ (conventionally $0.05$, a number Fisher chose for convenience (Fisher, 1925) and which has no deeper justification 1) and reject or fail to reject. Never accept — see the courtroom.

Now the observation that ties Part II together, and which most courses leave as a footnote.

A 95 % confidence interval is exactly the set of null values $\mu_0$ that a two-sided test at $\alpha = 0.05$ would fail to reject.

Proof.   The test rejects $\mu_0$ precisely when the statistic exceeds the critical value : $$ \left\lvert \frac{\bar{x} - \mu_0}{\operatorname{SE}} \right\rvert > 1.96 \iff \lvert \bar{x} - \mu_0 \rvert > 1.96\operatorname{SE} \iff \mu_0 \notin \big[\bar{x} - 1.96\operatorname{SE},\; \bar{x} + 1.96\operatorname{SE}\big].$$ So $\mu_0$ is not rejected exactly when it lies inside the interval $\eqref{eq:ci}$.

$\square$

In plain words : the interval already contains every test you could have run. Look at whether it covers the boring value and you have your answer — with the size of the effect and its uncertainty thrown in for free. This is why experienced people report intervals rather than p-values : the interval is strictly more informative, and it is much harder to misread.

7. The two ways to be wrong

A test makes a decision, and decisions can go wrong in two directions. Naming both, and reframing testing as a choice between two kinds of mistake rather than a verdict on truth, is Neyman and Pearson’s contribution (Neyman & Pearson, 1933). The whole design of an experiment is about choosing how much of each error to tolerate.

  $H_0$ is true (nothing happening) $H_0$ is false (something real)
reject $H_0$ Type I error — false alarm, rate $\alpha$ correct detection, rate $1-\beta$
fail to reject correct silence Type II error — missed it, rate $\beta$

In plain words : a Type I error is crying wolf. A Type II error is missing the wolf. A smoke alarm set to be very sensitive rarely misses a fire but goes off when you make toast ; set it less sensitive and the toast is fine but you might sleep through something real.

The power of a test is $1 - \beta$, the probability of detecting an effect that is genuinely there. Power depends on four things and nothing else : the size of the real effect $\Delta$, the noise $\sigma$, the sample size $n$, and the threshold $\alpha$.
decision threshold if nothing is happening if the effect is real β — the wolf you miss α — crying wolf here β = 0.15, so the power is 85% — the two shaded areas only shrink together if n grows
Figure 6 - The two errors, and why they fight. Everything right of the threshold is called a detection : that is the orange sliver of the null curve (false alarms, rate α) plus most of the green curve (correct detections). Everything left of it is called silence, including the blue region under the green curve — real effects you missed, rate β. Slide the threshold either way and one area grows as the other shrinks. Only more data, which narrows both curves, improves both at once.

Figure 6 shows why $\alpha$ and $\beta$ trade off. Move the decision threshold left and you catch more real effects but raise more false alarms ; move it right and the reverse. Only one thing improves both at once, and it is more data — more data narrows both curves, so they overlap less.

That relationship can be solved for $n$. To detect a difference $\Delta$ between two groups with power $1-\beta$ at level $\alpha$, you need about

\[\begin{equation} n \;\approx\; \frac{2\,\sigma^{2}\,\big(z_{1-\alpha/2} + z_{1-\beta}\big)^{2}}{\Delta^{2}} \quad\text{per group.} \label{eq:power} \end{equation}\]
import math

def n_per_group(delta, sd, z_alpha=1.96, z_beta=0.84):
    """z values shown are alpha = 0.05, power = 0.80."""
    z = z_alpha + z_beta
    return math.ceil(2 * (z * sd / delta) ** 2)

n_per_group(delta=0.5, sd=1.0)     # 63 per group
n_per_group(delta=0.1, sd=1.0)     # 1570 per group

Look at the $\Delta^{2}$ in the denominator : halving the effect you want to catch multiplies the required sample by four. Small effects are expensive, and this formula should be run before the experiment, not after.

Which leads to the most useful thing in this section :

A non-significant result from an underpowered study means nothing at all. If your design had a 30 % chance of detecting the effect even when it was real, failing to detect it is not evidence of absence — it is the expected outcome either way. "We found no significant difference" is only informative alongside "and we had 90 % power to find a difference of this size".
PART III — THE FINE PRINT
where the formulas come from, and where they break

8. Significance is not importance

A p-value answers whether there is an effect, never how big. And because the standard error shrinks like $1/\sqrt{n}$, a large enough study makes any non-zero difference significant. Test a landing page on ten million users and a conversion improvement of $0.001$ percentage points will come back with $p < 0.001$ — real, reproducible, and completely worthless.

So always report the effect size alongside. In raw units when they mean something (“4.10 euros more per order”), or standardised when they do not — Cohen’s $d = (\bar{x}1 - \bar{x}_2)/s{\text{pooled}}$, conventionally small at $0.2$, medium at $0.5$, large at $0.8$ (Cohen, 1988). Better still, report the confidence interval, which carries the size and the uncertainty in one object.

There is a matching trap in the other direction. A wide interval that includes zero does not mean “no effect” ; it means “we do not know”. Distinguishing “the effect is somewhere between $-0.1$ and $+0.1$, so it is negligible” from “the effect is somewhere between $-8$ and $+9$, so we learned nothing” is impossible from a p-value and immediate from an interval.

9. The multiple comparisons trap

Here is the mistake that quietly invalidates more analyses than any other. A test at $\alpha = 0.05$ raises a false alarm 5 % of the time per test. Run many tests and the alarms accumulate. If the tests are independent and nothing at all is going on,

\[\begin{equation} \mathbb{P}(\text{at least one false alarm}) = 1 - (1 - \alpha)^{m}. \label{eq:fwer} \end{equation}\]
for m in (1, 5, 20, 100):
    print(m, round(1 - 0.95 ** m, 3))
# 1 0.05   5 0.226   20 0.642   100 0.994

Twenty comparisons and it is more likely than not that you will find something. Segment an A/B test by country, device, browser and hour and you have run dozens of tests without ever intending to — and the one “significant” cell you then report is the expected output of pure noise.

The blunt fix is Bonferroni : test each hypothesis at $\alpha/m$ instead of $\alpha$, which keeps the overall false-alarm rate under $\alpha$. It is simple and provably safe, but conservative — with many tests it destroys your power. When you have hundreds or thousands of comparisons, the modern choice is to control the false discovery rate instead : rather than demanding no false alarms, accept that a controlled fraction (say 5 %) of your reported discoveries will be wrong. The Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995) does this and is far more powerful.

The subtler version of this trap needs no explicit tests at all. If you try a few outlier rules, a couple of transformations, and three ways of grouping, then report the analysis that “worked”, you have searched a large space of tests while performing only one — Gelman and Loken’s garden of forking paths (Gelman & Loken, 2014). The only real defence is to decide the analysis before seeing the outcome, and to label anything found afterwards as exploratory. This mechanism is a large part of why so much published research fails to replicate (Ioannidis, 2005).

10. Where the formulas come from : likelihood

So far every formula has been about means, because the CLT hands us means for free. But you will want intervals for medians, odds ratios, regression coefficients, model parameters. There is one principle underneath all of them.

Given data $x_1,\dots,x_n$ and a family of candidate distributions indexed by $\theta$, the likelihood is the probability of the data you actually saw, read as a function of $\theta$ : $$ L(\theta) = \prod_{i=1}^{n} f(x_i \mid \theta), \qquad \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n}\log f(x_i \mid \theta). $$ The maximum likelihood estimator $\hat\theta$ is the value making the observed data least surprising.

In plain words : you saw what you saw. Among all the possible worlds, pick the one in which what you saw was most expected. Applied to a normal population this returns the sample mean — so section 2 was a special case all along.

Two results make this the general engine of the subject. First, there is a floor on precision : no unbiased recipe, however clever, can wobble less than a certain amount.

Define the Fisher information $I(\theta) = -\mathbb{E}\left[\frac{\partial^{2}}{\partial\theta^{2}}\log f(X\mid\theta)\right]$, a measure of how sharply the likelihood peaks. Then any unbiased estimator $\hat\theta$ built from $n$ independent observations satisfies $$ \operatorname{Var}(\hat\theta) \;\ge\; \frac{1}{n\,I(\theta)} . $$

In plain words : a flat likelihood means the data barely distinguishes nearby values of $\theta$, and no amount of cleverness can extract precision that is not in the data. Information sets a hard limit (Cramér, 1946).

Second, the maximum likelihood estimator attains that floor, asymptotically, and is asymptotically normal :

\[\begin{equation} \sqrt{n}\,\big(\hat\theta - \theta\big) \;\xrightarrow{\;d\;}\; \mathcal{N}\!\left(0, \;\frac{1}{I(\theta)}\right). \label{eq:mle} \end{equation}\]

This is the licence for everything in Part II. It says that for any parameter estimated by maximum likelihood, with enough data, “estimate $\pm 1.96$ standard errors” is a valid 95 % interval — the same recipe, far beyond means.

If you have read the softmax and cross-entropy post, you have already done maximum likelihood without the name. Minimising cross-entropy is maximising the log-likelihood of a categorical model, and gradient descent is simply how you climb $\ell(\theta)$ when no closed form exists. Training a classifier and fitting a statistical model are the same operation.

11. When you cannot do the maths : the bootstrap

Sometimes there is no formula. What is the standard error of a median? Of the ratio of two correlations? Of the 90th percentile? Deriving these is hard or impossible — and since 1979 you have not needed to (Efron, 1979).

The idea is almost impudent. You do not have the population, but you do have the best available picture of it : your sample. So treat the sample as the population and resample from it, with replacement, thousands of times. The spread of the statistic across those resamples estimates the spread it would have had across real repeat samples.

import random, statistics

def bootstrap_ci(xs, stat=statistics.median, B=10000):
    """A 95 % interval for any statistic, no algebra."""
    n = len(xs)
    draws = (stat(random.choices(xs, k=n)) for _ in range(B))
    reps = sorted(draws)
    return reps[int(0.025 * B)], reps[int(0.975 * B)]

Four lines, and it works for essentially any statistic. It is not magic — it still needs independent observations, and it is unreliable for statistics that depend on the extreme tail (the maximum, for instance, since no resample can ever exceed the largest value you saw). But it converts an analytical problem into a computational one, and computation is cheap. The resamples are also completely independent of one another, which makes the bootstrap a textbook MapReduce job when $B$ needs to be large.

12. The assumptions that actually break things

Formulas fail quietly. These are the ones that matter, in the order they cause damage.

Independence. This is the big one, and $\eqref{eq:var}$ shows why : the derivation of $\sigma^2/n$ used independence to add the variances. If your observations are correlated — repeated measurements on the same users, sensor readings over time, pupils within the same classroom — then the effective sample size is much smaller than $n$. Your standard errors come out too small, your intervals too narrow, your p-values too impressive. Nothing in the output warns you. Clustered or time-series data needs methods that model the dependence.

Finite variance. The CLT requires $\sigma^2 < \infty$. Distributions with very heavy tails — some financial returns, city sizes, network traffic — can violate this or approach it so slowly that $n = 10{,}000$ is still not “enough”. A histogram and a look at the extreme values costs a minute and is worth it.

Representativeness, and why more data does not save you. Return to $\eqref{eq:mse}$ : $\operatorname{MSE} = \text{bias}^2 + \text{variance}$. Sample size only attacks the second term. If your sampling procedure is systematically skewed — you surveyed only the users who answered the phone, only the patients who came back, only the servers that did not crash — then collecting ten times the data shrinks the interval by a factor of three around the wrong number.

Bias does not shrink with $n$. Variance does. A biased sample plus a lot of data does not give you the truth — it gives you a narrow, confident interval that does not contain it.

That is the failure mode behind the famous polling disasters, and it is not fixable by any amount of arithmetic downstream. Go back and stir the pot.

13. A checklist for the future you

Reference material, in the order you will need it.

Before collecting anything.

  1. Write down the one question and the one number that answers it.
  2. Decide the smallest effect that would actually change a decision. Not the smallest detectable — the smallest useful.
  3. Run $\eqref{eq:power}$ with that effect to get the sample size. If the answer is unaffordable, say so now rather than running an underpowered study.
  4. Write down the analysis you will perform, before any data exists.

When the data is in.

  1. Plot it first. Always. Distribution shape, outliers, and impossible values are visible in ten seconds and invisible in a summary table.
  2. Check independence honestly. Ask what could make two observations resemble each other.
  3. Report the estimate, the confidence interval, and $n$. The interval is the deliverable ; the p-value is at best a footnote.

When reading someone else’s.

  1. Was the effect size reported, or only significance ?
  2. How many comparisons were actually made — including the ones not reported ?
  3. Was the analysis chosen before or after seeing the outcome ?
  4. For a null result : what was the power ? Without it, “no difference” means nothing.
  5. Where did the sample come from, and who could not have ended up in it ?

14. Conclusion

Strip away the vocabulary and this whole subject is one honest admission repeated in different registers : you saw a piece, not the whole, so anything you say must come with the size of your ignorance attached.

Estimation produces the guess. The standard error $\sigma/\sqrt{n}$ measures how much that guess wobbles, and its square root sets the price of every study ever run. The confidence interval reports the guess and the wobble together, and its 95 % is a promise about the method, never about the interval in your hand. The hypothesis test is the same interval read backwards, asking only whether one particular value survives — which is why the interval, carrying the effect and its uncertainty at once, is almost always the better thing to report. Underneath, maximum likelihood generalises the whole machinery beyond averages, and the bootstrap rescues you when the algebra runs out.

And one thing sits outside all of it. Every formula here assumes you stirred the pot. No amount of data, no correction, no cleverer estimator repairs a sample that was never representative — it only makes you more confident about the wrong answer. That is the part worth remembering when everything else has faded.

References

  1. Benjamini, Y., & Hochberg, Y. (1995). Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society, Series B, 57(1), 289–300.
    @article{Benjamini1995,
      author = {Benjamini, Yoav and Hochberg, Yosef},
      title = {Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing},
      journal = {Journal of the Royal Statistical Society, Series B},
      volume = {57},
      number = {1},
      pages = {289--300},
      year = {1995}
    }
    
  2. Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.). Lawrence Erlbaum Associates.
    @book{Cohen1988,
      author = {Cohen, Jacob},
      title = {Statistical Power Analysis for the Behavioral Sciences},
      edition = {2nd},
      year = {1988},
      publisher = {Lawrence Erlbaum Associates}
    }
    
  3. Cramér, H. (1946). Mathematical Methods of Statistics. Princeton University Press.
    @book{Cramer1946,
      author = {Cram{\'e}r, Harald},
      title = {Mathematical Methods of Statistics},
      year = {1946},
      publisher = {Princeton University Press}
    }
    
  4. Efron, B. (1979). Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics, 7(1), 1–26.
    @article{Efron1979,
      author = {Efron, Bradley},
      title = {Bootstrap Methods: Another Look at the Jackknife},
      journal = {The Annals of Statistics},
      volume = {7},
      number = {1},
      pages = {1--26},
      year = {1979}
    }
    
  5. Fisher, R. A. (1925). Statistical Methods for Research Workers. Oliver and Boyd.
    @book{Fisher1925,
      author = {Fisher, Ronald A.},
      title = {Statistical Methods for Research Workers},
      year = {1925},
      publisher = {Oliver and Boyd}
    }
    
  6. Gelman, A., & Loken, E. (2014). The Statistical Crisis in Science. American Scientist, 102(6), 460–465.
    @article{Gelman2014,
      author = {Gelman, Andrew and Loken, Eric},
      title = {The Statistical Crisis in Science},
      journal = {American Scientist},
      volume = {102},
      number = {6},
      pages = {460--465},
      year = {2014}
    }
    
  7. Ioannidis, J. P. A. (2005). Why Most Published Research Findings Are False. PLoS Medicine, 2(8), e124.
    @article{Ioannidis2005,
      author = {Ioannidis, John P. A.},
      title = {Why Most Published Research Findings Are False},
      journal = {PLoS Medicine},
      volume = {2},
      number = {8},
      pages = {e124},
      year = {2005}
    }
    
  8. Neyman, J. (1937). Outline of a Theory of Statistical Estimation Based on the Classical Theory of Probability. Philosophical Transactions of the Royal Society A, 236(767), 333–380.
    @article{Neyman1937,
      author = {Neyman, Jerzy},
      title = {Outline of a Theory of Statistical Estimation Based on the Classical Theory of Probability},
      journal = {Philosophical Transactions of the Royal Society A},
      volume = {236},
      number = {767},
      pages = {333--380},
      year = {1937}
    }
    
  9. Neyman, J., & Pearson, E. S. (1933). On the Problem of the Most Efficient Tests of Statistical Hypotheses. Philosophical Transactions of the Royal Society A, 231, 289–337.
    @article{NeymanPearson1933,
      author = {Neyman, Jerzy and Pearson, Egon S.},
      title = {On the Problem of the Most Efficient Tests of Statistical Hypotheses},
      journal = {Philosophical Transactions of the Royal Society A},
      volume = {231},
      pages = {289--337},
      year = {1933}
    }
    
  10. Wasserstein, R. L., & Lazar, N. A. (2016). The ASA Statement on p-Values: Context, Process, and Purpose. The American Statistician, 70(2), 129–133.
    @article{Wasserstein2016,
      author = {Wasserstein, Ronald L. and Lazar, Nicole A.},
      title = {The ASA Statement on p-Values: Context, Process, and Purpose},
      journal = {The American Statistician},
      volume = {70},
      number = {2},
      pages = {129--133},
      year = {2016}
    }
    
  11. Student. (1908). The Probable Error of a Mean. Biometrika, 6(1), 1–25.
    @article{Student1908,
      author = {{Student}},
      title = {The Probable Error of a Mean},
      journal = {Biometrika},
      volume = {6},
      number = {1},
      pages = {1--25},
      year = {1908}
    }
    
  1. The threshold $\alpha = 0.05$ has no mathematical standing whatsoever. Fisher suggested it as a convenience — “it is convenient to take this point as a limit in judging whether a deviation is to be considered significant or not” — and a convention hardened around it. Fields with more at stake pick their own : particle physics requires roughly $5\sigma$, about $p < 3 \times 10^{-7}$, precisely because it runs an enormous number of comparisons. ↩

who am i

Hi! I am a Data Scientist by profession, an Emacs devotee and an untalented bassist. I intend to use this space for writing about things that I think I have understood well in the hope that they may be helpful to others, including my future self.

what is this

OhMyKode is an opportunity to share knowledge about mathematics, computer science, machine learning and algorithmic beauty, which allows us to improve our skills and learn in depth. It is a sharing place to learn the how and the why.

© MMXVIII - MMXXVI by Maâmra Youcef - معامره يوسف
Content available under Creative Commons (BY-NC-SA) unless otherwise noted.
This site is hosted at Github Pages and powered by Jekyll & Papyrus.
“We can't skip Math forever !”