Chapter 35: Sampling & Estimation

The Population-Sample Distinction

A population is the complete collection of units about which we want to draw conclusions. A sample is the subset we actually observe. The population can be finite (all registered voters in a country) or conceptually infinite (all possible measurements a sensor could produce). The key insight is that we almost never observe the population directly — we make do with a sample, then reason backward.

This distinction shapes every subsequent decision. The population parameter $\theta$ (a mean $\mu$, proportion $p$, regression coefficient $\beta$, etc.) is a fixed, unknown constant. Our estimator $\hat{\theta}$ is a function of the sample; it is a random variable because it changes with each new sample we might draw.

Why Sampling Design Matters

The validity of inference depends critically on how the sample was collected. A biased sample produces biased inferences no matter how sophisticated the analysis.

Simple random sampling (SRS) gives every unit an equal probability of selection. With a sample of size $n$ from a population of size $N$, the inclusion probability is $n/N$. SRS is the baseline against which other designs are compared.

Stratified sampling partitions the population into homogeneous strata and draws a SRS within each. If stratum $h$ contains $N_h$ units and we draw $n_h$, the estimator of the population mean is the weighted average:

Stratification reduces variance when the strata differ substantially from each other (between-stratum variance is large) while being internally homogeneous (within-stratum variance is small).

Cluster sampling groups units into clusters, samples clusters, then observes all (or a subsample of) units within selected clusters. It is cost-efficient when clusters are geographically dispersed, but inflates variance because units within a cluster are often similar — captured by the intraclass correlation coefficient (ICC).

Systematic sampling selects every $k$-th unit after a random start. It is efficient in practice but can produce biased estimates if the population has periodicity that coincides with $k$.

The Convenience Sample Problem

Convenience samples — scraped web data, volunteers, users of a particular app — are the norm in industry machine learning but the nemesis of inference. They introduce selection bias: the mechanism that determines who appears in the sample is correlated with the outcome of interest. If you train a sentiment model on Amazon reviews, you oversample opinionated purchasers and your estimates of average sentiment will not generalize to the broader population.

Recognizing and documenting the sampling mechanism is not pedantry; it determines which conclusions are defensible.

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)

# Simulate a population: income (USD) stratified by education level
n_pop = 100_000
education = rng.choice(['HS', 'Bachelor', 'Graduate'], size=n_pop, p=[0.50, 0.35, 0.15])
mean_income = {'HS': 40_000, 'Bachelor': 65_000, 'Graduate': 95_000}
std_income  = {'HS': 12_000, 'Bachelor': 18_000, 'Graduate': 25_000}

pop_income = np.array([
    rng.normal(mean_income[e], std_income[e]) for e in education
])
pop_mean = pop_income.mean()  # ground truth
print(f"Population mean income: ${pop_mean:,.0f}")

# Simple random sample
srs_idx = rng.choice(n_pop, size=500, replace=False)
srs_mean = pop_income[srs_idx].mean()
print(f"SRS estimate:           ${srs_mean:,.0f}")

# Stratified sample (proportional allocation)
strata_sizes = {'HS': 250, 'Bachelor': 175, 'Graduate': 75}
strat_means = []
for stratum, n_h in strata_sizes.items():
    mask = education == stratum
    stratum_pop = pop_income[mask]
    sample = rng.choice(stratum_pop, size=n_h, replace=False)
    strat_means.append((mask.mean(), sample.mean()))  # (W_h, y_bar_h)

strat_estimate = sum(w * ybar for w, ybar in strat_means)
print(f"Stratified estimate:    ${strat_estimate:,.0f}")

Finite-Population Correction

When the sample is a non-negligible fraction of the population, the variance of the sample mean is reduced by the finite-population correction (FPC):

