Sources and References

Textbooks — Efron, B. and Tibshirani, R. J. (1993), An Introduction to the Bootstrap; Davison, A. C. and Hinkley, D. V. (1997), Bootstrap Methods and Their Application; Wasserman, L. (2004), All of Statistics Supplementary — Hastie, Tibshirani, and Friedman, The Elements of Statistical Learning; James, Witten, Hastie, and Tibshirani, An Introduction to Statistical Learning Online resources — Stanford STATS 305 materials; CMU 36-462 lecture notes

Chapter 1: Populations and Attributes

1.1 Finite Populations and Target versus Study Populations

In computational statistics we take a population-centric view of data analysis. A finite population is a well-defined collection of \(N\) units, written \(\mathcal{U} = \{u_1, u_2, \ldots, u_N\}\). Each unit carries one or more measured values called variates. If the variate of interest on unit \(u_i\) is \(y_i\), the population is fully described by the multiset \(\{y_1, y_2, \ldots, y_N\}\).

When the study population is a proper subset of the target population, any inference we make is valid only for the study population unless additional assumptions bridge the gap. A classic example: a telephone survey targets all adults in a country, but its study population is restricted to those who own telephones and answer calls. Awareness of this distinction is essential before any statistical analysis begins.

1.2 Population Attributes as Functions

An attribute of a population is any real-valued function of the population values. More precisely, if we denote the population data by the vector \(\mathbf{y} = (y_1, \ldots, y_N)\), then an attribute is any mapping \(A : \mathbb{R}^N \to \mathbb{R}\).

We distinguish between two kinds of attributes:

• Explicit attributes are defined by a closed-form formula applied directly to the population values, such as the mean, variance, or any moment.
• Implicit attributes are defined as the solution to an optimization problem or estimating equation. For instance, the population regression coefficients are the minimizers of \(\sum_{i=1}^{N}(y_i - \mathbf{x}_i^T\boldsymbol{\beta})^2\).

This functional perspective unifies seemingly disparate statistical quantities under one framework and makes it natural to study how attributes change when we observe only a sample from the population.

1.3 Graphical Attributes

Graphical attributes describe the shape of the distribution of population values. Although they are not single numbers, they are still functions of the population data and thus count as attributes in the broad sense.

Histograms

A histogram partitions the range of the data into bins and displays the frequency (or relative frequency) of observations in each bin. Given bin boundaries \(b_0 < b_1 < \cdots < b_K\) of equal width \(h = b_k - b_{k-1}\), the histogram estimator of the density at a point \(x \in [b_k, b_{k+1})\) is

\[ \hat{f}_{\text{hist}}(x) = \frac{n_k}{N \cdot h}, \]

where \(n_k\) is the number of population values falling in the \(k\)-th bin. The choice of bin width \(h\) has a profound effect on the appearance of the histogram. Too small and the histogram is noisy; too large and genuine features are smoothed away.

Kernel Density Estimation

Kernel density estimation (KDE) provides a smoother alternative to the histogram. Given population values \(y_1, \ldots, y_N\), the kernel density estimator is

\[ \hat{f}(x) = \frac{1}{Nh}\sum_{i=1}^{N} K\!\left(\frac{x - y_i}{h}\right), \]

where \(K\) is a kernel function (a symmetric, non-negative function integrating to one) and \(h > 0\) is the bandwidth. Common kernels include the Gaussian kernel \(K(u) = \frac{1}{\sqrt{2\pi}}e^{-u^2/2}\) and the Epanechnikov kernel \(K(u) = \frac{3}{4}(1-u^2)\) for \(|u| \leq 1\).

Silverman’s rule of thumb provides a quick bandwidth estimate for approximately normal data:

\[ h_{\text{Silverman}} = 0.9 \cdot \min\!\left(\hat{\sigma},\; \frac{\text{IQR}}{1.34}\right) \cdot N^{-1/5}, \]

where \(\hat{\sigma}\) is the sample standard deviation and IQR is the interquartile range. This rule should be used cautiously for multimodal distributions, where it tends to over-smooth.

More sophisticated methods include cross-validation bandwidth selection, which chooses \(h\) to minimize the integrated squared error estimated by leave-one-out cross-validation:

\[ \text{CV}(h) = \int \hat{f}(x)^2\,dx - \frac{2}{N}\sum_{i=1}^{N}\hat{f}_{-i}(y_i), \]

where \(\hat{f}_{-i}\) is the KDE constructed without observation \(y_i\).

1.4 Power Transforms

When population data are skewed or have non-constant variance, power transformations can symmetrize the distribution or stabilize variance. The general family of power transforms maps each value \(y > 0\) to \(y^{\lambda}\) for some parameter \(\lambda\).

The Box-Cox Transformation

The Box-Cox family, introduced by Box and Cox (1964), handles the discontinuity at \(\lambda = 0\) gracefully:

\[ y^{(\lambda)} = \begin{cases} \dfrac{y^{\lambda} - 1}{\lambda}, & \lambda \neq 0, \\[6pt] \ln y, & \lambda = 0. \end{cases} \]

This definition ensures that \(y^{(\lambda)}\) is a continuous function of \(\lambda\) for each fixed \(y > 0\). Special cases include:

\(\lambda\)	Transformation
1	Linear (shifted): \(y - 1\)
0.5	Square root (scaled): \(2(\sqrt{y} - 1)\)
0	Logarithm: \(\ln y\)
\(-0.5\)	Reciprocal square root (scaled)
\(-1\)	Reciprocal (scaled): \(-(1/y - 1)\)

Log and Square Root Transforms

The logarithmic transform (\(\lambda = 0\)) is the most commonly used power transform. It is appropriate when the data span several orders of magnitude or when the standard deviation is roughly proportional to the mean (coefficient of variation is approximately constant).

The square root transform (\(\lambda = 0.5\)) is commonly used for count data, as it stabilizes the variance of Poisson-distributed counts. If \(Y \sim \text{Poisson}(\mu)\), then \(\text{Var}(Y) = \mu\), and the delta method gives \(\text{Var}(\sqrt{Y}) \approx 1/4\), a constant independent of \(\mu\).

