02 — Statistics & Inference

Probability starts from a known distribution and predicts what data will look like. Statistics runs that arrow backwards: you are handed the data and must recover something about the process that produced it.

The recovery takes three forms: a best guess, an honest interval around that guess, and a ship/do-not-ship decision with a stated error rate. A fourth thing deserves equal attention — what has to be true of your data before any of the three mean anything.

By the end you should be able to derive the sample size an A/B test needs, state what a p-value is and is not without hedging, name the assumption behind every test you run, and recognise the handful of ways real experiments quietly violate those assumptions.

The applied half is A/B testing: splitting live traffic between an unchanged “control” experience and a changed “treatment” one, then measuring the difference. It is the most-asked topic in any data-facing interview, because it is where these ideas break in production.

Peeking at results early, testing many metrics at once, randomising at the wrong level, and Simpson’s paradox are all versions of the same failure: a procedure advertises an error rate that assumed something the analyst violated.

What goes in and what comes out

Fix the shape of the problem first, because everything below is a way of getting from the input to the output.

The input is always a finite list of numbers drawn from a process you cannot observe directly. It might be five measurements, [2, 4, 4, 4, 6]. It might be two hundred thousand rows of user_id, variant, converted. What makes it statistics rather than arithmetic is that you want to say something about the process, and all you have is one sample from it.

The output is three artefacts, and it is worth being able to name which one a stakeholder is actually asking for:

1. A point estimate — a single number that is your best guess at the unknown quantity. “The conversion rate is 5.2%.”
2. An interval — a range expressing how much that guess would move if you re-ran the whole experiment on fresh data. “5.2%, give or take 0.3 points.”
3. A decision — ship or do not ship, carrying an explicit rate at which the decision is wrong. “The treatment is better, and if it truly were not, we would only have said this 5% of the time.”

Everything in this chapter is either a way of producing one of those three, or a way of noticing when that machinery is misleading you.

Notation, read aloud once

Below is every symbol used in this chapter, with how to say it and what it means. You do not have to memorise it, but you should have somewhere to look up any symbol you hit later.

The single most important pattern in the table: Roman letters are things you compute from data, Greek letters are the unknown truths you are chasing. xbar and s are computed from data. mu and sigma are properties of a process you never observe directly.

n            "n"          how many observations are in the sample
x_1 ... x_n  "x one" ...  the observations themselves
xbar         "x-bar"      the sample mean: add the observations, divide by n
s            "s"          the sample standard deviation: typical distance
                          of an observation from xbar
s^2          "s squared"  the sample variance, the square of the above
mu           "mu"         the TRUE mean of the process; unknown, fixed
sigma^2      "sigma squared"  the TRUE variance of the process; unknown, fixed
sigma        "sigma"      the true standard deviation, the square root of sigma^2
theta        "theta"      a stand-in for whatever unknown number you are after
theta_hat    "theta-hat"  an estimate of it; a hat ALWAYS means "computed from data"
E[ · ]       "expected value of"   the long-run average over all possible samples
Var( · )     "variance of"         the average squared distance from that average
Cov( ·, · )  "covariance of"       how much two quantities move together
sum_i        "sum over i"          add the term once for every observation
sqrt( · )    "square root of"
~            "is distributed as"
->           "gives" or "converges to", depending on context
|a|          "absolute value of a", the size of a ignoring its sign

Two collisions to know about before they trip you up. theta is used in §1 and §2 as a generic name for “the unknown thing”, and again in §6 as the name of a specific tuning constant in CUPED; both are standard and the section says which is meant. And a standalone t in §3 onward always means Student’s t distribution, never theta.

Greek letters that appear later — alpha, beta, delta, rho — are each glossed at the point they are first used.

1. Estimators, and the three properties that matter

Every point estimate comes from some recipe for turning data into a number, and recipes can be compared. Three properties do the comparing — and “unbiased”, the one everybody leads with, is a far weaker recommendation than it sounds.

An estimator is a random variable

An estimator is any recipe that turns a sample into a number.

“Average the rows” is an estimator. “Take the middle row” is an estimator. “Ignore the data and always report 7” is also an estimator, just a bad one.

The sample is random: a different week of traffic hands you different rows. So the number the recipe produces is itself random. That makes an estimator a random variable — a quantity whose value depends on chance, and which therefore has a whole distribution of possible values rather than one fixed value.

The distribution of answers the recipe would give across all the samples you might have drawn is called the sampling distribution. It is the central object of this chapter. Every property below is a property of that distribution, not of the one number you happened to compute.

The three properties

In words, before symbols.

Bias is the systematic offset. Averaged over all the samples you might have drawn, does the recipe land on the truth or beside it?

Variance is how much the recipe moves when you swap one sample for another — its spread across samples.

Mean squared error, abbreviated MSE, is the average squared distance between the estimate and the truth. It is the only one of the three that answers the real question, because it charges you both for being wrong on average and for being erratic.

Here are the same three in symbols, plus a fourth property that is about what happens as the sample grows.

bias(theta_hat)   = E[theta_hat] - theta               systematic offset
Var(theta_hat)    = E[(theta_hat - E[theta_hat])^2]    sensitivity to the draw
MSE(theta_hat)    = bias^2 + Var                       the thing you actually care about
consistency       theta_hat -> theta in probability as n -> infinity

In words: the bias of theta-hat is the expected value of theta-hat minus the true theta. The variance of theta-hat is the expected squared distance between theta-hat and its own average. The mean squared error is the bias squared plus the variance.

Consistency means theta-hat converges to theta in probability as n grows without bound. “In probability” is the precise way of saying: name any tolerance you like, however small, and the chance of the estimate falling outside that tolerance goes to zero as the sample grows.

The diagram below shows how the four fit together. One estimator, viewed as a random variable, splits into a bias term and a variance term; the two recombine into MSE; consistency hangs off to the side as a separate question. The two orange boxes are the counterexamples explained just below.

flowchart TD
    E["Estimator theta_hat<br/>a random variable"] --> B["Bias<br/>E[theta_hat] - theta"]
    E --> V["Variance<br/>spread across samples"]
    B --> M["MSE = bias^2 + Var<br/>the only thing that matters"]
    V --> M
    E --> C["Consistency<br/>the estimate converges as n grows"]
    C --> N1["unbiased does NOT imply consistent<br/>(always use X_1: unbiased, never converges)"]
    C --> N2["consistent does NOT imply unbiased<br/>(MLE variance: biased at every n)"]

    style M fill:#2d6a4f,color:#fff
    style N1 fill:#bc6c25,color:#fff
    style N2 fill:#bc6c25,color:#fff

MSE is the only thing that matters when you have to choose between two estimators, because it is the sum and not a part.

That split — expected squared error equals squared bias plus variance — is the same algebra used in Total expectation and total variance for prediction error, applied to an estimator instead of a prediction. It is also the reason unbiasedness is oversold: unbiasedness zeroes one of the two terms and says nothing at all about the other.

Unbiased and consistent are independent properties

Both counterexamples in the diagram get asked in interviews, and each is one sentence.

• Unbiased but not consistent. Ignore the sample and report the first observation, X_1. Its expected value is exactly mu, the true mean, so it is perfectly unbiased at every sample size. It also never converges to anything, because collecting a million rows and reading only the first one is not learning.
• Consistent but biased. The maximum-likelihood estimator for variance, derived in the next section, sits below the true variance at every finite n, but the gap shrinks to nothing as n grows.

Neither property implies the other, and neither on its own tells you whether an estimator is good — only MSE does that.

A biased estimator that beats the unbiased one

The claim “MSE, not bias” is easy to accept in the abstract. Here it is in numbers you can check by hand.

Suppose the true mean is mu = 10, the true standard deviation is sigma = 10, and you collect n = 4 observations. Consider shrinking the sample mean toward zero by a factor c: the estimator is c · xbar.

E[c·xbar]   = c·mu = 10c        so  bias = 10c - 10 = 10(c - 1)
Var(c·xbar) = c^2·sigma^2/n = c^2·(100/4) = 25·c^2

MSE(c)      = 100·(c - 1)^2 + 25·c^2

c = 1.0  (the plain sample mean, unbiased):  100·0    + 25·1.00  = 25.0
c = 0.8:                                     100·0.04 + 25·0.64  =  4 + 16 = 20.0

