← Back to the Applied Science guide

Statistical Inference

Data is a noisy sample of the world. Inference is the craft of saying how much the sample tells you, and how sure you can be.

Every launch call, every model comparison, and every dashboard alert rests on inference. The interview Q&A page drills the common questions. This page builds the theory under them. Learn it once and you can derive the answers you forget.

Each section gives a plain definition, the math, the intuition, and where it bites in real work. Each ends with short interview checks. The code runs in plain Python 3 with no extra packages.

1. Estimators and their properties

Write the unknown parameter as θ and the estimator as θ̂. Because the sample is random, θ̂ is a random variable. It has its own distribution, called the sampling distribution. Every property below describes that distribution.

Bias, variance and MSE

• Bias. Bias(θ̂) = E[θ̂] − θ. It measures how far the average guess sits from the truth.
• Variance. Var(θ̂) = E[(θ̂ − E[θ̂])²]. It measures how much the guess jumps between samples. Its square root is the standard error.
• Mean squared error. MSE(θ̂) = E[(θ̂ − θ)²]. It measures total error.

The key identity follows from adding and subtracting E[θ̂] inside the square. The cross term has mean zero.

MSE(θ̂) = E[(θ̂ − Eθ̂ + Eθ̂ − θ)²]
       = E[(θ̂ − Eθ̂)²] + (Eθ̂ − θ)² + 2(Eθ̂ − θ)·E[θ̂ − Eθ̂]
       = Var(θ̂) + Bias(θ̂)² + 0

Worked case. The variance estimator that divides by n has E[σ̂²] = ((n − 1)/n)σ². It is biased low. Dividing by n − 1 removes the bias. That is Bessel's correction. But the biased version has lower variance, and for normal data it has lower MSE too. Dividing by n + 1 gives the lowest MSE of all. Unbiased is not the same as best.

Consistency

An estimator is consistent if θ̂_n → θ in probability as n grows. A handy sufficient condition is that bias and variance both go to zero. Then MSE goes to zero, and Chebyshev gives convergence. Consistency is the bare minimum. An estimator that stays wrong with infinite data is useless.

Bias and consistency are separate ideas. The divide-by-n variance is biased but consistent. The estimator "use only the first point" is unbiased for the mean but not consistent.

Efficiency, Fisher information and Cramér-Rao

Let f(x; θ) be the density and ℓ(θ) = log f. The score is the slope ∂ℓ/∂θ. Under regularity conditions the score has mean zero. Its variance is the Fisher information for one point.

I(θ) = E[ (∂ log f(X; θ) / ∂θ)² ] = −E[ ∂² log f(X; θ) / ∂θ² ]

For n iid points:  In(θ) = n · I(θ)

The Cramér-Rao lower bound says no unbiased estimator can beat the inverse information.