Standardization

Standardization (or z-scoring) transforms each value to have zero mean and unit variance:

\[ z_i = \frac{y_i - \bar{y}}{s}, \]

where \(\bar{y}\) and \(s\) are the sample mean and standard deviation. This is not a power transform but is frequently used alongside them. Standardization is essential when combining variables measured on different scales and is a prerequisite for many multivariate methods.

Chapter 2: Robust Estimation and Implicit Attributes

2.1 The Need for Robust Methods

Classical statistical estimators such as the sample mean and ordinary least squares are optimal under ideal conditions (e.g., normally distributed errors) but can break down catastrophically in the presence of outliers or heavy-tailed distributions. Robust statistics seeks estimators that perform well across a range of distributional assumptions.

2.2 Breakdown Point

The breakdown point of an estimator quantifies the proportion of arbitrarily corrupted observations the estimator can tolerate before giving an arbitrarily bad result.

The sample mean has a breakdown point of \(1/n\) (a single outlier can make it arbitrarily large), whereas the sample median has a breakdown point of approximately 0.5 (nearly half the data must be corrupted to make it unbounded). This makes the median far more robust than the mean.

2.3 Influence Function

The influence function measures the effect on an estimator of adding a single observation at a particular value. It was introduced by Hampel (1974) and provides a local measure of robustness.

For the mean functional \(T(F) = \int x\,dF(x)\), the influence function is \(\text{IF}(x; T, F) = x - \mu\), which is unbounded as \(x \to \pm\infty\). For the median, the influence function is bounded, confirming its robustness.

A robust estimator should have a bounded influence function, meaning that no single observation, no matter how extreme, can have an arbitrarily large effect on the estimate.

2.4 M-Estimators

M-estimators (for “maximum likelihood type”) were introduced by Huber (1964) as a generalization of maximum likelihood estimation. An M-estimator \(\hat{\theta}\) of a location parameter is defined as the minimizer of

\[ \hat{\theta} = \arg\min_{\theta} \sum_{i=1}^{n} \rho(y_i - \theta), \]

where \(\rho\) is a loss function. Equivalently, \(\hat{\theta}\) solves the estimating equation

\[ \sum_{i=1}^{n} \psi(y_i - \hat{\theta}) = 0, \]