The unbiased estimator scores 25. Shrinking it by 20% — which introduces a bias of 10(0.8 - 1) = -2 — scores 20. The biased estimator is 20% better on the loss you actually care about.

Differentiating MSE(c) and setting it to zero gives the best possible c:

d/dc = 200·(c - 1) + 50·c = 0   ->   250·c = 200   ->   c* = 0.8

So c = 0.8 is not a lucky guess, it is optimal. Nobody uses it, because c* works out to mu^2/(mu^2 + sigma^2/n), which depends on mu — the very thing you were trying to estimate.

Two lessons, and the second one comes back in §2. Unbiasedness is a constraint, and dropping it can lower the loss you actually care about. And the estimator that minimises MSE is usually a different estimator from the unbiased one, which is why §2 finds a third answer to “what should the variance divisor be”.

2. Maximum likelihood, worked twice

Judging estimators is one thing; where do they come from when you do not already have one? There is one general-purpose recipe, and applying it to the two distributions that come up most also settles where the mysterious n - 1 in the variance formula comes from.

Maximum likelihood estimation, universally abbreviated MLE, answers one question: which value of the unknown parameter would have made the data I actually saw as probable as possible?

The likelihood is that probability, written as a function of the parameter rather than of the data. You freeze the observations at the values you got and slide the parameter around, watching the probability rise and fall. The value that maximises it is the MLE.

The recipe is mechanical, five steps:

1. Write the probability of the data as a function of the parameter.
2. Take the logarithm. This turns a product into a sum, and it does not move the location of the maximum, because the logarithm is strictly increasing — whatever value of the parameter maximised the product also maximises its log.
3. Differentiate with respect to the parameter.
4. Set the derivative to zero and solve.
5. Check the second derivative is negative, so you found a maximum rather than a minimum.

Bernoulli

A Bernoulli trial is a single yes/no event with a fixed success probability p — a coin flip, or a visitor who either converts or does not. Suppose you observe n trials and k successes, and want to estimate p.

Steps 1 through 4 of the recipe, one per line:

L(p)      = p^k (1-p)^(n-k)
log L     = k·log p + (n-k)·log(1-p)
d/dp      = k/p - (n-k)/(1-p) = 0
          -> k(1-p) = (n-k)p  ->  k - kp = np - kp  ->  p_hat = k/n

Line by line:

• L(p), read “L of p”, is the likelihood: the probability of seeing exactly these k successes and n - k failures. For independent trials that is p multiplied by itself k times, times 1 - p multiplied by itself n - k times.
• Taking the logarithm turns that product into a sum, because log(a·b) = log a + log b.
• The third line is d/dp, read “the derivative with respect to p”, set equal to zero. The derivative of k·log p is k/p; the derivative of (n-k)·log(1-p) is -(n-k)/(1-p), with the minus sign coming from the chain rule on 1 - p.
• The fourth line multiplies out and cancels. k(1-p) = (n-k)p becomes k - kp = np - kp; the -kp appears on both sides and cancels, leaving k = np, so p_hat = k/n — read “p-hat equals k over n”.

Step 5 is the curvature check:

d2/dp2    = -k/p^2 - (n-k)/(1-p)^2  < 0            a maximum

The second derivative, read “d two by d p squared”, is a sum of two negative terms, so it is negative everywhere. That confirms the point where the first derivative vanishes is a maximum, not a minimum.

The answer is the sample proportion — 37 conversions out of 500 visitors gives p_hat = 0.074 — which is exactly what anyone would have guessed. That is the point: the “obvious” estimator is obvious because it is the MLE, not the other way round.

Gaussian, and where n-1 comes from

A Gaussian distribution, also called the normal distribution, is the familiar symmetric bell curve. It is pinned down entirely by two numbers: its mean mu and its variance sigma^2. Apply the same five-step recipe to n independent Gaussian observations and you get estimates for both parameters at once.

The block below starts from the log-likelihood and differentiates it twice — once with respect to mu, once with respect to sigma^2 — treating the other parameter as fixed each time.

log L = -(n/2)·log(2·pi·sigma^2) - sum_i (x_i - mu)^2 / (2·sigma^2)

d/dmu:      sum_i (x_i - mu) / sigma^2 = 0            ->  mu_hat = xbar
d/dsigma^2: -n/(2·sigma^2) + sum_i (x_i - mu)^2 / (2·sigma^4) = 0
                                                      ->  sigma^2_hat = (1/n) sum_i (x_i - xbar)^2

The first line is not magic. One Gaussian observation has probability density (1/sqrt(2·pi·sigma^2))·exp(-(x - mu)^2/(2·sigma^2)). Multiply n of those together for n independent observations, take the log, and the 1/sqrt(...) factors become -(n/2)·log(2·pi·sigma^2) while the exponents become the sum on the right. Nothing else happens.

Now the two derivatives.

Setting d/dmu to zero says sum_i (x_i - mu) = 0 — the deviations must add up to zero. That happens at exactly one value of mu, the sample mean, so mu_hat = xbar.

Setting d/dsigma^2 to zero and multiplying through by 2·sigma^4 gives -n·sigma^2 + sum_i (x_i - mu)^2 = 0, so sigma^2_hat is the average squared deviation from the fitted mean. Maximum likelihood tells you to divide the sum of squared deviations by n.

Why the divisor is n-1

Every statistics course then tells you to divide by n - 1 instead. Here is why, from the definition.

The obstacle is that you want to say something about spread around the true mean mu, but all you can measure is spread around the sample mean xbar. The trick is to write one in terms of the other, by adding and subtracting mu inside the bracket.

sum_i (x_i - xbar)^2 = sum_i [ (x_i - mu) - (xbar - mu) ]^2
                     = sum_i (x_i - mu)^2 - 2(xbar - mu)·sum_i (x_i - mu) + n(xbar - mu)^2

and  sum_i (x_i - mu) = n(xbar - mu),  so the middle term is -2n(xbar - mu)^2:

                     = sum_i (x_i - mu)^2 - n(xbar - mu)^2

take expectations:
E[ sum_i (x_i - xbar)^2 ] = n·sigma^2 - n·Var(xbar) = n·sigma^2 - n·(sigma^2/n)
                          = (n - 1)·sigma^2

Three steps in that block are worth slowing down on, because they are where people lose the thread.

The expansion. The second line is just (a - b)^2 = a^2 - 2ab + b^2 with a = (x_i - mu) and b = (xbar - mu), summed over i. The term b does not depend on i, so it comes out of the sum: sum_i 2ab becomes 2(xbar - mu)·sum_i (x_i - mu), and sum_i b^2 becomes n(xbar - mu)^2.

The middle term. sum_i (x_i - mu) = sum_i x_i - n·mu = n·xbar - n·mu = n(xbar - mu). Substituting that in makes the middle term -2n(xbar - mu)^2, which combines with the +n(xbar - mu)^2 on the end to leave -n(xbar - mu)^2.

The expectations. Each (x_i - mu)^2 has expected value sigma^2 — that is the definition of variance — so the first sum contributes n·sigma^2. And (xbar - mu)^2 has expected value Var(xbar) = sigma^2/n, because averaging n independent draws divides their variance by n. So n·sigma^2 - n·(sigma^2/n) = n·sigma^2 - sigma^2 = (n-1)·sigma^2.

The conclusion is E[sigma^2_MLE] = ((n-1)/n)·sigma^2. The maximum-likelihood variance is biased low by the factor (n-1)/n. Dividing by n - 1 instead of n corrects it exactly, at every sample size, not just asymptotically.

The mechanism in one sentence: xbar is by construction the value of c that minimizes sum (x_i - c)^2, so measuring spread around it necessarily understates spread around the true mu — by exactly Var(xbar) = sigma^2/n per observation, which is the one degree of freedom you spent estimating the center.

A degree of freedom is a piece of information the data spent on fitting something rather than on measuring spread. You fitted one number, the centre, so you have n - 1 pieces left.

The arithmetic on five numbers

Take x = [2, 4, 4, 4, 6] and follow it through.

The sample mean is (2+4+4+4+6)/5 = 20/5 = 4. The deviations from it are [-2, 0, 0, 0, 2]. The sum of their squares — written SS for sum of squares — is 4 + 0 + 0 + 0 + 4 = 8.

xbar = 4      deviations [-2, 0, 0, 0, 2]      SS = 8

MLE      (/n)     = 8/5 = 1.60          biased low by factor 4/5
unbiased (/(n-1)) = 8/4 = 2.00