where $S^2$ is the population variance. When $n/N \to 0$ (large populations), the FPC approaches 1 and is routinely ignored. For small populations (surveys of a specific company's employees, for instance), ignoring the FPC overstates uncertainty.

What Is a Sampling Distribution?

If we drew every possible sample of size $n$ from the population and computed a statistic (say, the sample mean $\bar{X}$) for each, the resulting distribution of those statistics is the sampling distribution. This is not the distribution of the data — it is the distribution of an estimator over hypothetical repeated experiments.

The sampling distribution is the conceptual engine of frequentist inference. Confidence intervals and p-values are both statements about where statistics fall in their sampling distributions.

The Sampling Distribution of the Mean

Let $X_1, X_2, \ldots, X_n$ be i.i.d. draws from a population with mean $\mu$ and variance $\sigma^2$. The sample mean is:

By linearity of expectation and independence:

The standard error of the mean (SEM) is $\text{SE}(\bar{X}) = \sigma / \sqrt{n}$. It captures how much the sample mean varies from sample to sample and shrinks at rate $1/\sqrt{n}$ as sample size grows.

The Central Limit Theorem

The Central Limit Theorem (CLT) guarantees that, regardless of the population's shape (provided it has finite variance), the standardized sample mean converges in distribution to a standard normal:

This is extraordinarily powerful. It means that for moderate $n$ (often $n \geq 30$ is cited as a rough rule, though heavy tails require larger $n$), we can use normal approximations for inference without knowing the population distribution.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

rng = np.random.default_rng(42)

# Population: exponential (heavily right-skewed)
lam = 0.2  # mean = 1/lam = 5
pop = rng.exponential(scale=1/lam, size=1_000_000)

fig, axes = plt.subplots(1, 4, figsize=(14, 4))
for ax, n in zip(axes, [1, 5, 30, 100]):
    xbars = [rng.choice(pop, n).mean() for _ in range(5000)]
    ax.hist(xbars, bins=50, density=True, color='steelblue', alpha=0.7)
    mu, se = 1/lam, (1/lam) / np.sqrt(n)
    x = np.linspace(mu - 4*se, mu + 4*se, 200)
    ax.plot(x, stats.norm.pdf(x, mu, se), 'r-', lw=2)
    ax.set_title(f'n = {n}')
    ax.set_xlabel('Sample mean')

plt.suptitle('CLT in action: Exponential population, sampling distribution of mean')
plt.tight_layout()
plt.savefig('clt_demo.png', dpi=120)
print("Figure saved.")

Beyond the Mean: Other Sampling Distributions

The CLT applies broadly, but some statistics have exact finite-sample distributions:

• Sample variance: If <!--MATHBLOCK13--> i.i.d., then <!--MATHBLOCK14-->.
• t-statistic: The ratio <!--MATHBLOCK15--> exactly under normality. The <!--MATHBLOCK16--> distribution has heavier tails than the normal, accounting for the additional uncertainty from estimating <!--MATHBLOCK17-->.
• F-statistic: Ratios of independent <!--MATHBLOCK18--> random variables, arising in ANOVA and regression.

Simulating Sampling Distributions for Any Statistic

For complex statistics (medians, quantiles, ratios, ML model coefficients), the CLT may not apply cleanly or may converge slowly. Simulation-based approaches let us construct sampling distributions empirically:

import numpy as np
from scipy import stats

rng = np.random.default_rng(99)

# Population: mixture of two normals
def draw_population(size, rng):
    mask = rng.random(size) < 0.4
    return np.where(mask,
        rng.normal(10, 2, size),
        rng.normal(25, 3, size)
    )

pop = draw_population(500_000, rng)
true_median = np.median(pop)

# Sampling distribution of the sample median, n=50
medians = [np.median(rng.choice(pop, 50, replace=False)) for _ in range(5000)]
print(f"True median: {true_median:.3f}")
print(f"Mean of sampling dist: {np.mean(medians):.3f}")
print(f"Std of sampling dist (SE of median): {np.std(medians):.3f}")

The Law of Large Numbers

While the CLT describes the shape of the sampling distribution, the Law of Large Numbers (LLN) guarantees that $\bar{X} \to \mu$ as $n \to \infty$. The weak LLN says convergence is in probability; the strong LLN says convergence is almost sure. Together, the LLN and CLT form the bedrock of why estimation works: the estimator converges to the truth, and we can quantify how spread out it is around that truth for any finite $n$.

What Makes a Good Estimator?

An estimator $\hat{\theta}$ is a function of the sample data that we use to approximate the unknown parameter $\theta$. Since $\hat{\theta}$ is a random variable (it depends on which sample we got), we evaluate it through the properties of its distribution.

Bias

The bias of an estimator is the expected difference between the estimate and the true parameter:

An estimator is unbiased if $\text{Bias}(\hat{\theta}) = 0$, meaning it is correct on average over repeated sampling. The sample mean $\bar{X}$ is an unbiased estimator of $\mu$:

However, the naive variance estimator $\tilde{S}^2 = \frac{1}{n}\sum (X_i - \bar{X})^2$ is biased:

This motivates the Bessel-corrected sample variance $S^2 = \frac{1}{n-1}\sum (X_i - \bar{X})^2$, which is unbiased. The intuition: $\bar{X}$ is computed from the same data, so the deviations $X_i - \bar{X}$ underestimate the true deviations $X_i - \mu$; dividing by $n-1$ corrects for this.

Variance

The variance of an estimator measures how much it fluctuates across repeated samples:

Lower variance means more consistent estimates. For the sample mean from an i.i.d. sample, $\text{Var}(\bar{X}) = \sigma^2/n$.

Mean Squared Error

Neither bias nor variance alone tells the full story. The mean squared error (MSE) combines both:

This decomposition is critical. A biased estimator can have lower MSE than an unbiased one if it compensates with substantially lower variance — this is the bias-variance tradeoff that also governs regularization in machine learning.

import numpy as np

rng = np.random.default_rng(42)
sigma_true = 5.0
n_vals = [5, 10, 30, 100]
B = 100_000

print(f"{'n':>5}  {'E[S^2]':>10}  {'E[S_biased^2]':>15}  {'Bias^2(biased)':>16}  {'MSE(biased)':>12}  {'MSE(unbiased)':>14}")
for n in n_vals:
    samples = rng.normal(0, sigma_true, (B, n))
    s2_unbiased = samples.var(axis=1, ddof=1)  # divide by n-1
    s2_biased   = samples.var(axis=1, ddof=0)  # divide by n

    bias_sq = (s2_biased.mean() - sigma_true**2)**2
    mse_b   = np.mean((s2_biased   - sigma_true**2)**2)
    mse_u   = np.mean((s2_unbiased - sigma_true**2)**2)

    print(f"{n:>5}  {s2_unbiased.mean():>10.4f}  {s2_biased.mean():>15.4f}  "
          f"{bias_sq:>16.4f}  {mse_b:>12.4f}  {mse_u:>14.4f}")

Consistency

An estimator is consistent if it converges in probability to the true parameter as $n \to \infty$: $\hat{\theta}_n \xrightarrow{p} \theta$. Equivalently, $\text{MSE}(\hat{\theta}_n) \to 0$. Consistency requires that both bias and variance vanish as the sample grows.

The sample mean is consistent for $\mu$ (this follows directly from the LLN). The biased sample variance $\tilde{S}^2$ is also consistent for $\sigma^2$ because its bias $(\frac{-\sigma^2}{n})$ vanishes as $n \to \infty$.

Efficiency and the Cramér-Rao Lower Bound

Among all unbiased estimators, is there a best one — one with the lowest possible variance? The Cramér-Rao Lower Bound (CRLB) answers this:

where $I(\theta)$ is the Fisher information:

An estimator that achieves this bound is called efficient. The sample mean is efficient for estimating the mean of a normal distribution — no other unbiased estimator can do better.

Asymptotic Normality

Many estimators satisfy the asymptotic normality property:

for some asymptotic variance $V$. This is the foundation for large-sample inference: even without knowing the exact sampling distribution, we can approximate it as normal and construct confidence intervals.

The Likelihood Principle

Maximum likelihood estimation (MLE) is the dominant method for deriving point estimators. The idea is intuitive: choose the parameter value that makes the observed data most probable.

Given i.i.d. observations $x_1, \ldots, x_n$ from a model with density/pmf $f(x; \theta)$, the likelihood function is:

The MLE is $\hat{\theta}_{\text{MLE}} = \arg\max_\theta L(\theta)$. Working with the log-likelihood $\ell(\theta) = \log L(\theta) = \sum_i \log f(x_i; \theta)$ is equivalent (log is monotone) and numerically superior since it converts products to sums.

Deriving MLEs Analytically

Normal Distribution

For $X_i \sim N(\mu, \sigma^2)$:

Setting $\partial \ell / \partial \mu = 0$ gives $\hat{\mu} = \bar{x}$ (the sample mean).

Setting $\partial \ell / \partial \sigma^2 = 0$ gives $\hat{\sigma}^2 = \frac{1}{n}\sum(x_i - \bar{x})^2$ — the biased sample variance. MLEs are not guaranteed to be unbiased, though they are consistent and asymptotically efficient.

Exponential Distribution

For $X_i \sim \text{Exponential}(\lambda)$ with density $f(x; \lambda) = \lambda e^{-\lambda x}$:

The MLE of the rate is the reciprocal of the sample mean — the reciprocal of the MLE of the mean.

The Invariance Property

If $\hat{\theta}$ is the MLE of $\theta$, then for any function $g$, the MLE of $g(\theta)$ is $g(\hat{\theta})$. This is the invariance property of MLEs, and it is extremely convenient: the MLE of $\sigma$ (not $\sigma^2$) is $\sqrt{\hat{\sigma}^2}$.

Numerical MLE

When the log-likelihood has no closed-form maximizer, we optimize numerically.

import numpy as np
from scipy import optimize, stats

rng = np.random.default_rng(42)

# Generate data from a Gamma distribution
alpha_true, beta_true = 3.5, 1.5  # shape, rate
data = rng.gamma(shape=alpha_true, scale=1/beta_true, size=300)

def neg_log_lik(params):
    alpha, beta = params
    if alpha <= 0 or beta <= 0:
        return np.inf
    return -np.sum(stats.gamma.logpdf(data, a=alpha, scale=1/beta))

result = optimize.minimize(
    neg_log_lik,
    x0=[2.0, 1.0],          # initial guess
    method='Nelder-Mead',
    options={'xatol': 1e-6, 'fatol': 1e-8, 'maxiter': 10000}
)
alpha_mle, beta_mle = result.x
print(f"True:  alpha={alpha_true}, beta={beta_true}")
print(f"MLE:   alpha={alpha_mle:.4f}, beta={beta_mle:.4f}")

# Compare with scipy's built-in fit (uses MLE internally)
a_fit, loc_fit, scale_fit = stats.gamma.fit(data, floc=0)
print(f"scipy: alpha={a_fit:.4f}, beta={1/scale_fit:.4f}")

Properties of MLEs

Under regularity conditions, MLEs enjoy three asymptotic properties:

1. Consistency: <!--MATHBLOCK24--> as <!--MATHBLOCK25-->.
2. Asymptotic normality: <!--MATHBLOCK26-->, where <!--MATHBLOCK27--> is the Fisher information.
3. Asymptotic efficiency: No consistent estimator has a smaller asymptotic variance.

This means that for large samples, we can approximate the distribution of any MLE as normal with variance equal to the inverse Fisher information — even without knowing the exact finite-sample distribution.

Observed vs. Expected Fisher Information

In practice, we estimate the Fisher information using the observed Fisher information — the negative Hessian of the log-likelihood evaluated at the MLE:

The estimated standard error of the MLE is then $\widehat{\text{SE}}(\hat{\theta}) = 1/\sqrt{\hat{I}(\hat{\theta})}$ for scalar $\theta$. For vector parameters, the covariance matrix is approximated by $[-\nabla^2 \ell(\hat{\theta})]^{-1}$.

import numpy as np
from scipy import optimize, stats

rng = np.random.default_rng(7)
data = rng.exponential(scale=4.0, size=200)  # true lambda = 0.25

def log_lik(lam):
    return np.sum(np.log(lam) - lam * data)

# MLE
lam_mle = 1.0 / data.mean()
print(f"MLE: lambda = {lam_mle:.5f}")

# Observed Fisher information (second derivative of log-lik)
# d^2/dlam^2 log L = -n / lam^2
obs_fisher = len(data) / lam_mle**2
se_mle = 1.0 / np.sqrt(obs_fisher)
print(f"SE(lambda_MLE) = {se_mle:.5f}")
print(f"Approx 95% CI: ({lam_mle - 1.96*se_mle:.4f}, {lam_mle + 1.96*se_mle:.4f})")

Common Pitfalls

• Local maxima: Numerical optimizers may converge to a local maximum. Always try multiple starting points.
• Boundary solutions: If the MLE lies on the boundary of the parameter space (e.g., variance estimate of zero), standard asymptotic theory breaks down.
• Model misspecification: MLE is consistent for <!--MATHBLOCK31--> only if the model is correctly specified. Under misspecification, MLE converges to the "pseudo-true" parameter that minimizes the Kullback-Leibler divergence — which may not be the true parameter.

Standard Errors

The standard error of an estimator is the standard deviation of its sampling distribution:

It answers: "How much would my estimate change if I drew a new sample?" A small SE means the estimator is tightly concentrated around its expected value; a large SE means estimates vary substantially across samples.

For the sample mean, $\text{SE}(\bar{X}) = \sigma/\sqrt{n}$, estimated by $\hat{\text{SE}} = S/\sqrt{n}$ where $S$ is the sample standard deviation. Note carefully: the SE shrinks as $\sqrt{n}$, so to halve the standard error you must quadruple the sample size.

The Bootstrap

For complex statistics (medians, correlation coefficients, AUC, model parameters), deriving analytic SEs is often intractable. The bootstrap, introduced by Efron (1979), provides a general simulation-based approach.

The plug-in principle: estimate the sampling distribution of $\hat{\theta}$ by treating the sample as a proxy population and resampling from it.

Algorithm — Nonparametric Bootstrap:

1. Compute <!--MATHBLOCK9--> from the original sample of size <!--MATHBLOCK10-->.
2. For <!--MATHBLOCK11-->: draw a bootstrap sample of size <!--MATHBLOCK12--> with replacement from the data, compute <!--MATHBLOCK13-->.
3. Use the distribution of <!--MATHBLOCK14--> to estimate the SE and construct confidence intervals.

The bootstrap SE is simply the standard deviation of the bootstrap replicates:

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
data = rng.lognormal(mean=2.0, sigma=0.8, size=80)

# Statistic of interest: ratio of mean to median
def statistic(x):
    return x.mean() / np.median(x)

theta_hat = statistic(data)

# Bootstrap
B = 5000
boot_reps = np.array([
    statistic(rng.choice(data, len(data), replace=True))
    for _ in range(B)
])

se_boot = boot_reps.std(ddof=1)
print(f"Estimate (mean/median):  {theta_hat:.4f}")
print(f"Bootstrap SE:            {se_boot:.4f}")
print(f"Bootstrap 95% CI (percentile): ({np.percentile(boot_reps, 2.5):.4f}, {np.percentile(boot_reps, 97.5):.4f})")

Confidence Intervals

A confidence interval (CI) is an interval estimator $[L(X), U(X)]$ constructed such that it covers the true parameter with a specified probability:

The Correct Interpretation

This is one of the most misunderstood concepts in statistics. The confidence level $1-\alpha$ is a property of the procedure, not of a particular interval. Before the data are collected, 95% of intervals constructed this way will contain $\theta$. After the data are observed and the interval is computed, the interval either contains $\theta$ or it does not — there is no probability statement about a specific realized interval.

Wrong: "There is a 95% probability that $\theta$ is between 2.1 and 3.7."

Right: "This interval was constructed using a procedure that covers the true parameter in 95% of repeated applications."

Constructing Confidence Intervals

Pivot method (exact, when available): Find a function of the data and the parameter — a pivot — whose distribution does not depend on unknown parameters, then invert to isolate the parameter.

For the normal mean with unknown variance, the pivot is $T = (\bar{X} - \mu)/(S/\sqrt{n}) \sim t_{n-1}$, giving the t-interval:

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
data = rng.normal(loc=10.0, scale=3.0, size=25)
n = len(data)
xbar, s = data.mean(), data.std(ddof=1)
t_crit = stats.t.ppf(0.975, df=n-1)
ci = (xbar - t_crit * s/np.sqrt(n), xbar + t_crit * s/np.sqrt(n))
print(f"Sample mean: {xbar:.3f}, S: {s:.3f}")
print(f"95% t-interval: ({ci[0]:.3f}, {ci[1]:.3f})")
print(f"True mean 10.0 is {'inside' if ci[0] <= 10.0 <= ci[1] else 'outside'} the CI")

Bootstrap confidence intervals: Three common types:

• Percentile method: Use the <!--MATHBLOCK21--> and <!--MATHBLOCK22--> quantiles of the bootstrap distribution directly.
• Basic (reflected) method: <!--MATHBLOCK23-->. Corrects for asymmetry.
• BCa (bias-corrected and accelerated): Adjusts for both bias and skewness in the bootstrap distribution — the most reliable but computationally intensive.

Coverage Probability Simulation

The true test of a CI method is its coverage probability — what fraction of constructed intervals actually contain the true parameter:

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
mu_true = 5.0
n = 20
alpha = 0.05
B_boot = 500
reps = 2000

t_cover, boot_cover = 0, 0

for _ in range(reps):
    data = rng.exponential(scale=mu_true, size=n)  # skewed data
    xbar, s = data.mean(), data.std(ddof=1)
    t_crit = stats.t.ppf(1 - alpha/2, df=n-1)
    t_lo, t_hi = xbar - t_crit*s/np.sqrt(n), xbar + t_crit*s/np.sqrt(n)
    t_cover += (t_lo <= mu_true <= t_hi)

    boot_reps = [rng.choice(data, n, replace=True).mean() for _ in range(B_boot)]
    b_lo, b_hi = np.percentile(boot_reps, [2.5, 97.5])
    boot_cover += (b_lo <= mu_true <= b_hi)

print(f"t-interval  coverage: {t_cover/reps:.3f} (nominal: 0.950)")
print(f"Boot pctile coverage: {boot_cover/reps:.3f} (nominal: 0.950)")
print("Note: Both under-cover for skewed data at n=20.")

Pitfalls and Practical Guidance

• Multiple testing: Constructing many CIs inflates the chance of at least one failing to cover the true value. Use Bonferroni or simultaneous CIs when comparing many parameters.
• Bootstrap assumptions: The bootstrap relies on the sample adequately representing the population. It fails badly for statistics that depend on extreme order statistics (min, max) and for heavily dependent data (use block bootstrap for time series).
• Interpreting width: A wide CI reflects genuine uncertainty, not a flaw in the analysis. Reporting a point estimate without an SE or CI obscures this uncertainty and is bad practice.