Var(θ̂) ≥ 1 / (n · I(θ))                   for unbiased θ̂
Var(θ̂) ≥ (1 + b'(θ))² / (n · I(θ))     if Bias = b(θ)

The proof is one line of Cauchy-Schwarz. Cov(θ̂, score) = 1 for an unbiased estimator. So 1 ≤ Var(θ̂) · n I(θ).

• Bernoulli(p). ℓ = x log p + (1 − x) log(1 − p). I(p) = 1 / (p(1 − p)). The bound is p(1 − p)/n. The sample proportion hits it exactly.
• Normal mean, known σ. I(μ) = 1/σ². The bound is σ²/n. The sample mean hits it.

An unbiased estimator that reaches the bound is efficient. The ratio of the bound to the actual variance is its efficiency. The sample median of normal data has asymptotic efficiency 2/π ≈ 0.64. You need about 57% more data to match the mean. But the median wins on heavy tails.

Why it matters in practice

• Shrinkage beats unbiased. Ridge, empirical Bayes, and CUPED all add a little bias to cut variance. MSE is the target.
• Experiment sizing. Standard errors scale as 1/√n. Halving the error needs four times the traffic.
• Metric choice. A capped or trimmed revenue metric is biased but far less noisy. It often detects real effects sooner.
• Design of experiments. Fisher information guides which features to log and which arms to sample in adaptive tests.

Interview check

• Is an unbiased estimator always better? No. Compare by MSE. A biased estimator with much lower variance often wins.
• What does Fisher information measure? The expected curvature of the log-likelihood. It sets the best possible variance through Cramér-Rao.
• Biased but consistent example? The divide-by-n variance. Its bias is −σ²/n, which vanishes as n grows.

2. Method of moments vs MLE

Method of moments

Match as many moments as you have parameters. Take a Gamma with shape k and scale θ. Its mean is kθ and its variance is kθ². Solve with the sample mean x̄ and variance s².

k̂ = x̄² / s²        θ̂ = s² / x̄

MoM is fast and needs no optimizer. It is consistent under mild conditions. But it ignores much of the shape of the data, so it is often inefficient. It also gives good starting points for MLE.

Maximum likelihood

Maximize ℓ(θ) = ∑_i log f(x_i; θ). For Bernoulli data with k successes in n trials:

ℓ(p) = k log p + (n − k) log(1 − p)
dℓ/dp = k/p − (n − k)/(1 − p) = 0   ⇒   p̂ = k/n

Asymptotic normality of the MLE

Under regularity conditions the MLE is consistent and asymptotically normal. Its variance reaches the Cramér-Rao bound in the limit.

√n (θ̂ − θ) → N(0, I(θ)−1)

Sketch. Expand the score around the truth. 0 = ℓ'(θ̂) ≈ ℓ'(θ) + (θ̂ − θ) ℓ''(θ). So θ̂ − θ ≈ −ℓ'(θ) / ℓ''(θ). The numerator is a sum of iid scores. By the CLT it is about N(0, n I). The denominator is about −n I by the law of large numbers. Divide and the result follows.

Three more facts make MLE the default.

• Invariance. If θ̂ is the MLE of θ, then g(θ̂) is the MLE of g(θ).
• Standard errors for free. Invert the observed information, the negative Hessian at θ̂. That gives the covariance matrix.
• Misspecification. If the model is wrong, the MLE converges to the KL-closest model. Its variance is then the sandwich form A⁻¹BA⁻¹, not I⁻¹.

Why it matters in practice

• Log loss is MLE. Training a classifier with cross-entropy is maximum likelihood for a Bernoulli or categorical model.
• Squared loss is Gaussian MLE. MSE regression assumes constant-variance normal noise. Heavy tails argue for Huber or quantile loss.
• Uncertainty from the Hessian. Laplace approximations in Bayesian deep learning reuse this exact idea.
• Quick fits. MoM works well for streaming dashboards where you only store running sums.

Interview check

• Why prefer MLE to MoM? MLE is asymptotically efficient and invariant. MoM is simpler but often has higher variance.
• What is the asymptotic variance of the MLE? The inverse Fisher information, 1/(n I(θ)), under regularity conditions.
• Give a case where MLE is biased. The normal variance MLE divides by n. The Uniform(0, θ) MLE is the sample max.

3. Hypothesis testing

The framework

• Type I error. Reject a true H₀. Its probability is the size α.
• Type II error. Keep a false H₀. Its probability is β. Power is 1 − β.
• p-value. The chance, under H₀, of a statistic at least as extreme as the one seen. It is not the chance that H₀ is true.
• Duality. A level-α test rejects θ₀ exactly when a 1 − α confidence interval excludes θ₀.

Neyman-Pearson lemma

Take two simple hypotheses, H₀: θ = θ₀ and H₁: θ = θ₁. Among all tests of size α, the most powerful one rejects when the likelihood ratio is large.

Λ(x) = L(θ1; x) / L(θ0; x) > k,     with k set so P0(Λ > k) = α

For normal data the ratio is monotone in x̄. So "reject when x̄ is large" is most powerful. Because that test does not depend on θ₁, it is uniformly most powerful for all θ₁ > θ₀. Two-sided UMP tests do not exist in general.

Generalized likelihood ratio tests

Real hypotheses are composite. Maximize the likelihood under each and compare.

λ = supθ ∈ Θ0 L(θ) / supθ ∈ Θ L(θ)

Wilks:  −2 log λ → χ²d,   d = (free params in Θ) − (free params in Θ0)

The LR test has two cousins. The Wald test uses only the full fit, (θ̂ − θ₀)/SE. The score test uses only the null fit, through the slope at θ₀. All three agree in large samples. In small samples LR is usually the most reliable. Wald can behave badly near a boundary.

Which test when

• z-test. Variance known, or n large enough for the CLT. The standard tool for A/B tests on means and proportions with large traffic. Statistic: (x̄ − μ₀)/(σ/√n).
• t-test. Variance unknown and estimated. Exact for normal data, and robust for moderate n. Use Welch's version by default. It does not assume equal variances. Use the paired t-test for before and after on the same units.
• Chi-square test. Counts in categories. Goodness of fit, or independence in a contingency grid. Statistic: ∑ (O − E)²/E. Needs expected counts of about 5 or more per cell. Otherwise use Fisher's exact test. A 2 × 2 chi-square equals a two-proportion z-test squared.
• Mann-Whitney U. Two independent samples, no normality needed. It ranks the pooled data. It tests whether P(X > Y) = 1/2. It is a test of medians only if the two shapes match. Use it for skewed, ordinal, or outlier-heavy data.
• Wilcoxon signed-rank. The paired version of Mann-Whitney.

Why it matters in practice

• Revenue metrics. Heavy tails make t-tests slow to converge. Winsorize, or use Mann-Whitney, or bootstrap. Know which question each answers.
• Sample size. For a two-sample z-test, n per arm ≈ 2(z_1−α/2 + z_1−β)²σ²/δ². That is about 16σ²/δ² for 80% power at α = 0.05.
• Model comparison. Nested models compare with an LR test. Likelihood ratio is also the logic behind sequential tests like SPRT.
• Sample ratio mismatch. A chi-square test on arm counts catches broken randomization before you read any metric.

Interview check

• State the Neyman-Pearson lemma. For simple vs simple, the likelihood ratio test is most powerful at its size.
• t-test or Mann-Whitney for session length? It depends on the question. Mean impact calls for t or bootstrap. A shift in the typical user calls for Mann-Whitney.
• What does Wilks' theorem give you? −2 log λ is about chi-square, with df equal to the number of restricted parameters.

4. Confidence intervals

The randomness lives in the interval, not in θ. Once you compute [0.12, 0.18], it either covers the truth or it does not. "95% chance θ is in this interval" is a Bayesian statement. A credible interval makes that claim. A confidence interval does not.

Wald intervals

Take the estimate and add plus or minus z times its standard error. For a proportion:

p̂ ± z1−α/2 √( p̂(1 − p̂) / n )

It is simple and fine for large n with p near 0.5. It fails near 0 or 1. At p̂ = 0 it gives the absurd interval [0, 0]. Coverage can fall far below the nominal level. The coverage also swings with n because the data are discrete.

Wilson score interval

Invert the score test instead. Solve |p̂ − p| ≤ z√(p(1 − p)/n) for p. The standard error uses the hypothesized p, not p̂.

        p̂ + z²/(2n)  ±  z √( p̂(1 − p̂)/n + z²/(4n²) )
CI  =  ──────────────────────────────────────
                         1 + z²/n

The center shrinks toward 1/2. The interval never leaves [0, 1]. Coverage stays close to nominal even for small n and rare events. Agresti-Coull is a close cousin: add z²/2 successes and failures, then use Wald.

import math, random

def wald(k, n, z=1.96):
    p = k / n
    h = z * math.sqrt(p * (1 - p) / n)
    return p - h, p + h

def wilson(k, n, z=1.96):
    p = k / n
    d = 1 + z * z / n
    c = (p + z * z / (2 * n)) / d
    h = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d
    return c - h, c + h

def coverage(ci, p, n, reps=20000, seed=0):
    rng = random.Random(seed)
    hit = 0
    for _ in range(reps):
        k = sum(rng.random() < p for _ in range(n))
        lo, hi = ci(k, n)
        hit += lo <= p <= hi
    return hit / reps

for p in (0.5, 0.05, 0.01):
    print(f"p={p:<5} n=50  Wald={coverage(wald, p, 50):.3f}  "
          f"Wilson={coverage(wilson, p, 50):.3f}")

p=0.5   n=50  Wald=0.933  Wilson=0.933
p=0.05  n=50  Wald=0.920  Wilson=0.961
p=0.01  n=50  Wald=0.398  Wilson=0.911

At p = 0.01 the Wald interval covers only 40% of the time. Most samples have zero successes, and then Wald has zero width. Wilson stays near 95%.

Bootstrap intervals

Resample the data with replacement B times. Recompute the statistic each time. The spread of those replicates stands in for the sampling distribution.

• Percentile. Take the α/2 and 1 − α/2 quantiles of the replicates. Simple, and it respects the range of the statistic. It is transformation invariant. But it assumes the replicate distribution is centered right and has the right shape. With bias or skew, it under-covers.
• Basic. Reflect the percentiles around θ̂: [2θ̂ − q_hi, 2θ̂ − q_lo]. Handles bias but not skew.
• BCa (bias-corrected and accelerated). Adjusts which quantiles you take. Two constants do the work.

z0 = Φ−1( share of replicates below θ̂ )                       bias correction
a  = ∑(θ̄(·) − θ(i))³ / ( 6 [∑(θ̄(·) − θ(i))²]3/2 )      acceleration, from the jackknife

α1 = Φ( z0 + (z0 + zα/2) / (1 − a(z0 + zα/2)) )      lower quantile to use
α2 = Φ( z0 + (z0 + z1−α/2) / (1 − a(z0 + z1−α/2)) )  upper quantile to use

z₀ fixes median bias. The acceleration a fixes a standard error that changes with θ. BCa is second-order accurate. Its coverage error shrinks like 1/n, against 1/√n for percentile. The cost is a jackknife pass and more replicates, often 2,000 or more.

Why it matters in practice

• Rare events. Click, fraud, and crash rates sit near zero. Report Wilson or exact intervals, never Wald.
• Odd statistics. Medians, percentiles, AUC, and NDCG have no simple SE. Use a bootstrap, with BCa when skewed.
• Offline eval. Put a bootstrap CI on every model metric. A 0.2% AUC gain inside the noise is not a win.
• Clustered data. Resample users, not events. Event-level resampling gives intervals that are far too narrow.

Interview check

• Interpret a 95% CI. The method covers the truth in 95% of repeated samples. It is not a 95% probability for this interval.
• Why Wilson over Wald? Wald collapses near 0 and 1 and under-covers. Wilson inverts the score test and stays close to nominal.
• When does the percentile bootstrap fail? With biased or skewed statistics, and small samples. BCa corrects both effects.

5. The delta method

The math

Suppose √n(θ̂ − θ) → N(0, σ²) and g is differentiable with g'(θ) ≠ 0. A first-order Taylor expansion gives g(θ̂) ≈ g(θ) + g'(θ)(θ̂ − θ). So:

√n ( g(θ̂) − g(θ) ) → N( 0, g'(θ)² σ² )

Multivariate:  Var( g(θ̂) ) ≈ ∇g(θ)T Σ ∇g(θ)

Example. The log odds of a proportion. g(p) = log(p/(1 − p)) and g'(p) = 1/(p(1 − p)). So Var(logit p̂) ≈ 1/(n p(1 − p)).

Ratio metrics

Many product metrics are ratios of two sums. Click-through rate is clicks over impressions. Revenue per session is revenue over sessions. The randomization unit is the user, but the denominator counts something else. Treating each session as independent ignores the correlation within a user. The naive SE is then too small.

Let X_i and Y_i be user i's numerator and denominator. The metric is R = X̄ / Ȳ. Apply the delta method with g(x, y) = x/y. The gradient is (1/μ_Y, −μ_X/μ_Y²).

Var(R) ≈ (1/n) · [ σX²/μY²  −  2 μX σXY/μY³  +  μX² σY²/μY⁴ ]

Plug in sample moments. The demo checks the result against a user-level bootstrap.

import random, statistics as st

# Delta method vs bootstrap for a ratio metric: clicks per session.
rng = random.Random(1)
users = []
for _ in range(2000):
    s = 1 + int(rng.expovariate(1 / 4))          # sessions per user
    c = sum(rng.random() < 0.3 for _ in range(s)) # clicks per user
    users.append((c, s))

X = [c for c, _ in users]; Y = [s for _, s in users]
n = len(users); mx, my = st.mean(X), st.mean(Y)
r = mx / my
vx, vy = st.variance(X), st.variance(Y)
cxy = sum((x - mx) * (y - my) for x, y in users) / (n - 1)
var_r = (vx / my**2 - 2 * mx * cxy / my**3 + mx**2 * vy / my**4) / n
print(f"ratio = {r:.4f}, delta-method SE = {var_r ** 0.5:.5f}")

boots = []
for _ in range(2000):
    smp = [users[rng.randrange(n)] for _ in range(n)]
    boots.append(sum(c for c, _ in smp) / sum(s for _, s in smp))
print(f"bootstrap SE         = {st.stdev(boots):.5f}")

ratio = 0.3018, delta-method SE = 0.00504
bootstrap SE         = 0.00495

The two agree to within 2%. The delta method costs one pass over the data. The bootstrap costs 2,000 passes.

Why it matters in practice

• A/B platforms. Most large experiment platforms use the delta method for every ratio metric. It is cheap and works on aggregated sums.
• Relative lift. Lift is (μ_T − μ_C)/μ_C. Its CI needs the delta method, or Fieller's method when μ_C is noisy.
• Log transforms. The SE of log(x̄) is about SE(x̄)/x̄, the coefficient of variation.
• Failure mode. When the denominator mean is near zero, or n is small, linearization breaks. Fall back to the bootstrap.

Interview check

• Why not compute CTR variance from per-impression Bernoullis? Impressions from one user are correlated. The user is the unit, so the naive SE is too small.
• State the delta method. Var(g(θ̂)) ≈ g'(θ)² Var(θ̂), from a first-order Taylor expansion.
• When does it fail? When g'(θ) = 0, or g is far from linear at the noise scale. Also when the denominator is near zero.

6. Linear regression theory

Derivation

X is n × p with full column rank. Minimize the residual sum of squares.

RSS(β) = (y − Xβ)T(y − Xβ)
∂RSS/∂β = −2 XT(y − Xβ) = 0      (the normal equations)
β̂ = (XTX)−1 XTy

Var(β̂ | X) = σ² (XTX)−1        σ̂² = RSS / (n − p)

The normal equations say X^Te = 0. The residuals are orthogonal to every column of X. So OLS is a projection. The hat matrix H = X(X^TX)⁻¹X^T projects y onto the column space of X. Its diagonal entries are the leverages.

Gauss-Markov

Under these assumptions, OLS is the best linear unbiased estimator (BLUE). No other linear unbiased estimator has smaller variance for any coefficient combination.

1. Linearity. y = Xβ + ε in the parameters. Features can be nonlinear transforms.
2. Full rank. No exact linear dependence among columns.
3. Exogeneity. E[ε | X] = 0. This gives unbiasedness. It fails with omitted confounders or reverse causality.
4. Spherical errors. Var(ε | X) = σ²I. Constant variance and no correlation between errors.

Normality is not on the list. It only buys exact t and F tests in small samples. With large n the CLT makes the usual tests valid anyway. Note that "best" is only among linear unbiased estimators. Ridge is biased and can have lower MSE.

Heteroskedasticity-robust standard errors

If the error variance changes across rows, β̂ stays unbiased. But σ²(X^TX)⁻¹ is now wrong. The fix keeps β̂ and changes the variance to a sandwich.

HC0:  V̂ = (XTX)−1 ( ∑i ei² xi xiT ) (XTX)−1
HC1:  HC0 · n/(n − p)          HC3:  weight each ei² by 1/(1 − hii)²

HC3 behaves best in small samples. For grouped data, sum the score within clusters first. That gives cluster-robust SEs, which you need when rows from one user are correlated.

Multicollinearity

When columns are nearly collinear, X^TX is near singular. Coefficients stay unbiased, but their variances blow up. The variance inflation factor measures this.

VIFj = 1 / (1 − Rj²)        Rj² from regressing xj on the other features

A VIF over 10 means the SE of β_j is over 3 times what it would be alone. Predictions can still be fine. Only the split of credit between correlated features is unstable. Fixes include dropping or combining features, ridge, or simply not interpreting those coefficients.

Interpreting coefficients

• Level-level. One unit more of x_j, others held fixed, goes with β_j units more of y.
• Log-level. With log y, one unit of x_j goes with a 100(e^β_j − 1)% change in y. That is about 100β_j% when β_j is small.
• Log-log. β_j is an elasticity. 1% more x goes with β_j% more y.
• Interactions. With x₁x₂ in the model, the effect of x₁ depends on x₂. The main effect is the effect at x₂ = 0. Center features to make it meaningful.
• Held fixed. The Frisch-Waugh-Lovell theorem makes this precise. β_j equals the slope from regressing residualized y on residualized x_j.

A coefficient is causal only if exogeneity holds for it. Otherwise it is a conditional association.

Why it matters in practice

• CUPED is OLS. Regressing the metric on a pre-period covariate cuts variance by a factor 1 − ρ².
• Experiment analysis. Regress the outcome on treatment plus covariates. Use robust or cluster SEs.
• Feature attribution. Correlated features make linear weights meaningless as importance scores.
• Calibration and stacking. Linear models on top of model scores remain a common, well-understood final layer.

Interview check

• Derive the OLS estimator. Set the gradient of RSS to zero. The normal equations give β̂ = (X^TX)⁻¹X^Ty.
• What does heteroskedasticity break? Not the estimate. It breaks the SEs and tests. Use the sandwich estimator.
• Is normality needed for Gauss-Markov? No. Only linearity, full rank, exogeneity, and spherical errors.

7. Generalized linear models

Three parts

• Random part. Y follows an exponential family: normal, Bernoulli, binomial, Poisson, gamma.
• Linear predictor. η = Xβ.
• Link. g(μ) = η, where μ = E[Y | X]. The canonical link makes η the natural parameter of the family.

The family fixes the mean-variance link. Var(Y) = φ V(μ). Normal has V = 1. Bernoulli has V = μ(1 − μ). Poisson has V = μ.

Logistic regression

log( p / (1 − p) ) = xTβ        p = 1 / (1 + e−xTβ)

ℓ(β) = ∑i [ yi log pi + (1 − yi) log(1 − pi) ]
∇ℓ  = XT(y − p)          Hessian = −XTWX,   W = diag(pi(1 − pi))

The log-likelihood is concave, so there is one optimum. e^β_j is an odds ratio. One unit of x_j multiplies the odds by e^β_j. It does not multiply the probability. Perfect separation drives β to infinity. Regularization or Firth's correction fixes it.

Poisson regression

log μ = xTβ + log(exposure)        ∇ℓ = XT(y − μ)

e^β_j is a rate ratio. The offset log(exposure) turns counts into rates, such as clicks per hour of use. Real counts are often overdispersed, with variance above the mean. Then Poisson SEs are too small. Use quasi-Poisson, robust SEs, or a negative binomial model.

Fitting: IRLS

Newton's method on a GLM becomes a loop of weighted least squares. Each step solves β = (X^TWX)⁻¹X^TWz. Here z is a working response, the linearized outcome. With the canonical link, Newton and Fisher scoring are the same.

Deviance

D = 2 ( ℓsaturated − ℓmodel )

Normal:     D = ∑ (yi − μi)²                      (the RSS)
Poisson:    D = 2 ∑ [ yi log(yi/μi) − (yi − μi) ]
Bernoulli:  D = −2 ∑ [ yi log pi + (1 − yi) log(1 − pi) ]   (2 × log loss)

Deviance generalizes RSS. For nested models, the drop in deviance is an LR statistic. It is about χ² with df equal to the extra parameters. Deviance divided by its residual df estimates dispersion. A ratio well above 1 flags overdispersion in Poisson and binomial counts.

Why it matters in practice

• Click and conversion models. Logistic regression is still the baseline and often the calibration layer.
• Count outcomes. Messages sent, sessions, and crashes are counts. Poisson or negative binomial fits them better than OLS.
• Deep nets are GLMs at the top. A sigmoid output with log loss is a logistic GLM on learned features.
• Odds vs probability. Stakeholders hear "doubles the odds" as "doubles the chance". Translate to probability at a typical baseline.

Interview check

• Interpret a logistic coefficient of 0.7. One unit more of x multiplies the odds by e^0.7 ≈ 2.0.
• What is deviance? Twice the log-likelihood gap to the saturated model. It plays the role of RSS for GLMs.
• Your Poisson model has deviance/df = 4. Now what? Overdispersion. Use negative binomial, quasi-Poisson, or robust SEs.

8. Bayesian inference

p(θ | x) = p(x | θ) p(θ) / p(x)   ∝   likelihood × prior

Beta-Binomial

Prior p ~ Beta(a, b). Observe k successes in n trials. The Beta density is p^a−1(1 − p)^b−1. Multiply by p^k(1 − p)^n−k and the form stays Beta.

Posterior:  p | k ~ Beta(a + k, b + n − k)
Mean:       (a + k) / (a + b + n)  =  w · (k/n) + (1 − w) · a/(a + b),    w = n/(a + b + n)

The prior acts like a + b pseudo-observations. The posterior mean is a weighted mix of the data rate and the prior rate. As n grows, the data win.

Normal-Normal

Data x_i ~ N(μ, σ²) with σ known. Prior μ ~ N(μ₀, τ₀²). Precisions add.

1/τn² = 1/τ0² + n/σ²
μn   = τn² ( μ0/τ0² + n x̄/σ² )

The posterior mean is a precision-weighted average of prior mean and sample mean. This is the engine of shrinkage. Empirical Bayes estimates μ₀ and τ₀ from many groups, then shrinks each group toward the pool.

Posterior predictive

To predict new data, average the model over the posterior.

p(x̃ | x) = ∫ p(x̃ | θ) p(θ | x) dθ

Beta-Binomial:  P(next trial succeeds) = (a + k) / (a + b + n)
Normal-Normal:  x̃ | x ~ N( μn, σ² + τn² )

The predictive variance adds two terms. One is noise in the data, σ². The other is doubt about the parameter, τ_n². Plug-in prediction drops the second and is overconfident. Posterior predictive checks simulate fake data and compare it with the real data. They are the Bayesian way to test model fit.

Credible vs confidence intervals

• Credible interval. P(θ ∈ [L, U] | data) = 0.95. A direct statement about θ, given the prior and model.
• Confidence interval. P([L, U] covers θ) = 0.95 over repeated samples. A statement about the procedure.
• When they agree. With flat priors and large n, they are numerically close. The Bernstein-von Mises theorem says the posterior becomes N(θ̂_MLE, I⁻¹/n).
• When they differ. Small n, strong priors, or parameters near a boundary.

MCMC intuition

Most posteriors have no closed form, because p(x) is an intractable integral. Markov chain Monte Carlo avoids it. It builds a random walk whose long-run distribution is the posterior. Then it treats the visited points as samples.

Metropolis-Hastings. From θ, propose θ' from q(θ' | θ). Accept with this probability:

A = min( 1,  p(θ' | x) q(θ | θ') / ( p(θ | x) q(θ' | θ) ) )

The unknown p(x) cancels in the ratio. That is the whole trick. Uphill moves are always accepted. Downhill moves are sometimes accepted, so the chain explores the tails in proportion.

• Gibbs sampling. Update one block at a time from its full conditional. Works well with conjugate pieces.
• HMC and NUTS. Use gradients to make long, well-aimed moves. They are the default in Stan and PyMC.
• Diagnostics. Discard burn-in. Run several chains and check R̂ is near 1.00. Report effective sample size, since draws are autocorrelated.
• Alternative. Variational inference fits a simple distribution by optimization. It is faster but tends to understate the variance.

Why it matters in practice

• Small segments. Shrink conversion rates for new ads or small markets toward a group prior. Raw rates from 3 clicks mislead.
• Bandits. Thompson sampling draws from Beta posteriors to balance explore and exploit.
• Decision framing. "There is a 92% chance B beats A" is easier for stakeholders than a p-value. It still needs a sound prior.
• Regularization. MAP with a Gaussian prior is L2. With a Laplace prior it is L1.

Interview check

• Prior Beta(1, 1), 3 of 10 convert. Posterior? Beta(4, 8). Mean 4/12 ≈ 0.33. P(next converts) is also 1/3.
• Credible vs confidence? Credible is a probability about θ given the data. Confidence is a coverage rate of the method.
• Why does MCMC not need the evidence p(x)? The acceptance ratio uses a ratio of posteriors, so p(x) cancels.

9. Multiple comparisons

With m independent true nulls at α = 0.05, P(at least one false positive) = 1 − 0.95^m. For m = 20 that is 64%. Twenty metrics, ten segments, or five arms all trigger this.

Two error rates

Let V be the number of false rejections and R the total rejections.

• FWER (family-wise error rate). P(V ≥ 1). The chance of even one false discovery. Use it when one false claim is costly, like a launch decision.
• FDR (false discovery rate). E[V / max(R, 1)]. The expected share of discoveries that are false. Use it for screening many candidates.

FDR control is weaker than FWER control, so it has more power. When every null is true, the two coincide.

The procedures

Sort the p-values p₍₁₎ ≤ … ≤ p_(m).

• Bonferroni. Reject if p_i ≤ α/m. Controls FWER under any dependence by the union bound. Very conservative for large m.
• Holm. Step down. Compare p_(k) with α/(m − k + 1) for k = 1, 2, …. Stop at the first failure. Controls FWER under any dependence. It is never less powerful than Bonferroni, so there is no reason to use plain Bonferroni.
• Benjamini-Hochberg. Step up. Find the largest k with p_(k) ≤ kq/m. Reject hypotheses 1 through k. Controls FDR at level q under independence or positive dependence (PRDS).
• Benjamini-Yekutieli. BH with q replaced by q / ∑_j=1^m(1/j). Valid under any dependence, at a cost in power.

def benjamini_hochberg(pvals, q=0.05):
    m = len(pvals)
    order = sorted(range(m), key=lambda i: pvals[i])
    cutoff = 0
    for rank, i in enumerate(order, start=1):
        if pvals[i] <= rank * q / m:
            cutoff = rank          # largest rank that passes
    return sorted(order[:cutoff])

def holm(pvals, alpha=0.05):
    m = len(pvals)
    order = sorted(range(m), key=lambda i: pvals[i])
    keep = []
    for rank, i in enumerate(order):
        if pvals[i] > alpha / (m - rank):
            break                  # stop at first failure
        keep.append(i)
    return sorted(keep)

p = [0.001, 0.008, 0.012, 0.030, 0.041, 0.20, 0.45, 0.62, 0.77, 0.91]
print("Bonferroni:", [i for i, x in enumerate(p) if x <= 0.05 / len(p)])
print("Holm:      ", holm(p))
print("BH:        ", benjamini_hochberg(p))

Bonferroni: [0]
Holm:       [0]
BH:         [0, 1, 2]

Note the step-up rule in BH. A p-value above its own line still counts if a later one falls below the line.

Why it matters in practice

• Experiment scorecards. Pre-register one or two primary metrics. Apply FDR or Holm to the long tail of secondary metrics.
• Segment dives. "It works for iOS users in Brazil" after 50 cuts is noise until replicated.
• Feature screening. Testing thousands of features or genes calls for BH, not Bonferroni.
• Peeking. Checking a test daily is also multiple testing. Use sequential methods or alpha spending.

Interview check

• FWER or FDR for 15 launch metrics? FWER for the few metrics that gate the decision. FDR for exploratory ones.
• Why is Holm always at least as good as Bonferroni? Its first threshold equals Bonferroni's, and later thresholds are looser.
• Does BH work under dependence? Under positive dependence, yes. For arbitrary dependence use Benjamini-Yekutieli.

10. Nonparametric and resampling methods

Permutation tests

Under the null of no treatment effect, the group labels are arbitrary. Shuffle them many times and recompute the statistic. The p-value is the share of shuffles at least as extreme as the observed value.

• Exactness. If units are exchangeable under H₀, the test has exact size for any n. No normality is needed.
• The null. It tests that the two distributions are identical, a sharp null. A difference in variance alone can trigger it.
• Counting. Use (extreme + 1)/(B + 1). The +1 counts the observed labeling and keeps p above zero.

The bootstrap

Treat the sample as the population. Draw n points with replacement, compute the statistic, and repeat B times. The spread of the replicates estimates the standard error. Efron's plug-in idea is that the empirical distribution F̂_n approximates F.

import random, statistics as st

rng = random.Random(7)
a = [rng.expovariate(1 / 10.0) for _ in range(200)]   # control
b = [rng.expovariate(1 / 12.5) for _ in range(200)]   # treatment
obs = st.mean(b) - st.mean(a)

# Permutation test: shuffle labels under H0 (no difference).
pool, n, B = a + b, len(a), 5000
extreme = 0
for _ in range(B):
    rng.shuffle(pool)
    d = st.mean(pool[n:]) - st.mean(pool[:n])
    extreme += abs(d) >= abs(obs)
p_perm = (extreme + 1) / (B + 1)

# Percentile bootstrap CI for the difference in means.
boots = []
for _ in range(B):
    ra = [rng.choice(a) for _ in a]
    rb = [rng.choice(b) for _ in b]
    boots.append(st.mean(rb) - st.mean(ra))
boots.sort()
lo, hi = boots[int(0.025 * B)], boots[int(0.975 * B) - 1]

print(f"observed diff = {obs:.2f}")
print(f"permutation p = {p_perm:.4f}")
print(f"95% bootstrap CI = ({lo:.2f}, {hi:.2f})")

observed diff = 2.99
permutation p = 0.0058
95% bootstrap CI = (0.91, 5.17)

The two methods agree here. The permutation test answers "is there any effect?" The bootstrap answers "how big is it, and how sure are we?"

When the bootstrap fails

• Extremes. The sample max cannot exceed the observed max. The bootstrap is inconsistent for it.
• Infinite variance. Very heavy tails break the bootstrap of the mean.
• Dependence. Time series and clustered data need block or cluster bootstrap. Resample whole blocks or whole users.
• Tiny samples. With n = 10, F̂_n is a poor stand-in for F.

The jackknife

Leave out one point at a time. Compute θ̂_(i) on each of the n subsets, and let θ̄_(·) be their mean.

Biaŝ   = (n − 1)(θ̄(·) − θ̂)
Var̂    = ((n − 1)/n) ∑i (θ̂(i) − θ̄(·))²
θ̂jack = θ̂ − Biaŝ           removes O(1/n) bias

The jackknife is deterministic and needs only n refits. It is a linear approximation to the bootstrap. It fails for non-smooth statistics like the median. It also supplies the acceleration constant for BCa.

Why it matters in practice

• Complex metrics. Bootstrap CIs for NDCG, recall at k, or calibration error, resampling queries or users.
• Small experiments. Permutation tests give exact p-values for small B2B or geo tests with few units.
• Randomization inference. Re-randomize the real assignment process. This respects stratification and switchback designs.
• Cost. Poisson bootstrap weights let you bootstrap in one streaming pass over billions of rows.

Interview check

• Permutation test vs bootstrap? Permutation tests a null by shuffling labels. The bootstrap estimates sampling spread by resampling.
• How do you bootstrap clustered data? Resample whole clusters, such as users, with replacement. Keep each cluster intact.
• Name a statistic the bootstrap gets wrong. The sample maximum, or the mean under infinite variance.

11. Common pitfalls

p-hacking

Analysts try many choices until p < 0.05 appears. They switch metrics, drop outliers, add covariates, cut segments, or stop early. Each choice alone looks harmless. Together they inflate the false positive rate far past α. This is the "garden of forking paths".

• Fix. Pre-register the metric, test, and stopping rule. Report all analyses run. Replicate surprising wins.
• Tell. A pile of p-values just under 0.05 is a warning sign.

Regression to the mean

Units picked for extreme values will look less extreme next time. Part of the extreme was luck, and luck does not repeat. If test and retest correlate at ρ, the expected next value is only ρ times as far from the mean.

• Example. The worst-performing ad campaigns get a fix and "improve". The worst ones would have improved anyway.
• Fix. Use a randomized control group drawn from the same extreme set.

Survivorship bias

You only see units that passed a filter. Conclusions then describe the survivors, not the population.

• Example. Long-tenure users love feature X. But users who disliked it already churned.
• Example. A model trained only on approved loans never sees the outcome of rejected ones. That is selection bias in the labels.
• Fix. Define the cohort at the start, before the filter. Follow everyone, and count dropouts as outcomes.

Ecological fallacy

A pattern across groups need not hold for individuals inside them. Countries with more internet use may have higher income. That does not prove a given user earns more. Simpson's paradox is the extreme case, where the sign flips on aggregation.

• Fix. Analyze at the level you want to make claims about. Use multilevel models when you need both levels.

Other traps worth naming

• Base rate neglect. A precise classifier on a rare event still gives mostly false alarms.
• Absence of evidence. A non-significant result with low power is not proof of no effect. Report the CI.
• Winner's curse. The measured lift of a chosen winner is biased upward. Expect shrinkage after launch.
• Unit mismatch. Randomize by user but analyze by event, and your SEs are far too small.

Why it matters in practice

• Launch reviews. Reviewers look for these errors first. Naming them yourself builds trust.
• Offline to online gaps. Survivorship in logged data explains many models that win offline and lose online.
• Long-term effects. Novelty effects and regression to the mean both make early lift fade.
• Training data. Feedback loops filter which items get exposure, so labels carry selection bias.

Interview check

• Top 10% of stores by sales got a new manager and sales fell. Did the manager hurt? Not shown. Regression to the mean predicts a drop. You need a control group.
• How do you guard against p-hacking on a team? Pre-registered metrics, fixed stopping rules, multiplicity correction, and replication.
• What is the ecological fallacy? Inferring individual behavior from group-level correlations.