Dividing by n = 5 gives 1.60. Dividing by n - 1 = 4 gives 2.00. The ratio is 1.60/2.00 = 0.8 = 4/5, which is (n-1)/n exactly as the derivation predicted.

Two facts that separate a good answer from a great one

Unbiasedness does not survive a nonlinear transform. The sample variance s^2 is unbiased for sigma^2, but the sample standard deviation s is not unbiased for sigma.

Taking a square root is a concave operation. Jensen’s inequality — for a concave function, the function of the average is at least the average of the function — forces E[s] < sigma. So the “unbiased standard deviation” everybody quotes is not unbiased. In the five-number example, s = sqrt(2.00) = 1.414, and on average across resamples that figure sits below the true sigma.

The correction is unnecessary for the sample mean. xbar is already unbiased at any n, because averaging is a linear operation and expectation passes straight through linear operations. Only the nonlinear step — squaring, then dividing — created the problem in the first place.

Unbiased is not the same as best

Unbiasedness is a constraint, and constraints cost something. Ask instead which divisor minimises MSE, and a third answer appears.

Consider the whole family of estimators T_c = c · sum_i (x_i - xbar)^2, read “T sub c”, where c is any constant you are free to choose. Setting c = 1/n gives the MLE; setting c = 1/(n-1) gives the unbiased estimator. The question is which c minimises MSE.

To answer it you need the mean and the variance of SS = sum_i (x_i - xbar)^2. Under normality the rescaled quantity SS/sigma^2 follows a chi-square distribution with n - 1 degrees of freedom. A chi-square is the distribution of a sum of squared independent standard normal variables, which is exactly what the deviations are once you divide them by sigma. That distribution hands you both moments:

E[SS]   = (n-1)·sigma^2
Var(SS) = 2(n-1)·sigma^4

Now assemble the MSE. Since T_c = c·SS, expectation and variance scale in the usual way — E[c·SS] = c·E[SS] and Var(c·SS) = c^2·Var(SS):

bias(T_c) = c(n-1)·sigma^2 - sigma^2 = sigma^2·( c(n-1) - 1 )
bias^2    = sigma^4·( c(n-1) - 1 )^2
Var(T_c)  = c^2·2(n-1)·sigma^4

MSE(T_c)/sigma^4 = ( c(n-1) - 1 )^2  +  2·c^2·(n-1)

Dividing through by sigma^4 on the last line is deliberate: it makes the answer independent of the unknown scale, so the best c does not depend on a quantity you do not know. (Contrast that with the shrinkage example in §1, where the optimal c did depend on mu.)

Differentiate with respect to c and set to zero:

d/dc = 2(n-1)( c(n-1) - 1 ) + 4c(n-1) = 0

divide both sides by 2(n-1):
       ( c(n-1) - 1 ) + 2c = 0
       cn - c - 1 + 2c = 0
       cn + c = 1
       c·(n+1) = 1        ->    c* = 1/(n+1)

Read c* as “c-star”. The minimum-MSE estimator divides by n + 1, not n - 1 and not n.

On the five-number sample from above, that is SS/(n+1) = 8/6 = 1.33, sitting below both the MLE’s 1.60 and the unbiased 2.00.

Nobody uses it. Unbiasedness composes nicely through downstream algebra — the average of unbiased estimates is still unbiased, which is not true of minimum-MSE ones — and a small bias is easier to reason about than a small MSE gain. But knowing that the “correct” divisor depends on which loss you chose is what shows you understand that n - 1 is a choice, not a law.

What maximum likelihood assumes

Every result above is conditional. The table lists the four conditions, what each one means, and the exact form each takes when it breaks in a real product. Read the right-hand column first: those are the situations you will actually meet.

The n - 1 correction in particular is exact only under independence. If your rows are correlated, dividing by n - 1 does not make the variance estimate unbiased — it makes it wrong by a smaller-looking amount, which is the bug described in Randomization unit and the variance bug it causes.

3. Confidence intervals, and the interpretation nobody gets right

A point estimate alone is incomplete: it names a number without saying how far that number would move on fresh data. The interval that says so is also the most misread object in statistics, and the two standard misreadings cause real damage.

A confidence interval — abbreviated CI — is a range computed from the data that is designed to contain the unknown parameter a specified fraction of the time. For a mean it takes this form:

95% CI  =  xbar  +/-  t_{0.975, n-1} · s / sqrt(n)

Read aloud: the interval is the sample mean, plus and minus a multiplier times s over the square root of n. Three pieces need naming.

s is the sample standard deviation: the typical distance of one observation from xbar.

s / sqrt(n) is the standard error, abbreviated SE. It is the standard deviation of the estimate itself across repeated samples. s describes how spread out individual rows are; SE describes how much xbar would jump around if you re-ran the whole study. The two are constantly confused, and keeping them apart is what this section is for. Note that SE shrinks as n grows while s does not — more data does not make your users more alike, it makes your estimate of their average steadier.

t_{0.975, n-1}, read “t at the 97.5th percentile with n minus 1 degrees of freedom”, is the multiplier. It comes from Student’s t distribution: a bell curve with heavier tails than the Gaussian, used because you estimated s from the same data instead of knowing it, and the extra tail weight pays for that uncertainty. It is 0.975 rather than 0.95 because the 5% of missed cases is split between the two tails, 2.5% each.

The interpretation, and the two misreadings

Correct interpretation: the procedure produces intervals that contain the true parameter in 95% of repeated experiments. It is a statement about the method’s long-run coverage, not about this interval.

The common wrong one: “there is a 95% probability the true value lies in [50.0, 54.0].”

Under frequentist assumptions the parameter is a fixed unknown number and the interval is what is random, so that probability is 0 or 1 — you just do not know which. (“Frequentist” means the school of statistics that treats probability as long-run frequency over repeated experiments; the parameter has no distribution because it is not random.)

A Bayesian credible interval does support that statement, because Bayesian methods treat the parameter as having a distribution. But it requires a prior — a distribution over the parameter that you supply before seeing the data — and the two kinds of interval coincide only under specific priors.

The second wrong one, which does real damage: confusing a confidence interval with a prediction interval. A prediction interval, abbreviated PI, is a range for where the next single observation will fall, not for where the mean is.

The two are worth computing side by side, because the gap between them is far larger than most people expect. Take n = 100 observations with sample mean xbar = 52 and sample standard deviation s = 10. Then SE = 10 / sqrt(100) = 10/10 = 1.0, and the multiplier is t_{0.975, 99} = 1.984.

CI for the MEAN:            52 +/- 1.984 × 1.0                    = [50.02, 53.98]   width  3.97
PI for a NEW OBSERVATION:   52 +/- 1.984 × 10 × sqrt(1 + 1/100)   = [32.06, 71.94]   width 39.9

Check the second line by hand: sqrt(1 + 0.01) = 1.005, so the half-width is 1.984 × 10 × 1.005 = 19.94, giving 52 - 19.94 = 32.06 and 52 + 19.94 = 71.94.

The prediction interval uses s rather than s / sqrt(n). That single change is what does the damage: the mean of 100 observations is ten times steadier than one observation. The extra sqrt(1 + 1/n) factor then adds the uncertainty about where the mean sits on top of the spread of individuals.

Ten times wider, because the CI describes uncertainty about a mean and the PI must also absorb the spread of individuals. When a stakeholder asks “so a user will convert between 50% and 54%?” they are asking for the second interval and you have handed them the first. That confusion is behind most over-confident forecasts shipped to a business audience.

A confidence interval and a hypothesis test are the same object

The 95% CI is exactly the set of null values that a two-sided test at alpha = 0.05 would fail to reject.

Two terms in that sentence get defined properly in Hypothesis testing: alpha, read “alpha”, is the error rate you agree to tolerate, and “two-sided” means you would be surprised by a difference in either direction.

The practical consequence: if the interval excludes zero, then p < 0.05 for the corresponding test. You never need to compute both.

The word corresponding is doing real work, though. The duality is exact only when the interval and the test are built from the same pivot — a quantity whose distribution does not depend on the unknown parameter, which is what lets a single cutoff be valid for every possible truth. Build them from different pivots and they can disagree.

The classic disagreement is on a proportion near 0 or 1. Three different pivots exist for the same question:

• A Wald interval takes the estimate plus or minus a multiple of its standard error.
• A likelihood-ratio test compares how well the data fit the alternative against how well they fit the null.
• A score test uses the slope of the log-likelihood at the null instead.

All three answer “is p equal to this value?” and they agree only in large samples. Near the boundary they can give different verdicts on the same data.

The same mismatch appears between a robust-standard-error interval, which estimates the standard error from the observed spread rather than from an assumed model, and a permutation test, which builds the null distribution by repeatedly reshuffling the group labels and recomputing the statistic.

What a confidence interval assumes

The interval is only as honest as four conditions. The middle column says which part of the formula would break; the right-hand column is the version of that break you will actually meet in an A/B test.

The normality assumption is about the mean, not about the data. Individual conversions are 0/1 and could not be less bell-shaped; that is fine. What must be approximately Gaussian is xbar, and for a rare binary event the usual rule of thumb is at least about 10 successes and 10 failures per arm before you trust the approximation.

4. Hypothesis testing

Estimates and intervals describe; at some point somebody has to decide. The decision half of statistics collapses the data into ship or do-not-ship, with an error rate agreed on in advance.

A hypothesis test is a decision rule with a controlled error rate. It has three moving parts.

The null hypothesis, written H0 and read “H-nought”, is the boring default you are trying to disprove. Usually: “the treatment changed nothing.”

The alternative hypothesis, written H1, is what you would conclude instead. Usually: “the treatment changed something.”

The test statistic is a single number summarising the data, chosen so that you know its distribution if H0 were true. You compute it, then ask how surprising that value would be in a world where H0 holds.

The p-value, stated precisely

The p-value is the probability, computed assuming H0 is true, of observing a test statistic at least as extreme as the one observed — under the sampling and stopping rules you actually specified in advance.

Every clause in that sentence is load-bearing.

“Assuming H0 is true” means the whole calculation lives inside a hypothetical world where there is no effect. It never leaves that world.

“At least as extreme” means you count not only the result you got but everything further from the null. That is why the p-value depends on what you decided counts as extreme before you looked.

“The stopping rules you specified in advance” means the number is only valid for the experiment you committed to running — not the one you improvised once the numbers came in.

Three things it is not. Each row is a sentence you will hear in a real meeting, followed by what is wrong with it.

The first row is the one that matters most, so read both quantities in words and put them side by side.

P(H0 | data), read “the probability of H-nought given the data”, is the probability the null hypothesis is true now that you have seen the evidence. That is a posterior — a belief updated by data — and computing one requires a prior.

P(data-or-worse | H0), read “the probability of data this extreme or worse given that H-nought is true”, is the p-value. Same two symbols, opposite sides of the bar.

The p-value is not the probability the hypothesis is true, and no arithmetic on it alone will produce that probability.

What bridges the two is Bayes’ theorem from Bayes theorem and the question everyone gets wrong: the probability of a hypothesis given the evidence equals the probability of the evidence given the hypothesis, scaled by how plausible the hypothesis was to begin with.

Put numbers on it. Suppose 10% of product ideas genuinely work — an optimistic hit rate — so the prior is P(H1) = 0.10. You run every test at alpha = 0.05 with 80% power, meaning a genuine effect is detected 80% of the time. Now ask: of the results that come back significant, what fraction are real?

P(H1 | significant) = (0.10)(0.80) / [ (0.10)(0.80) + (0.90)(0.05) ]
                    = 0.080 / (0.080 + 0.045)  =  0.080 / 0.125  =  0.64

The numerator counts true findings: 10% of ideas work, and 80% of those are caught, so 0.10 × 0.80 = 0.080 of all ideas end up as real wins.

The denominator adds the false ones: 90% of ideas do not work, and 5% of those pass anyway, so 0.90 × 0.05 = 0.045 of all ideas end up as fake wins. Total significant results: 0.080 + 0.045 = 0.125.

The ratio 0.080 / 0.125 = 0.64 is the fraction of your “wins” that are real.

A “statistically significant at p < 0.05” result is wrong 36% of the time under those assumptions, and the p-value knows nothing about it.

Turn the dials and watch what moves:

The p-value is fixed at 0.05 in all three rows. Only the prior and the power moved. This is why underpowered studies are worse than useless: low power does not just miss real effects, it makes the effects you do find less likely to be real.

Errors, power, and the four cells

The null is either true or false, and you either reject it or you do not. Two choices times two truths gives exactly four outcomes. The two on the diagonal are correct; the two off it are the mistakes, and each has a name and a rate.

In plain words:

• A type I error is a false alarm. You declare an effect when there is none. Its rate is alpha, also called the significance level.
• A type II error is a miss. A real effect goes undetected. Its rate is beta, read “beta”.
• Power is 1 - beta: the probability of catching an effect that is genuinely there.

Note that the bottom-left cell says “fail to reject”, not “accept”. A test that does not reject has not shown the null is true. It has shown only that the data were not surprising enough to rule it out.

The asymmetry between the two dials matters. alpha is something you set — you pick 0.05 and that is that. beta is a consequence of four things: alpha, the sample size, the noise, and the true effect size.

You cannot pin down beta without committing to an effect size worth detecting: the smallest difference that would actually change your decision. In experimentation that committed number is called the minimum detectable effect, abbreviated MDE.

The sample-size formula, derived

Everything above becomes actionable once you can compute how many users a test needs.

Set the scene. Two groups of n users each. Each user’s outcome has variance sigma^2. The true difference between the groups is delta, read “delta”.

Your estimate of that difference, written d_hat and read “d-hat”, is the difference of the two group means. The variance of a difference of independent quantities is the sum of their variances, so Var(d_hat) = sigma^2/n + sigma^2/n = 2·sigma^2/n, and the standard error is its square root:

SE = sqrt( 2·sigma^2 / n )

Now the derivation. Five lines, each explained underneath.

Under H0:  d_hat ~ Normal(0, SE^2).      Reject when |d_hat|/SE > z_{1-alpha/2}.
Under H1:  d_hat ~ Normal(delta, SE^2).

Power = P( d_hat/SE > z_{1-alpha/2} | true effect delta )
      = P( Z > z_{1-alpha/2} - delta/SE )

For power 1 - beta we need     delta/SE - z_{1-alpha/2}  >=  z_{1-beta}

      delta  >=  ( z_{1-alpha/2} + z_{1-beta} ) · sqrt( 2·sigma^2 / n )

      n  =  2·sigma^2·( z_{1-alpha/2} + z_{1-beta} )^2 / delta^2      per arm

Line 1. If the null is true, the estimated difference is a bell curve centred on zero. You reject when its size in standard errors exceeds a cutoff z_{1-alpha/2}, read “z at one minus alpha over two” — the point on the standard Normal curve with alpha/2 of the probability beyond it in each tail. At alpha = 0.05 that cutoff is 1.960.

Line 2. If the alternative is true, the same bell curve is centred on delta instead of on zero. Nothing else about it changes; the whole curve slides right by delta/SE standard errors.

Line 3. Power is the chance that shifted curve still lands past the cutoff. To turn it into something you can look up, subtract the centre: if d_hat ~ Normal(delta, SE^2) then Z = (d_hat - delta)/SE is a standard Normal. So d_hat/SE > z_{1-alpha/2} is the same event as Z > z_{1-alpha/2} - delta/SE. That substitution is the only real step in the derivation.

Line 4. You want that probability to be at least 1 - beta. A standard Normal exceeds z_{1-beta} with probability exactly beta, and exceeds anything smaller with higher probability, so the requirement is z_{1-alpha/2} - delta/SE <= -z_{1-beta}, which rearranges to the line shown.

Line 5. Substitute SE = sqrt(2·sigma^2/n), square both sides, and solve for n.

The number worth memorising

Plug in the standard settings, alpha = 0.05 and power = 0.80:

z_{1-alpha/2} = z_{0.975} = 1.960
z_{1-beta}    = z_{0.80}  = 0.8416
sum                       = 2.8016
squared                   = 7.849

Since 7.849 × 2 = 15.70, which rounds to 16, the formula collapses to n ≈ 16·sigma^2/delta^2 per arm. That is the version to carry in your head.

The rounding is for your memory, not for your planning doc: the table below is computed with the full-precision (z_a + z_b)^2 = 7.84888, not with 7.849 and not with 16. Substituting 16 would inflate every row by about 2%.

Working it on a conversion rate