where \(\psi = \rho'\) is the derivative of the loss function.

Different choices of \(\rho\) yield different estimators:

Loss function \(\rho(r)\)	\(\psi(r) = \rho'(r)\)	Estimator
\(r^2\)	\(2r\)	Sample mean
(	r	)
Huber loss	Huber \(\psi\)	Huber M-estimator

The Huber Loss

The Huber loss function combines the best properties of squared and absolute loss:

\[ \rho_k(r) = \begin{cases} \frac{1}{2}r^2, & |r| \leq k, \\ k|r| - \frac{1}{2}k^2, & |r| > k, \end{cases} \]

where \(k > 0\) is a tuning parameter. The corresponding \(\psi\) function is

\[ \psi_k(r) = \begin{cases} r, & |r| \leq k, \\ k\,\text{sign}(r), & |r| > k. \end{cases} \]

For observations close to the center (residual \(|r| \leq k\)), the Huber loss behaves like squared-error loss, giving the efficiency of the mean. For outlying observations (\(|r| > k\)), it behaves like absolute-error loss, bounding the influence of extreme values. A common choice is \(k = 1.345\sigma\), which achieves 95% efficiency at the normal distribution while providing substantial protection against outliers.

2.5 Median and MAD

The median is the simplest robust location estimator. For robust scale estimation, the median absolute deviation (MAD) is defined as

\[ \text{MAD} = \text{median}_{i}\,|y_i - \text{median}(y_1,\ldots,y_n)|. \]

To make the MAD consistent for the standard deviation at the normal distribution, we use the scaled version:

\[ \hat{\sigma}_{\text{MAD}} = 1.4826 \cdot \text{MAD}. \]

The factor 1.4826 arises because for \(Z \sim \mathcal{N}(0,1)\), \(\text{median}(|Z|) = \Phi^{-1}(0.75) \approx 0.6745\), and \(1/0.6745 \approx 1.4826\).

The MAD has a breakdown point of 0.5, making it much more robust than the sample standard deviation (which has a breakdown point of \(1/n\)).

2.6 Implicit Attributes via Optimization

An implicitly defined attribute is one that is specified as the solution to an optimization problem rather than as an explicit formula. This concept is central to the STAT 341 framework.

This viewpoint subsumes many familiar quantities. Ordinary least squares regression coefficients are the minimizers of the sum of squared residuals. The population median minimizes the sum of absolute deviations. Quantile regression coefficients minimize an asymmetrically weighted absolute deviation loss.

Least Absolute Deviations (LAD)

Least absolute deviations regression finds the coefficients that minimize the sum of absolute residuals:

\[ \hat{\boldsymbol{\beta}}_{\text{LAD}} = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^{N} |y_i - \mathbf{x}_i^T\boldsymbol{\beta}|. \]

Unlike OLS, LAD does not have a closed-form solution and must be computed iteratively (e.g., by linear programming or iteratively reweighted least squares). LAD regression is more robust to outliers in the response variable because the absolute value loss grows linearly rather than quadratically. The fitted line passes through at least \(p\) data points (where \(p\) is the number of predictors), which makes LAD regression equivalent to median regression (the 0.5 quantile regression).

2.7 Robust Regression

Beyond LAD, several robust regression methods exist. M-estimation in regression minimizes \(\sum_{i=1}^{n}\rho(r_i/\hat{\sigma})\) where \(r_i = y_i - \mathbf{x}_i^T\boldsymbol{\beta}\) and \(\hat{\sigma}\) is a robust scale estimate. Using the Huber loss for \(\rho\) gives a method that is efficient for normal errors and resistant to outliers.

For even higher breakdown, least trimmed squares (LTS) minimizes the sum of the \(h\) smallest squared residuals:

\[ \hat{\boldsymbol{\beta}}_{\text{LTS}} = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^{h} r_{(i)}^2(\boldsymbol{\beta}), \]

where \(r_{(1)}^2 \leq r_{(2)}^2 \leq \cdots \leq r_{(n)}^2\) are the ordered squared residuals and \(h = \lfloor n/2 \rfloor + \lfloor(p+1)/2\rfloor\). LTS achieves a breakdown point of approximately 0.5, the maximum possible, at the cost of lower efficiency.

Chapter 3: Numerical Optimization for Statistics

3.1 Why Optimization Matters in Statistics

Many statistical problems reduce to optimization: maximum likelihood estimation, M-estimation, penalized regression, and Bayesian MAP estimation all require minimizing (or maximizing) an objective function. When the objective has no closed-form optimum, we must resort to iterative numerical methods.

The three methods covered in this chapter, gradient descent, Newton-Raphson, and iteratively reweighted least squares, form a hierarchy of increasing sophistication and convergence speed.

3.2 Gradient Descent

The Basic Algorithm

Gradient descent is the simplest iterative optimization method. To minimize a differentiable function \(f(\boldsymbol{\theta})\), we start from an initial guess \(\boldsymbol{\theta}^{(0)}\) and iterate:

\[ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} - \alpha_t \nabla f(\boldsymbol{\theta}^{(t)}), \]

where \(\alpha_t > 0\) is the step size (or learning rate) and \(\nabla f\) is the gradient of \(f\).

The intuition is straightforward: the negative gradient \(-\nabla f(\boldsymbol{\theta})\) points in the direction of steepest decrease of \(f\) at \(\boldsymbol{\theta}\). Each step moves in this direction, with the step size controlling how far we go.

Fixed Step Size

With a fixed step size \(\alpha_t = \alpha\), convergence is guaranteed under certain conditions. If \(f\) is convex and its gradient is Lipschitz continuous with constant \(L\) (meaning \(\|\nabla f(\mathbf{x}) - \nabla f(\mathbf{y})\| \leq L\|\mathbf{x} - \mathbf{y}\|\) for all \(\mathbf{x}, \mathbf{y}\)), then gradient descent with step size \(\alpha \leq 1/L\) converges.

Line Search

Instead of fixing \(\alpha\), we can choose it adaptively at each step. Exact line search sets

\[ \alpha_t = \arg\min_{\alpha > 0} f\!\left(\boldsymbol{\theta}^{(t)} - \alpha\,\nabla f(\boldsymbol{\theta}^{(t)})\right). \]

This is optimal in the sense that each step makes the largest possible decrease along the gradient direction, but the one-dimensional optimization may itself be expensive.

Backtracking line search is a practical compromise. Given parameters \(0 < \beta < 1\) and \(0 < c < 1\), start with \(\alpha = 1\) and repeatedly set \(\alpha \leftarrow \beta\alpha\) until the Armijo condition is satisfied:

\[ f(\boldsymbol{\theta} - \alpha\nabla f) \leq f(\boldsymbol{\theta}) - c\,\alpha\,\|\nabla f\|^2. \]

This guarantees sufficient decrease at each step while avoiding the cost of exact line search.

3.3 Newton-Raphson Method

Derivation

Newton-Raphson (or simply Newton’s method) uses second-order information to achieve much faster convergence. The idea is to approximate \(f\) locally by a quadratic (its second-order Taylor expansion) and minimize that quadratic exactly.

At the current iterate \(\boldsymbol{\theta}^{(t)}\), the second-order Taylor approximation of \(f\) is

\[ f(\boldsymbol{\theta}) \approx f(\boldsymbol{\theta}^{(t)}) + \nabla f(\boldsymbol{\theta}^{(t)})^T(\boldsymbol{\theta} - \boldsymbol{\theta}^{(t)}) + \frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}^{(t)})^T \mathbf{H}(\boldsymbol{\theta}^{(t)})(\boldsymbol{\theta} - \boldsymbol{\theta}^{(t)}), \]

where \(\mathbf{H}(\boldsymbol{\theta}) = \nabla^2 f(\boldsymbol{\theta})\) is the Hessian matrix of second partial derivatives. Minimizing this quadratic by setting the gradient to zero gives the Newton-Raphson update:

\[ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} - \bigl[\mathbf{H}(\boldsymbol{\theta}^{(t)})\bigr]^{-1}\nabla f(\boldsymbol{\theta}^{(t)}). \]

Application to Maximum Likelihood

For MLE, we maximize the log-likelihood \(\ell(\boldsymbol{\theta})\), which is equivalent to minimizing \(f = -\ell\). The Newton-Raphson update becomes

\[ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + \bigl[-\nabla^2\ell(\boldsymbol{\theta}^{(t)})\bigr]^{-1}\nabla\ell(\boldsymbol{\theta}^{(t)}). \]

Writing \(\mathbf{U}(\boldsymbol{\theta}) = \nabla\ell(\boldsymbol{\theta})\) for the score function and \(\mathcal{J}(\boldsymbol{\theta}) = -\nabla^2\ell(\boldsymbol{\theta})\) for the observed information matrix, this becomes

\[ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + \bigl[\mathcal{J}(\boldsymbol{\theta}^{(t)})\bigr]^{-1}\mathbf{U}(\boldsymbol{\theta}^{(t)}). \]

3.4 Fisher Scoring

Fisher scoring replaces the observed information matrix \(\mathcal{J}(\boldsymbol{\theta})\) with the expected (Fisher) information matrix \(\mathcal{I}(\boldsymbol{\theta}) = E[\mathcal{J}(\boldsymbol{\theta})]\):

\[ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + \bigl[\mathcal{I}(\boldsymbol{\theta}^{(t)})\bigr]^{-1}\mathbf{U}(\boldsymbol{\theta}^{(t)}). \]