Take a baseline conversion rate of p1 = 0.05. A “relative lift” of 20% means the rate moves from 5% to 6%. So the absolute effect is delta = p1 × relative lift = 0.05 × 0.20 = 0.0100, and the treatment arm sits at p2 = p1 + delta = 0.06.

Here is the trap that makes most sample-size tables quietly under-deliver.

For a yes/no outcome, the per-unit variance is p(1-p). That is a function of the rate, so the two arms cannot have the same variance whenever the effect is real. Write q = 1 - p. The control arm sits at p1·q1 = 0.05 × 0.95 = 0.0475. The treatment arm, under H1, sits at p2·q2 = 0.06 × 0.94 = 0.0564.

The lazy formula uses 2·sigma^2 = 2 × 0.0475 = 0.0950, plugging in the smaller of the two variances twice. The honest one uses their sum, 0.0475 + 0.0564 = 0.1039:

n  =  ( p1·q1 + p2·q2 )·( z_{1-alpha/2} + z_{1-beta} )^2 / delta^2      per arm

where q = 1 - p. Equivalently, plug the pooled rate pbar = (p1 + p2)/2 into
n = 2·pbar·qbar·(z_a + z_b)^2 / delta^2, which lands within a handful of users.

Substitute the 20%-lift numbers all the way through:

n = 0.1039 × 7.84888 / (0.0100)^2
  = 0.81550 / 0.0001
  = 8,155 per arm        (the lazy version gives 0.0950 × 7.84888 / 0.0001 = 7,457)

The same arithmetic at four smaller lifts gives the rest of the table. Notice how fast the right-hand column grows as the lift you want to detect shrinks.

Say the consequence out loud, because it is the whole reason to care: a test run to 7,457 per arm instead of 8,155 has 76.4% power. It delivers 76% where you asked for 80%.

The shortfall narrows as the lift narrows, because p2 approaches p1 and the two variances converge. At a 10% lift the lazy formula asks for 29,826 instead of 31,231; at 5%, 119,303 instead of 122,121; at 2%, 745,644 instead of 752,700; at 1%, 2,982,575 instead of 2,996,695.

Type I error is unaffected in either version. Only the power falls short, which is exactly the kind of failure that does not announce itself. This is the “finite, roughly equal variance” row of the assumption table below, applied to your own planning arithmetic rather than to someone else’s test.

The one scaling law that runs the roadmap

n scales with 1/delta^2, so halving the effect you want to detect quadruples the traffic.

Check it against the table: going from a 10% lift to a 5% lift halves delta from 0.0050 to 0.0025, and n goes from 31,231 to 122,121. That is a factor of 3.91, just under the clean 4, because the variance term p1·q1 + p2·q2 also shrinks slightly as p2 moves back toward p1.

Put that in weeks. At 100k users per week, the 5% lift needs 244,242 users, so 2.4 weeks. The 1% lift needs 5,993,390 users, so about 60 weeks — a year and change.

The answer to “can we detect a 1% lift?” is almost never “yes”. It is “not with this metric,” which is why variance reduction (Cuped derived) is the main lever an experimentation platform can offer.

What the two-sample test assumes

The formula above is a chain of assumptions, and each link has a well-known way of snapping in an A/B test. Read the table one row at a time: what must hold, how a real experiment breaks it, and what you would see and do about it. The last column is the only one you can act on.

None of these violations produce an error message. The analysis runs, the number renders, the decision is made. That is precisely why the assumption list is worth memorising alongside the formula.

5. Multiple testing

A guaranteed error rate is a per-test guarantee, and almost nobody runs one test. Run one test at alpha = 0.05 and you accept a 5% false-positive rate. Run m independent tests and the probability of at least one false positive follows the complement rule from The setup in the only form that matters — the probability of at least one event is one minus the probability of none:

FWER = 1 - (1 - alpha)^m

FWER stands for family-wise error rate: the probability of making at least one false rejection anywhere in the family of tests you ran. Read the formula as “one minus the chance that all m tests individually avoid a false positive”.

Substituting alpha = 0.05 for a few values of m gives the table below. Look at the row for 20, the number of metrics a typical experiment dashboard shows.

That row is 1 - 0.95^20 = 1 - 0.358 = 0.642.

Twenty metrics on a null experiment yields a false positive about two times in three, not one time in twenty. This is not a subtle statistical concern; it is the default outcome of every experiment dashboard that reports twenty metrics. And the count is worse than you think, because it includes every segment slice, every metric variant, and every mid-flight look.

Bonferroni, and what it costs

Bonferroni correction tests each hypothesis at alpha/m instead of alpha. You have a 5% error budget; divide it evenly across the tests.

Its great virtue is that it controls FWER for any dependence structure whatsoever — the tests can be correlated in any way at all and the guarantee holds. When they are heavily correlated it is genuinely conservative, charging you as though they were independent.

The cost is exact and computable. At m = 20, alpha drops to 0.05/20 = 0.0025, so the two-sided cutoff z_{1-alpha/2} rises from 1.960 to 3.023. Since n scales with (z_a + z_b)^2 from the sample-size formula in §4, the sample-size penalty is the ratio of those squared sums:

(3.023 + 0.8416)^2 / (1.960 + 0.8416)^2  =  14.94 / 7.85  =  1.90

Correcting for 20 metrics costs you 1.9× the sample size at the same power. That is the real reason teams “forget” to correct.

Benjamini-Hochberg, and what it controls instead

Benjamini-Hochberg, usually written BH, controls a different quantity: the false discovery rate, or FDR. That is the expected fraction of your rejections that are false, rather than the probability of any false rejection at all.

The distinction is the whole point. FWER asks “did I make any mistake?” FDR asks “what share of my claimed findings are mistakes?”

The procedure has three steps:

1. Sort the p-values from smallest to largest, writing them p_(1) <= p_(2) <= ... <= p_(m). The parenthesised subscript means “the k-th smallest”, not “the k-th metric you happened to list”.
2. Find the largest rank k whose p-value still satisfies p_(k) <= (k/m)·q, where q is the FDR you are willing to accept.
3. Reject everything from rank 1 through k — including any ranks along the way that failed their own threshold.

Here it is on ten p-values with q = 0.10. The threshold column is (k/10) × 0.10, so it climbs from 0.010 at rank 1 to 0.100 at rank 10.

m = 10, q = 0.10

rank k    p_(k)      threshold (k/10)(0.10)    reject?
   1      0.001            0.010                yes
   2      0.008            0.020                yes
   3      0.020            0.030                yes
   4      0.040            0.040                yes   <- largest k that passes
   5      0.060            0.050                no
   6      0.110            0.060                no
   7      0.170            0.070                no
   8      0.240            0.080                no
   9      0.510            0.090                no
  10      0.880            0.100                no

BH rejects the top 4.        Bonferroni (0.10/10 = 0.01) rejects only the top 2.

BH is a step-up procedure. You scan for the largest k that passes, not the first k that fails.

So rank 5 failing does not end the search. You still have to check 6 through 10, and if any of them had passed, BH would reject everything up to that rank — including the ranks that failed on the way there.

Here none of them passes: 0.110 > 0.060, 0.170 > 0.070, 0.240 > 0.080, 0.510 > 0.090, 0.880 > 0.100. That is what licenses the claim “BH rejects the top 4”. You could not know it from a table that stopped at rank 6.

Note also that BH rejects rank 3 and 4 even though their p-values exceed the flat Bonferroni cutoff of 0.01, because the threshold rises with rank: a p-value only has to beat (k/m)·q, which grows as you move down the list.

Which one to use

The choice is a question about consequences, not statistics.

If one false positive is a disaster — a launch decision, a safety claim — control FWER with Bonferroni or Holm. (Holm is a step-down variant that gives the same FWER guarantee while rejecting at least as much, so it dominates plain Bonferroni.)

If you are screening thousands of candidates and will validate the survivors — feature selection, genomics, an offline metric sweep — control FDR and accept that a known fraction of your hits are noise.

The one thing that is never defensible is testing twenty metrics, reporting the one that moved, and quoting an uncorrected p-value for it.

What the corrections assume

Each correction buys its guarantee under different conditions. The right-hand column is how each one actually goes wrong.

The hardest part in practice is not the arithmetic, it is agreeing on what counts as the family. Twenty metrics, five segments and three mid-flight looks is not twenty tests; it is closer to three hundred, and no correction can rescue a family defined after the fact.

6. A/B testing

All of the theory above earns its keep the day you have to run a real experiment. The diagram below is the whole chapter applied in order, as a pipeline you walk once per experiment. Every box is a commitment you make before any data exists. The three shaded boxes are the ones that deal in sample size; the four unshaded ones are about correctness. Each box is expanded in the subsections that follow.

flowchart TD
    H["Hypothesis + ONE primary metric<br/>chosen before launch"] --> P["Power analysis<br/>n = 16·sigma^2/MDE^2 per arm"]
    P --> R["Randomize<br/>pick the unit deliberately"]
    R --> AA["A/A test<br/>validates the pipeline, not the idea"]
    AA --> RUN["Run to the pre-committed n<br/>full weeks, no peeking"]
    RUN --> V["Variance reduction<br/>CUPED on a pre-period covariate"]
    V --> AN["Analyze: effect + CI<br/>guardrails, novelty check, segments"]
    AN --> D{"Ship?"}
    LEG["shaded — the sample-size budget<br/>P sets n · RUN spends it · V buys it back"]

    style P fill:#2d6a4f,color:#fff
    style RUN fill:#2d6a4f,color:#fff
    style V fill:#2d6a4f,color:#fff
    style LEG fill:#2d6a4f,color:#fff

Sample size is the currency this pipeline spends. The power analysis computes n, running to the pre-committed n is where you spend it, and variance reduction is the only box that buys any of it back — which is why those three are shaded.

Walk the boxes in order:

1. Hypothesis and ONE primary metric, chosen before launch. A single pre-declared metric is what makes the alpha in Multiple testing mean anything. Declare twenty and you are in the 64% row of the FWER table.
2. Power analysis. Use n = 16·sigma^2/MDE^2 per arm from Hypothesis testing to learn how much traffic the question costs, before you find out the hard way.
3. Randomize, picking the unit deliberately. The randomisation unit decides what independence you are allowed to assume, and therefore what n means in every formula.
4. A/A test. Serve the same experience to both arms. This validates the pipeline, not the idea: if an experiment with no difference at all still reports significant differences, your assignment, logging or analysis is broken. You want to find that out before a real launch, not during one.
5. Run to the pre-committed n, in full weeks, without peeking. Full weeks so that weekday and weekend traffic are balanced in both arms.
6. Variance reduction. CUPED on a pre-period covariate buys back sample size for free.
7. Analyze. Report the effect with a confidence interval, check the guardrail metrics you promised not to harm, do a novelty check against exposure age, and look at segments. Only then decide whether to ship.

Randomization unit, and the variance bug it causes

The randomisation unit is the thing you flip a coin for: a user, a session, a request, or a whole geographic market. It determines what independence you are allowed to assume, and therefore what n means in every formula above.

Three rules cover almost every case:

• Randomise by user when the treatment is visible across sessions — any user-interface change. The same person must not see the old design on Monday and the new one on Tuesday.
• Randomise by session or request only when the change genuinely is invisible across visits, such as a backend ranking tweak.
• Randomise by cluster — a geography, a market, a social group — when units interfere with each other.

The design effect

The common failure mode is randomising at one level and analysing at another.

Randomise by user, analyse by page-view, and your rows are not independent. Each user contributes k correlated observations, so the sample looks k times larger than the amount of independent information it actually holds.

The design effect quantifies exactly how much you overcounted. Here is the formula and a worked instance.

effective n = n / [ 1 + (k - 1)·rho_intra ]

n = 100,000 page-view rows from 10,000 users, so k = 10 page-views per user;
rho_intra = 0.3:
  design effect = 1 + 9(0.3)     =    3.7
  effective n   = 100,000 / 3.7  = 27,027
  naive SE is understated by sqrt(3.7) = 1.92

rho_intra, read “rho intra”, is the intra-class correlation: how alike two page-views from the same user are, compared with two page-views from different users. It runs from 0 (no relationship at all) to 1 (page-views from one user are identical).

Follow the arithmetic. With k = 10 page-views per user and rho_intra = 0.3, the design effect is 1 + (10 - 1)(0.3) = 1 + 2.7 = 3.7. Divide the row count by that: 100,000 / 3.7 = 27,027. Those 100,000 rows carry about as much information as 27,000 independent ones.

Standard errors scale with 1/sqrt(n), so the naive standard error is too small by a factor of sqrt(3.7) = 1.92.

Now run that backward on a result you were about to ship. A test statistic is an effect divided by a standard error, so a standard error that is 1.92× too small makes the statistic 1.92× too big:

naive:      z = 2.17   ->  p = 0.030      "significant, ship it"
corrected:  z = 2.17 / 1.92 = 1.13  ->  p = 0.259      not remotely significant

A p = 0.03 becomes p = 0.26 from a units mismatch alone, and nothing in the analysis errors out.

Three fixes, all of which treat the user as the unit:

• Cluster-robust standard errors estimate the variance by treating each user as one block rather than each row as one observation.
• The delta method handles ratio metrics: a first-order Taylor expansion that gives the standard error of a ratio such as clicks per impression.
• A bootstrap that resamples users rather than rows, covered in Bootstrap.

Peeking, and the type-I inflation computed

The sample-size formula assumed you fixed n in advance. Peeking — checking the result repeatedly and stopping the moment it looks good — breaks that assumption, and the damage is larger than most people guess.

Under H0, the running test statistic is a random walk: a path that wanders up and down by small random steps as data accumulates, with no tendency to go anywhere in particular.

Checking that path repeatedly and stopping the first time it crosses +/-1.96 is not a 5% test. It is “the probability that a random walk ever crosses the boundary during k looks”, which is strictly larger, because you have given the walk k separate chances to wander across.

The middle column below is the real error rate. The right-hand column is what you would get if the looks were independent tests — it is there as an upper bound, to show that the truth sits between 5% and that.

That table supports two readings and you need both.

The looks are positively correlated — each one contains all the data of the previous ones — which is why 5 looks gives 14% rather than the 23% you would get from independent tests. The second column is genuinely below the third; peeking is not as bad as running independent tests.

And 14% is still nearly triple the rate you advertised.

Continuous monitoring is worse than any row in the table. By the law of the iterated logarithm — a theorem describing exactly how far a random walk wanders in its extremes — the standardized statistic has limsup = +infinity almost surely, meaning its running maximum has no ceiling. So under a true null, a dashboard checked continuously will eventually show p < 0.05 with probability 1. The only question is when.

The behavioural version is worse still, because peeking is not symmetric in practice. Teams stop early on wins and let losses “run a bit longer.”

That does two kinds of damage. It converts a 14% error rate into something much higher, and it biases the measured effect size upward, because you stopped precisely on a favourable random fluctuation. Even the tests you would have won report inflated lifts.

The table above is not folklore, and running the simulation is the fastest way to convince a sceptical stakeholder. The function below generates data with no effect whatsoever — every observation is drawn from a Normal with mean exactly 0 — then applies the stopping rule and counts how often it declares a winner. Watch the break: that is the peeking.

import math
import random


def peeking_type_i(n_looks, n_total=4000, trials=40_000, z_crit=1.96, seed=1):
    """False-positive rate under a TRUE NULL when you stop at the first
    look that crosses z_crit. n_looks=1 recovers the nominal alpha."""
    rng = random.Random(seed)
    step, hits = n_total // n_looks, 0
    for _ in range(trials):
        running_sum = n = 0.0
        for _ in range(n_looks):
            for _ in range(step):
                running_sum += rng.gauss(0.0, 1.0)      # H0: mean is exactly 0
                n += 1
            if abs(running_sum / math.sqrt(n)) > z_crit:
                hits += 1
                break                                   # stop on first "win"
    return hits / trials


# Actual returns at the defaults above (seed=1, trials=40,000):
#   1 look -> 0.0498   2 -> 0.0806   3 -> 0.1068   5 -> 0.1410   10 -> 0.1900

The comment on the last two lines records what the function actually returns at the seed and trial count written into the signature. Paste the block and you can check it against a fixed target rather than against a number that moves every run.

Those five results reproduce the table above to within simulation noise. At 40,000 trials that noise is about +/- 0.003 at two standard errors on a rate near 0.10, so do not read a disagreement in the third decimal place as a bug.

Sequential testing

If you want to look, use a method that priced the looking in. Three families exist, and they all trade the same currency: validity at every moment is bought with sample size. The right-hand column is that price.