The Fisher information matrix also connects to the asymptotic variance of the MLE via the Cramer-Rao bound: \(\text{Var}(\hat{\boldsymbol{\theta}}) \approx \mathcal{I}(\hat{\boldsymbol{\theta}})^{-1}\) for large samples. Thus the Fisher scoring iteration naturally produces both the MLE and an estimate of its precision.

3.5 Iteratively Reweighted Least Squares (IRLS)

IRLS is a powerful technique that reduces certain nonlinear optimization problems to a sequence of weighted least squares problems. It arises naturally in the fitting of generalized linear models (GLMs) and M-estimation.

IRLS for GLMs

Consider a GLM with response \(\mathbf{y}\), linear predictor \(\boldsymbol{\eta} = \mathbf{X}\boldsymbol{\beta}\), link function \(g(\mu_i) = \eta_i\), and variance function \(V(\mu_i)\). The Fisher scoring update can be rewritten as a weighted least squares problem:

\[ \boldsymbol{\beta}^{(t+1)} = (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}, \]

where the working weights are

\[ w_i^{(t)} = \frac{1}{V(\mu_i^{(t)})\,(g'(\mu_i^{(t)}))^2} \]

and the working response (or adjusted dependent variate) is

\[ z_i^{(t)} = \eta_i^{(t)} + (y_i - \mu_i^{(t)})\,g'(\mu_i^{(t)}). \]

Each iteration of IRLS thus solves a standard weighted least squares problem with updated weights and working responses.

IRLS for M-Estimation

For M-estimation with loss \(\rho\), the IRLS algorithm proceeds as follows. Define weights

\[ w_i^{(t)} = \frac{\psi(r_i^{(t)}/\hat{\sigma})}{r_i^{(t)}/\hat{\sigma}}, \]

where \(r_i^{(t)} = y_i - \mathbf{x}_i^T\boldsymbol{\beta}^{(t)}\) and \(\psi = \rho'\). Then update:

\[ \boldsymbol{\beta}^{(t+1)} = (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{y}. \]

Observations with large residuals receive small weights, naturally down-weighting outliers. The iterations continue until the coefficients stabilize.

3.6 Convergence Diagnostics

In practice, we need criteria to decide when an iterative algorithm has converged. Common convergence diagnostics include:

1. Absolute change in parameters: Stop when \(\|\boldsymbol{\theta}^{(t+1)} - \boldsymbol{\theta}^{(t)}\| < \varepsilon\).
2. Relative change in parameters: Stop when \(\|\boldsymbol{\theta}^{(t+1)} - \boldsymbol{\theta}^{(t)}\|/\|\boldsymbol{\theta}^{(t)}\| < \varepsilon\).
3. Change in objective: Stop when \(|f(\boldsymbol{\theta}^{(t+1)}) - f(\boldsymbol{\theta}^{(t)})| < \varepsilon\).
4. Gradient norm: Stop when \(\|\nabla f(\boldsymbol{\theta}^{(t)})\| < \varepsilon\).

It is wise to monitor multiple criteria simultaneously and to set a maximum iteration count to prevent infinite loops. One should also inspect the Hessian at convergence to confirm that the solution is a local minimum (positive definite Hessian) rather than a saddle point or maximum.

Comparison of Methods

Method	Per-iteration cost	Convergence rate	Storage
Gradient descent	\(O(np)\)	Linear, \(O(1/t)\)	\(O(p)\)
Newton-Raphson	\(O(np^2 + p^3)\)	Quadratic	\(O(p^2)\)
Fisher scoring	\(O(np^2 + p^3)\)	Quadratic	\(O(p^2)\)
IRLS	\(O(np^2 + p^3)\)	Quadratic (canonical)	\(O(np + p^2)\)

Gradient descent is cheapest per iteration but slowest overall. Newton-Raphson and Fisher scoring converge quadratically but require computing and inverting the Hessian, which is expensive when \(p\) is large. IRLS inherits the fast convergence of Newton-Raphson for GLMs while being conceptually simpler.

Chapter 4: Sampling Theory and Design

4.1 Finite Population Inference

In finite population inference, randomness arises from the sampling design rather than from a probability model for the data. The population values \(y_1, \ldots, y_N\) are regarded as fixed but unknown; what is random is which units are selected into the sample.

This design-based framework, pioneered by Neyman (1934) and developed by Horvitz and Thompson (1952), provides inferences that are valid regardless of the underlying distribution of the population values.

4.2 Simple Random Sampling

Simple random sampling without replacement (SRSWOR) is the most fundamental design. A sample of fixed size \(n\) is drawn such that every subset of size \(n\) is equally likely. There are \(\binom{N}{n}\) possible samples, each with probability \(p(s) = 1/\binom{N}{n}\).

Under SRSWOR, the sample mean \(\bar{y}_s = \frac{1}{n}\sum_{i \in s}y_i\) is an unbiased estimator of the population mean:

\[ E[\bar{y}_s] = \bar{y}_U. \]

The variance of the sample mean under SRSWOR is

\[ \text{Var}(\bar{y}_s) = \frac{S^2_U}{n}\left(1 - \frac{n}{N}\right), \]

where \(S^2_U = \frac{1}{N-1}\sum_{i=1}^{N}(y_i - \bar{y}_U)^2\) and the factor \((1 - n/N)\) is the finite population correction (fpc). When the sampling fraction \(n/N\) is small, the fpc is close to 1 and can be ignored.

4.3 Stratified Sampling

In stratified sampling, the population is divided into \(H\) non-overlapping subgroups (strata) \(\mathcal{U}_1, \ldots, \mathcal{U}_H\) with sizes \(N_1, \ldots, N_H\) and \(\sum N_h = N\). An independent simple random sample of size \(n_h\) is drawn from each stratum.

The stratified estimator of the population mean is

\[ \bar{y}_{\text{st}} = \sum_{h=1}^{H} W_h\,\bar{y}_h, \]

where \(W_h = N_h/N\) is the stratum weight and \(\bar{y}_h\) is the sample mean in stratum \(h\).

Allocation Methods

The question of how to distribute the total sample size \(n\) among the strata is called allocation:

• Proportional allocation: \(n_h = n \cdot W_h\). This ensures that each stratum is represented in proportion to its size.
• Neyman (optimal) allocation: \(n_h \propto N_h S_h\). This minimizes the variance of \(\bar{y}_{\text{st}}\) for a fixed total sample size. Strata with greater variability receive larger samples.
• Cost-optimal allocation: \(n_h \propto N_h S_h / \sqrt{c_h}\), where \(c_h\) is the cost per observation in stratum \(h\).

4.4 Inclusion Probabilities

The first-order inclusion probability of unit \(i\) is

\[ \pi_i = P(i \in s) = \sum_{s \ni i} p(s). \]

The second-order inclusion probability of units \(i\) and \(j\) is

\[ \pi_{ij} = P(i \in s \text{ and } j \in s) = \sum_{s \ni i, j} p(s). \]

These probabilities play a central role in the Horvitz-Thompson framework.

Under SRSWOR of size \(n\) from a population of \(N\):

\[ \pi_i = \frac{n}{N}, \qquad \pi_{ij} = \frac{n(n-1)}{N(N-1)} \quad (i \neq j). \]

Under stratified SRSWOR, the inclusion probabilities are stratum-specific: \(\pi_i = n_h/N_h\) for unit \(i\) in stratum \(h\).

In general, a probability proportional to size (PPS) design sets \(\pi_i \propto z_i\) where \(z_i\) is a size measure for unit \(i\). This is efficient when larger units contribute more to the population total.

4.5 The Horvitz-Thompson Estimator

The Horvitz-Thompson (HT) estimator is a general-purpose unbiased estimator of the population total that applies to any probability sampling design with known inclusion probabilities.

The HT estimator of the population mean is simply \(\hat{\bar{y}}_{\text{HT}} = \hat{\tau}_{\text{HT}}/N\).

Variance of the HT Estimator

An unbiased estimator of this variance, due to Horvitz and Thompson (1952) and Sen (1953) and Yates and Grundy (1953), is

\[ \widehat{\text{Var}}(\hat{\tau}_{\text{HT}}) = \sum_{i \in s}\sum_{j \in s}\frac{\pi_{ij} - \pi_i\pi_j}{\pi_{ij}}\cdot\frac{y_i}{\pi_i}\cdot\frac{y_j}{\pi_j}, \]

provided \(\pi_{ij} > 0\) for all pairs \(i,j\).

4.6 Design Effect

The design effect (deff) compares the variance of an estimator under a complex design to the variance that would be obtained under SRSWOR of the same sample size:

\[ \text{deff} = \frac{\text{Var}_{\text{complex}}(\hat{\bar{y}})}{\text{Var}_{\text{SRS}}(\hat{\bar{y}})}. \]

A design effect less than 1 indicates that the complex design is more efficient than SRS; a design effect greater than 1 indicates it is less efficient. Stratified sampling typically has deff < 1 (more efficient), while cluster sampling typically has deff > 1 (less efficient because units within clusters tend to be similar).

The effective sample size is \(n_{\text{eff}} = n/\text{deff}\), representing the size of an SRSWOR sample that would achieve the same precision.

4.7 Sample Size Determination

Given a desired margin of error \(e\) and confidence level \(1 - \alpha\), the required sample size under SRSWOR is approximately

\[ n = \frac{z_{\alpha/2}^2\,S_U^2}{e^2 + z_{\alpha/2}^2\,S_U^2/N}, \]

where \(z_{\alpha/2}\) is the standard normal critical value. When \(N\) is large relative to \(n\), this simplifies to \(n \approx z_{\alpha/2}^2 S_U^2 / e^2\).

Chapter 5: Statistical Inference and Hypothesis Testing

5.1 Test Statistics as Functions of Data

A test statistic is a function of the observed data that measures evidence against a null hypothesis. In the computational statistics framework, a test statistic is simply another attribute, this time of the sample rather than the population.

The p-value is the probability, under \(H_0\), of observing a test statistic as extreme as or more extreme than the one actually observed:

\[ p = P_{H_0}(T \geq t_{\text{obs}}) \]

for a one-sided test, or \(p = P_{H_0}(|T| \geq |t_{\text{obs}}|)\) for a two-sided test. Small p-values provide evidence against \(H_0\).

5.2 T-Like Measures

A T-like measure (or T-like statistic) is a test statistic that has the general form

\[ T = \frac{\text{signal}}{\text{noise}} = \frac{\hat{\theta} - \theta_0}{\widehat{\text{se}}(\hat{\theta})}, \]

where \(\hat{\theta}\) is an estimate of the parameter of interest, \(\theta_0\) is the null hypothesis value, and \(\widehat{\text{se}}(\hat{\theta})\) is an estimated standard error.

The classical Student \(t\)-statistic is the prototypical example. Under the null \(H_0: \mu = \mu_0\),

\[ T = \frac{\bar{y} - \mu_0}{s/\sqrt{n}}, \]

which has a \(t_{n-1}\) distribution if the data are normally distributed. The T-like framework generalizes this to any estimator and any standard error estimate.

T-like measures are useful because they are pivotal or approximately pivotal: their distribution under \(H_0\) does not depend (or depends only weakly) on unknown parameters. This makes it possible to compute p-values without knowing the exact distribution of the data.

5.3 Permutation Tests

Permutation tests (also called randomization tests) provide an exact, nonparametric way to test hypotheses. The key idea is that under the null hypothesis of no treatment effect, the assignment of labels to observations is exchangeable.

The procedure for a permutation test is as follows. We compute the test statistic \(T_{\text{obs}}\) on the observed data. Then, for each possible permutation of the data labels (or a large random subset of permutations), we recompute the test statistic. The permutation p-value is the proportion of permutation test statistics that are as extreme as or more extreme than \(T_{\text{obs}}\):

\[ p_{\text{perm}} = \frac{\#\{T_{\pi} \geq T_{\text{obs}}\}}{M}, \]

where \(M\) is the number of permutations considered and \(T_\pi\) denotes the test statistic under permutation \(\pi\).

5.4 The Multiple Testing Problem

When we perform \(m\) simultaneous hypothesis tests, the chance of at least one false rejection increases rapidly. If each test is conducted at level \(\alpha = 0.05\) and the tests are independent, the probability of at least one false positive is \(1 - (1-\alpha)^m\), which approaches 1 as \(m\) grows.

The FDR is a less conservative criterion than the FWER. When many hypotheses are tested, controlling the FDR at level \(q\) tolerates some false discoveries while maintaining that, on average, the proportion of false discoveries among all discoveries is at most \(q\).

5.5 Bonferroni Correction

The Bonferroni correction is the simplest method for controlling the FWER. It rejects the \(i\)-th hypothesis if \(p_i \leq \alpha/m\), where \(m\) is the number of tests.

The Bonferroni correction is valid regardless of the dependence structure among the tests, making it very general. However, it is also very conservative, especially when \(m\) is large, leading to low power.

5.6 Holm’s Method

Holm’s method (1979) is a stepwise improvement over Bonferroni that is uniformly more powerful while still controlling the FWER at level \(\alpha\). The procedure is:

1. Order the p-values: \(p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}\).
2. Find the smallest index \(k\) such that \(p_{(k)} > \alpha/(m - k + 1)\).
3. Reject hypotheses corresponding to \(p_{(1)}, \ldots, p_{(k-1)}\).

5.7 False Discovery Rate: Benjamini-Hochberg Procedure

The Benjamini-Hochberg (BH) procedure (1995) controls the FDR at level \(q\). The procedure is:

1. Order the p-values: \(p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}\).
2. Find the largest index \(k\) such that \(p_{(k)} \leq \frac{k}{m}q\).
3. Reject all hypotheses corresponding to \(p_{(1)}, \ldots, p_{(k)}\).

The BH procedure is substantially more powerful than Bonferroni or Holm when many hypotheses are tested, because it allows a controlled proportion of false discoveries rather than demanding that the probability of even one false discovery be small.

Chapter 6: Bootstrap Methods

6.1 The Bootstrap Principle

The bootstrap, introduced by Efron (1979), is a general-purpose resampling method for estimating the sampling distribution of a statistic. The fundamental insight is that the empirical distribution of the observed data can serve as a stand-in for the unknown population distribution.

The logic is as follows. In the real world, the sample \(\mathbf{y}\) is drawn from an unknown distribution \(F\), and we want to know the distribution of \(T(\mathbf{y})\). In the bootstrap world, the bootstrap sample \(\mathbf{y}^*\) is drawn from \(\hat{F}_n\), and the distribution of \(T(\mathbf{y}^*)\) approximates that of \(T(\mathbf{y})\).

The Nonparametric Bootstrap Algorithm

The algorithm proceeds as follows.

1. Compute the statistic of interest \(\hat{\theta} = T(\mathbf{y})\) on the original data.
2. For \(b = 1, 2, \ldots, B\):
   • Draw a bootstrap sample \(\mathbf{y}^{*b}\) by sampling \(n\) values with replacement from \(\mathbf{y}\).
   • Compute \(\hat{\theta}^{*b} = T(\mathbf{y}^{*b})\).
3. The collection \(\{\hat{\theta}^{*1}, \ldots, \hat{\theta}^{*B}\}\) approximates the sampling distribution of \(\hat{\theta}\).

Typical values of \(B\) are 200 to 1000 for estimating standard errors and 1000 to 10000 for constructing confidence intervals.

6.2 Bootstrap Standard Errors

The bootstrap estimate of the standard error of \(\hat{\theta}\) is the sample standard deviation of the bootstrap replications:

\[ \widehat{\text{se}}_{\text{boot}} = \sqrt{\frac{1}{B-1}\sum_{b=1}^{B}\bigl(\hat{\theta}^{*b} - \bar{\theta}^*\bigr)^2}, \]

where \(\bar{\theta}^* = \frac{1}{B}\sum_{b=1}^{B}\hat{\theta}^{*b}\).

6.3 Bootstrap Bias Estimation

The bootstrap can also estimate the bias of an estimator. The bias of \(\hat{\theta}\) for estimating a parameter \(\theta\) is \(\text{Bias}(\hat{\theta}) = E[\hat{\theta}] - \theta\).

A bias-corrected estimate is \(\hat{\theta}_{\text{bc}} = \hat{\theta} - \widehat{\text{Bias}}_{\text{boot}} = 2\hat{\theta} - \bar{\theta}^*\). However, bias correction increases variance, and it should be used judiciously. A useful rule of thumb from Efron and Tibshirani (1993) is to apply bias correction only when the estimated bias is more than 25% of the estimated standard error.

6.4 Bootstrap Confidence Intervals

Bootstrap confidence intervals are one of the most important applications of the bootstrap. Several methods exist, each with different theoretical properties.

Percentile Method

The simplest bootstrap confidence interval is the percentile interval. A \(100(1-\alpha)\%\) percentile confidence interval is

\[ \bigl[\hat{\theta}^*_{(\alpha/2)},\;\hat{\theta}^*_{(1-\alpha/2)}\bigr], \]

where \(\hat{\theta}^*_{(\alpha/2)}\) is the \(\alpha/2\) quantile of the bootstrap distribution. For example, with \(B = 2000\) and \(\alpha = 0.05\), the 95% percentile interval uses the 50th and 1950th ordered bootstrap values.

The percentile method is simple and transformation-respecting (if we apply a monotone transformation \(\phi = g(\theta)\), the percentile interval for \(\phi\) is just \([g(\hat{\theta}^*_{(\alpha/2)}), g(\hat{\theta}^*_{(1-\alpha/2)})]\)). However, it can have poor coverage when the bootstrap distribution is biased or skewed.

Basic (Pivotal) Bootstrap Interval

The basic bootstrap interval (also called the pivotal or reflection method) is based on the idea that the distribution of \(\hat{\theta}^* - \hat{\theta}\) approximates the distribution of \(\hat{\theta} - \theta\). This leads to

\[ \bigl[2\hat{\theta} - \hat{\theta}^*_{(1-\alpha/2)},\; 2\hat{\theta} - \hat{\theta}^*_{(\alpha/2)}\bigr]. \]

Note the reversal of quantiles: the upper quantile of the bootstrap distribution is subtracted from \(2\hat{\theta}\) to get the lower confidence bound, and vice versa.

BCa (Bias-Corrected and Accelerated) Interval

The BCa interval, introduced by Efron (1987), is a second-order accurate method that corrects for both bias and skewness in the bootstrap distribution. The BCa interval adjusts the percentile endpoints:

\[ \bigl[\hat{\theta}^*_{(\alpha_1)},\;\hat{\theta}^*_{(\alpha_2)}\bigr], \]

where

\[ \alpha_1 = \Phi\!\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z_{\alpha/2})}\right), \qquad \alpha_2 = \Phi\!\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{1-\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z_{1-\alpha/2})}\right). \]

The two correction parameters are:

• Bias correction \(\hat{z}_0\): estimated as \(\hat{z}_0 = \Phi^{-1}\!\left(\frac{\#\{\hat{\theta}^{*b} < \hat{\theta}\}}{B}\right)\), the proportion of bootstrap estimates below the original estimate, transformed to the normal scale.

• \[ \hat{a} = \frac{\sum_{i=1}^{n}(\bar{\theta}_{(\cdot)} - \hat{\theta}_{(-i)})^3}{6\left[\sum_{i=1}^{n}(\bar{\theta}_{(\cdot)} - \hat{\theta}_{(-i)})^2\right]^{3/2}}. \]

Bootstrap-t Interval

The bootstrap-t interval (or studentized bootstrap) is based on the bootstrap distribution of the studentized statistic \(T^* = (\hat{\theta}^* - \hat{\theta})/\widehat{\text{se}}^*\), where \(\widehat{\text{se}}^*\) is the standard error estimate computed within each bootstrap sample. Let \(t^*_{(\alpha/2)}\) and \(t^*_{(1-\alpha/2)}\) be the quantiles of the bootstrap-t distribution. The interval is

\[ \bigl[\hat{\theta} - t^*_{(1-\alpha/2)}\,\widehat{\text{se}},\;\hat{\theta} - t^*_{(\alpha/2)}\,\widehat{\text{se}}\bigr]. \]

This method is also second-order accurate. Its main drawback is that it requires computing a standard error estimate within each bootstrap sample, which may itself require an inner bootstrap loop (a “bootstrap within a bootstrap”), making it computationally expensive.

6.5 Parametric Bootstrap

The parametric bootstrap replaces the empirical distribution \(\hat{F}_n\) with a fitted parametric model \(\hat{F}_{\hat{\theta}}\). Bootstrap samples are generated from the parametric model rather than by resampling the data.

The parametric bootstrap algorithm is:

1. Fit a parametric model to obtain \(\hat{\theta}\).
2. For \(b = 1, \ldots, B\): generate \(\mathbf{y}^{*b} \sim \hat{F}_{\hat{\theta}}\) and compute \(\hat{\theta}^{*b} = T(\mathbf{y}^{*b})\).
3. Use the distribution of \(\{\hat{\theta}^{*b}\}\) for inference.

6.6 Bootstrap Hypothesis Testing

The bootstrap can also be used for hypothesis testing. To test \(H_0: \theta = \theta_0\), we construct a bootstrap version of the null distribution.

One approach is to generate bootstrap samples from a distribution that satisfies \(H_0\). For example, to test \(H_0: \mu_1 = \mu_2\) for two samples, we can center both samples (subtract their respective means and add the pooled mean), then bootstrap from the centered data. The p-value is the proportion of bootstrap test statistics exceeding the observed test statistic.

6.7 When the Bootstrap Fails

The bootstrap is remarkably versatile but not universally applicable. It can fail in several situations:

1. Extreme order statistics: For the maximum of the sample, the bootstrap distribution is inconsistent. If \(\hat{\theta} = y_{(n)}\), the bootstrap distribution of \(y^*_{(n)} - y_{(n)}\) does not converge to the distribution of \(y_{(n)} - \theta\).

2. Non-smooth statistics: Statistics that are not smooth functions of the empirical distribution (such as the mode) may have unreliable bootstrap distributions.

3. Heavy tails without finite variance: If the underlying distribution has infinite variance, the bootstrap for the sample mean is inconsistent because the central limit theorem does not apply.

4. Dependent data: The standard iid bootstrap is invalid for dependent data (time series, spatial data). Special bootstrap variants are needed, such as the block bootstrap or the stationary bootstrap.

5. Small samples: With very small \(n\), the empirical distribution is a poor approximation to \(F\), and bootstrap inference may be unreliable.

Chapter 7: Prediction and Model Assessment

7.1 Prediction Error

In supervised learning and regression, the primary goal is often prediction: given a new observation \(\mathbf{x}_{\text{new}}\), we want to predict the response \(y_{\text{new}}\). The quality of a prediction rule \(\hat{f}\) is measured by its prediction error (or generalization error).

7.2 Training Error vs. Test Error

The training error is the average loss on the data used to fit the model:

\[ \overline{\text{err}} = \frac{1}{n}\sum_{i=1}^{n}L(y_i, \hat{f}(\mathbf{x}_i)). \]

The test error (or generalization error) is the expected loss on new, independent data:

\[ \text{Err} = E[L(Y_0, \hat{f}(\mathbf{X}_0))]. \]

Training error systematically underestimates test error because the model was fit to minimize loss on the training data. More flexible models achieve lower training error but may not generalize well; this phenomenon is called overfitting.

7.3 Bias-Variance Decomposition

The expected prediction error can be decomposed into three components. For a regression problem with squared-error loss:

This decomposition reveals the fundamental bias-variance trade-off: simple models have high bias and low variance; complex models have low bias and high variance. The optimal model balances these two sources of error. The irreducible error \(\sigma^2\) sets a lower bound on prediction error that no model can beat.

7.4 Cross-Validation

Cross-validation (CV) is the most widely used method for estimating prediction error from the training data alone. Instead of fitting the model once, we repeatedly hold out part of the data for testing and train on the rest.

Leave-One-Out Cross-Validation (LOOCV)

In LOOCV, each observation is left out one at a time. For \(i = 1, \ldots, n\), fit the model on all data except observation \(i\) to get \(\hat{f}^{(-i)}\), and predict the left-out observation. The LOOCV estimate of prediction error is

\[ \text{CV}_{(n)} = \frac{1}{n}\sum_{i=1}^{n}L(y_i, \hat{f}^{(-i)}(\mathbf{x}_i)). \]

LOOCV is approximately unbiased for the true prediction error because each model is trained on \(n-1\) observations, which is almost the full sample. However, the \(n\) estimates are highly correlated (they share \(n-2\) observations), which can lead to high variance in the CV estimate.

K-Fold Cross-Validation

In \(K\)-fold CV, the data are randomly partitioned into \(K\) roughly equal-sized groups (folds) \(C_1, \ldots, C_K\). For \(k = 1, \ldots, K\), fit the model on all data outside fold \(C_k\) and predict the observations in \(C_k\):

\[ \text{CV}_{(K)} = \frac{1}{n}\sum_{k=1}^{K}\sum_{i \in C_k}L(y_i, \hat{f}^{(-C_k)}(\mathbf{x}_i)). \]

Common choices are \(K = 5\) or \(K = 10\).

Repeated Cross-Validation

To further reduce the variance of the CV estimate, we can repeat \(K\)-fold CV multiple times with different random partitions and average the results. Repeated \(K\)-fold CV typically uses 10 to 50 repetitions. While computationally more expensive, it provides more stable estimates of prediction error.

7.5 Information Criteria

Information criteria provide an alternative to cross-validation for model selection. They add a penalty to the training error to account for model complexity.

Akaike Information Criterion (AIC)

The AIC is defined as

\[ \text{AIC} = -2\ell(\hat{\boldsymbol{\theta}}) + 2p, \]

where \(\ell(\hat{\boldsymbol{\theta}})\) is the maximized log-likelihood and \(p\) is the number of estimated parameters. The first term measures fit (smaller is better) and the second penalizes complexity.

For linear regression with normal errors, the AIC simplifies to

\[ \text{AIC} = n\ln(\text{RSS}/n) + 2p, \]

where RSS is the residual sum of squares.

Bayesian Information Criterion (BIC)

The BIC (or Schwarz criterion) replaces the penalty \(2p\) with \(p\ln(n)\):

\[ \text{BIC} = -2\ell(\hat{\boldsymbol{\theta}}) + p\ln(n). \]

Since \(\ln(n) > 2\) for \(n \geq 8\), the BIC imposes a heavier penalty on model complexity than AIC for moderate to large sample sizes. This makes BIC more conservative, tending to select simpler models.

7.6 Model Selection via Cross-Validation

Cross-validation is often used directly for model selection. The procedure is:

1. For each candidate model \(m \in \mathcal{M}\), compute the CV estimate of prediction error \(\widehat{\text{Err}}_{\text{CV}}(m)\).
2. Select the model \(\hat{m} = \arg\min_{m \in \mathcal{M}} \widehat{\text{Err}}_{\text{CV}}(m)\).

The One-Standard-Error Rule

The one-standard-error rule, proposed by Breiman et al. (1984) and popularized by Hastie, Tibshirani, and Friedman (2009), selects the simplest model whose CV error is within one standard error of the minimum:

\[ \hat{m}_{1\text{SE}} = \text{simplest } m \text{ such that } \widehat{\text{Err}}_{\text{CV}}(m) \leq \widehat{\text{Err}}_{\text{CV}}(\hat{m}) + \text{SE}(\widehat{\text{Err}}_{\text{CV}}(\hat{m})). \]

This rule incorporates a preference for parsimony, recognizing that CV estimates have uncertainty and that a simpler model that performs nearly as well is often preferable.

7.7 Optimism and the Training-Test Gap

The optimism of the training error is the expected difference between the test error and the training error:

\[ \text{op} = E[\text{Err}] - E[\overline{\text{err}}]. \]

For squared-error loss and linear models, the optimism is

\[ \text{op} = \frac{2}{n}\sum_{i=1}^{n}\text{Cov}(y_i, \hat{y}_i) = \frac{2p\sigma^2}{n}, \]

where \(p\) is the number of parameters and \(\sigma^2\) is the error variance. This yields Mallows’ \(C_p\) statistic:

\[ C_p = \overline{\text{err}} + \frac{2p\hat{\sigma}^2}{n}, \]

which estimates the expected prediction error. The \(C_p\) statistic, AIC, and LOOCV are all asymptotically equivalent for linear models.

7.8 Summary of Model Assessment Strategies

The following table summarizes the key properties of model assessment methods:

Method	Type	Computational cost	Best suited for
LOOCV	Resampling	\(n\) model fits	Small datasets, linear models (shortcut)
\(K\)-fold CV	Resampling	\(K\) model fits	General purpose, moderate-large datasets
AIC	Analytic	1 model fit per candidate	Large datasets, parametric models
BIC	Analytic	1 model fit per candidate	Model identification, large \(n\)
Bootstrap (.632)	Resampling	\(B\) model fits	Small datasets, complex models

The choice among these methods depends on the sample size, computational budget, and the goal of the analysis (prediction accuracy vs. model identification). In most practical settings, 10-fold cross-validation provides a reliable, general-purpose approach to model assessment.