The middle row carries the most jargon, so unpack it.

mSPRT stands for mixture sequential probability ratio test. The likelihood ratio it accumulates is the same one behind the likelihood-ratio test of Confidence intervals and the interpretation nobody gets right.

A martingale is a running quantity whose expected future value equals its current value — a fair game, in the gambling sense. That property is exactly what stops repeated looking from accumulating error: because the statistic is a martingale under the null, a threshold crossing stays rare no matter how often you look. That is what “always valid” means.

Confidence sequences are the interval version of the same idea: a sequence of intervals that simultaneously cover the truth at every sample size, so you may read them at any moment without penalty.

There is no free lunch: continuous validity is bought with power, and the price is a longer test for the same detectable effect. The right framing in an interview is that this is a trade you make deliberately — a platform serving hundreds of experiments should pay 15% more traffic to remove peeking as an entire class of bug.

Novelty, primacy, and network effects

Two behavioural effects break the assumption that a unit’s outcome is a stable draw from a fixed distribution, and they break it in ways no statistical correction can repair.

Novelty and primacy are the same phenomenon with opposite signs: users react to change, not to the new state.

A redesign gets clicked because it is different. That is novelty, and the measured effect decays toward zero as the novelty wears off. Or the redesign gets ignored because it is unfamiliar. That is primacy, and the effect grows as people learn the new layout.

Diagnose either by plotting the treatment effect against days-since-first-exposure per cohort, not against calendar day. This is the “novelty check” box in the pipeline diagram. Calendar day mixes users who joined the experiment on different days, which smears the pattern away; exposure age lines them up.

A curve that decays toward zero over two weeks is a novelty effect, and shipping on the week-1 number is a mistake.

Three cures: run longer; restrict the analysis to users exposed for at least k days; keep a small long-term holdout that never receives the treatment.

Network effects and interference break a different assumption — that a unit’s outcome depends only on its own assignment.

Marketplaces are the canonical case. Showing treatment users better inventory takes that inventory away from control. Control degrades, and the measured lift is doubly overstated: treatment gained what control lost, so the gap between them counts the same transfer twice. Social products have the same structure through the feed.

The fixes here are structural, not statistical:

• Cluster randomization by geography or social component, so that interacting units land in the same arm.
• Switchback designs, which flip the whole market between treatment and control over time windows, so every unit experiences both.
• Budget-split designs, which give each arm its own supply.

CUPED, derived

CUPED stands for Controlled-experiment Using Pre-Experiment Data. It is the highest-return trick in experimentation and it is four lines of algebra.

The idea in one sentence: a large part of the variation in any user metric is just the user being who they already were, and if you measured that part before the experiment you can subtract it out without touching the treatment effect.

Concretely: heavy spenders spend a lot in both arms. That variation is noise as far as your test is concerned, and it is noise you can predict from last month’s data.

Let X be a covariate — any variable you measured on the same units — recorded before randomization. Usually it is the same metric in the pre-period. Let Y be the outcome metric you care about.

Note that this section’s theta is a tuning constant, not the generic unknown parameter of §1. The name is standard in the CUPED literature.

The four lines below define the adjusted metric, show that adjusting cannot bias it, write down its variance, and then pick the theta that makes that variance as small as possible.

Y_cuped = Y - theta·( X - E[X] )

E[Y_cuped] = E[Y]                                 unbiased for ANY theta, since E[X - E[X]] = 0

Var(Y_cuped) = Var(Y) + theta^2·Var(X) - 2·theta·Cov(X, Y)

d/dtheta = 2·theta·Var(X) - 2·Cov(X, Y) = 0   ->   theta* = Cov(X, Y) / Var(X)

Var(Y_cuped) at theta* = Var(Y) - Cov(X,Y)^2/Var(X) = Var(Y)·( 1 - rho^2 )

Line by line:

The adjustment. Y_cuped subtracts theta times each unit’s deviation from the average covariate. A user whose pre-period value was above average gets a deduction; a user below average gets a credit.

Why it is safe. The mean is untouched for any theta, because X - E[X] averages to zero by construction, and theta · 0 = 0. So no choice of theta can bias your effect estimate. That is what licenses tuning theta freely.

The variance. This is Var(A - B) = Var(A) + Var(B) - 2·Cov(A, B) with A = Y and B = theta·(X - E[X]). The constant theta comes out squared in the variance term and linearly in the covariance term. Cov(X, Y), read “covariance of X and Y”, measures how much the two move together.

The minimisation. Differentiate that variance with respect to theta, set to zero, and you get theta* = Cov(X, Y) / Var(X).

The payoff. Substituting theta* back gives Var(Y)·(1 - rho^2), where rho, read “rho”, is the correlation between X and Y — a number between −1 and 1 measuring how strongly they track each other.

The variance falls by exactly rho^2, the squared correlation between the metric and its pre-period value.

theta* is also the OLS slope of Y on X, which is why CUPED and regression adjustment (ANCOVA) are the same estimator wearing different names. OLS is ordinary least squares, the standard way of fitting a straight line by minimising squared errors. ANCOVA is analysis of covariance, the classical name for a comparison of groups that first adjusts for a covariate.

The table turns rho into the two numbers you care about: how much extra sample the adjustment is worth, and how small an effect you can now detect.

Read the bolded row to fix the pattern. At rho = 0.6, the variance retained is 1 - 0.6^2 = 1 - 0.36 = 0.64. Since n scales with variance, the effective sample multiplier is 1/0.64 = 1.56×. And the MDE scales with the square root of the variance, so it shrinks to sqrt(0.64) = 0.80, or 80% of what it was.

In calendar terms: at rho = 0.6 a four-week test becomes a 4 × 0.64 = 2.6-week test, for the cost of one join against last month’s data.

The one non-negotiable condition is that X must be measured strictly before randomization. Adjust on anything the treatment could have influenced and you are conditioning on a post-treatment variable — the mediator trap of Simpsons paradox — and the estimate becomes biased rather than merely noisy.

That is also why the pre-period value of the same metric is the standard choice. It is the highest-rho covariate available and it is unambiguously pre-treatment.

A second, quieter condition: theta* is estimated from the same data it is applied to. On a large experiment the resulting bias is negligible; on a few hundred units it is not. The honest fix is to estimate theta on units held out from the analysis, or by cross-fitting — split the data in half, estimate theta on one half, apply it to the other, then swap.

7. Bootstrap

Everything so far leaned on a formula for the sampling distribution. When no formula exists, there is a general-purpose escape hatch: simulate the sampling distribution instead of deriving it.

The procedure is one sentence. Resample n rows from your data with replacement — meaning the same row can be drawn more than once — repeat that B times to get B synthetic datasets, recompute the statistic on each, and use the spread of those B values as the sampling distribution.

The logic is that the sample is your best available stand-in for the population, so drawing from the sample imitates drawing from the population.

The code below is that procedure with nothing added. The line doing the real work is x[rng.integers(0, n, n)], which draws n random indices with replacement and takes those rows.

import numpy as np


def bootstrap_ci(x, stat=np.median, B=10_000, alpha=0.05, seed=0):
    """Percentile bootstrap CI for any statistic."""
    rng = np.random.default_rng(seed)
    x = np.asarray(x)
    n = len(x)
    reps = np.empty(B)
    for b in range(B):
        reps[b] = stat(x[rng.integers(0, n, n)])
    lo, hi = np.quantile(reps, [alpha / 2, 1 - alpha / 2])
    return float(reps.std(ddof=1)), (float(lo), float(hi))

The function returns two things.

First, the standard deviation of the B recomputed statistics — which is the standard error, by definition, since the standard error is the spread of the estimate across repeated samples.

Second, the 2.5th and 97.5th percentile of those B values, a percentile being the value below which that share of the numbers fall. That pair is the interval, which is why this is called the percentile bootstrap.

Checking it against the closed form

Work it through on x = [1, 2, 3, 4, 5, 6, 7, 8, 100] — nine points with one heavy outlier. Compare the bootstrap against the textbook formulas where a textbook formula exists.

                     closed form                 bootstrap (B -> infinity)
mean   = 15.11       SE = s/sqrt(n) = 10.64      SE = 10.03
median =  5.00       SE = 1 / (2·f(m)·sqrt(n))   SE =  3.87
                     requires a density estimate

The bootstrap column is quoted at the B -> infinity limit rather than at some finite B, because a finite B with no stated seed is not a number anyone can check. Call bootstrap_ci(x, stat=np.mean, B=20_000, seed=s) from the code above and you get 9.98, 10.10 and 10.04 for s = 0, 1, 2. That spread is the simulation talking, not the statistic.

Both limits are exact numbers, not approximations:

• The mean’s limit is s/sqrt(n) × sqrt((n-1)/n) = 10.0301.
• The median’s limit is 3.8682. It is available in closed form because with n = 9 the resampled median can only ever be one of the nine observed values, so its distribution is a finite sum of binomial terms rather than something you have to simulate.

Two things to read off that comparison.

For the mean, the bootstrap reproduces the closed form. The entire gap is a divisor: the bootstrap’s internal variance uses 1/n where the sample standard deviation uses 1/(n-1). Check it — 10.64 × sqrt(8/9) = 10.03, exactly the bootstrap figure. That agreement is the validation that the method is not magic.

For the median, there is no usable closed form. The asymptotic standard error 1/(2·f(m)·sqrt(n)) depends on f(m), the probability density at the median — how tightly the data pile up right at the middle value. You would have to estimate f(m) yourself, which means choosing a smoothing width (how wide a window of neighbouring points to average over) and carrying its assumptions. The bootstrap needs none of that.

The bootstrap replaces “derive the sampling distribution analytically” with “simulate it,” and the only cost is compute.

When to reach for it

The bootstrap beats a closed form in three situations, and all three are common:

1. The statistic has no clean standard error. Medians and other percentiles qualify. So do ratios of metrics, such as revenue per session. So does a difference in AUC — the area under the receiver-operating-characteristic curve, a single number summarising how well a classifier ranks positives above negatives. So does a Gini coefficient, which measures how unevenly a total is spread across a population. And so does anything you computed with a pipeline rather than a formula.

2. The metric is skewed and the central limit theorem has not yet arrived. The central limit theorem says the average of many independent draws becomes Gaussian. How many “many” has to be depends on the skew — that is the Berry-Esseen bound from The distributions and what each one models, which puts an explicit rate on how slowly a skewed distribution converges.

3. You can compute the point estimate but nobody has done the delta-method algebra for its standard error.

When it silently fails

It fails in four situations, and every one of them gets asked in interviews.

Every row is the same three assumptions failing in a different way. The bootstrap requires that rows be independent and identically distributed, that the statistic vary smoothly with the data, and that the sample be large enough to resemble the population.

Two pieces of vocabulary in that table are worth pinning down.

The empirical distribution is your sample treated as if it were the whole population: every observed row equally likely, nothing else possible. That is exactly what resampling draws from, and it is why a small sample poisons the method — the bootstrap can only ever tell you about the population it was shown.

Asymptotics means the guarantees that hold only as n grows without bound. A parameter pinned at the edge of its range — a variance, which cannot go below zero — breaks them, because the estimate cannot wander symmetrically in both directions.

The two repairs named in the table follow the same logic as each other. A block bootstrap resamples contiguous stretches rather than single rows, preserving the time structure. A cluster bootstrap resamples whole groups, preserving the within-group correlation.

The compute budget

B = 1,000 is enough for a standard error. Percentile intervals need B = 10,000, because there you are estimating the tail quantiles — the cut points far out in the distribution, a quantile being the same idea as a percentile on a 0-to-1 scale — and the tails are where simulation noise is worst.

For A/B tests, resample the randomization unit. Bootstrap users, not events, or you have reintroduced the design-effect bug from Randomization unit and the variance bug it causes inside a method you adopted to avoid it.

8. Simpson’s paradox

All the machinery so far assumes you are running the right comparison in the first place. One result shows why no amount of statistics can substitute for knowing what causes what: a treatment can win in every subgroup and lose overall, and which number is correct is not a question the data can answer.

Start with the classic table: the real kidney-stone treatment comparison of Charig, Webb, Payne and Wickham, British Medical Journal 292:879-882 (1986). Each cell shows a success rate with the raw counts behind it. Read across each row, then compare the bolded numbers.

A wins in both subgroups and loses overall.

That is not an arithmetic error, and it is worth checking rather than believing. For A: (81 + 192)/350 = 273/350 = 0.780. For B: (234 + 55)/350 = 289/350 = 0.826.

The mechanism is the case mix.

Large stones are harder to treat regardless of method — 73% and 69% success, against 93% and 87% for small ones. And the two treatments did not get the same mix of cases. A was given overwhelmingly to large stones (263 of its 350 cases); B was given overwhelmingly to small ones (270 of its 350).

The overall rates are weighted averages with different weights, so each treatment’s total is dragged toward the difficulty of the cases it happened to receive.

Give both treatments the same case mix — weight the two subgroups equally — and the reversal disappears. A gets (0.93 + 0.73)/2 = 0.83. B gets (0.87 + 0.69)/2 = 0.78. A wins.

So which number is right, the subgroup one or the overall one? That is a causal question, and the answer depends entirely on the role of the third variable.

The diagram below draws the three roles a third variable can play. An arrow means “causes”. The green node is the one case where you must adjust; the two red nodes are cases where adjusting makes things worse.

flowchart LR
    subgraph CONF["CONFOUNDER -- you must adjust"]
        S1(("stone<br/>size")) --> T1(("treatment"))
        S1 --> O1(("outcome"))
        T1 --> O1
    end
    subgraph MED["MEDIATOR -- never adjust"]
        T2(("drug")) --> M2(("blood<br/>pressure"))
        M2 --> O2(("heart<br/>attack"))
    end
    subgraph COL["COLLIDER -- never adjust"]
        A3(("talent")) --> H3(("hired"))
        B3(("interview<br/>prep")) --> H3
    end

    style S1 fill:#2d6a4f,color:#fff
    style M2 fill:#9d0208,color:#fff
    style H3 fill:#9d0208,color:#fff

Taking the three panels in order:

• Confounder — adjust. Stone size causes both the treatment choice and the outcome, so it is a common cause of the two things you are comparing. You must adjust, because otherwise the difference you measure is partly a difference in case mix. In the kidney-stone table the subgroup numbers are the correct ones, and treatment A is better.
• Mediator — never adjust. A drug lowers blood pressure, which in turn lowers heart attacks, so blood pressure sits on the causal path between the treatment and the outcome. Conditioning on it blocks part of the very effect you are trying to measure. Here the aggregate number is the correct one. This is exactly why CUPED requires a strictly pre-randomization covariate: a post-treatment covariate is a mediator.
• Collider — never adjust. Talent and interview prep both cause getting hired, so “hired” is a common effect rather than a common cause. Conditioning on it creates a spurious association out of nothing. Among people who were hired, weak preparation implies high talent, so the two look negatively related even though they are independent in the population at large. This is the same selection mechanism as Conditional probability and independence.

No statistical test tells you which diagram you are in. You need domain knowledge about what causes what, which is why “should I control for this variable?” is never answerable from the data alone.

The diagram you believe is called a DAG — a directed acyclic graph, meaning arrows that never form a loop. It is an assumption you assert, not a result you compute.

Where this bites in machine-learning work

Three concrete versions of the same trap:

• A model that improves overall AUC while degrading it in every segment, or the reverse, because the segments have different base rates — different underlying frequencies of the positive outcome before any model touches them — and different volumes.
• A launch that looks positive in aggregate because the treatment shifted the traffic composition rather than improving anything.
• Any before/after comparison across a period when the user mix changed.

The prophylactic is cheap: always report the primary metric sliced by the two or three dimensions you know affect both assignment and outcome, and reconcile the slices against the total before anyone sees the number.

9. The assumption ledger

What each procedure above requires of your data, and the specific way production experiments tend to break it, belongs in one place: the previous eight sections indexed by failure rather than by topic, because failure is the order in which problems arrive.

Use it as a checklist. Find the row for the procedure you are about to run, and check the third column against what you know about your data.

The pattern across every row is the same: the violation does not raise an exception, it changes the number. That is the whole argument for knowing the assumption list by heart rather than looking it up when something goes wrong, because nothing will tell you that something has gone wrong.

Cheat sheet

One line per idea in the chapter. If you can reconstruct the derivation behind a row from the middle column alone, you know that row.

Next: 03 — Linear Algebra — the other half of the mathematics, and the one that decides what runs fast.