Resampling & Permutation Methods: Bootstrap and Jackknife

Let $\mathbf{X} = (X_1, X_2, \dots, X_n)$ be an independent and identically distributed (i.i.d.) random sample drawn from an unknown cumulative distribution function $F$:

We desire to estimate an unknown population parameter $\theta = t(F)$, such as the population mean $\mu = \int x\,dF(x)$, population median $F^{-1}(0.5)$, or variance $\sigma^2$.

Let $\hat{\theta} = s(\mathbf{X})$ be a sample estimator of $\theta$. To quantify the statistical uncertainty of $\hat{\theta}$, we need its sampling distribution $F_{\hat{\theta}}$, standard error $\text{se}_F(\hat{\theta})$, and bias $\text{bias}_F(\hat{\theta}) = \mathbb{E}_F[\hat{\theta}] - \theta$.

1. The Empirical Distribution Function (EDF)

The sample data $\mathbf{X}$ define the Empirical Distribution Function (EDF) $\hat{F}_n$:

The EDF assigns discrete probability mass $1/n$ to each of the $n$ observed data points:

By the Glivenko-Cantelli Theorem (The Fundamental Theorem of Statistics):

As $n \to \infty$, the empirical distribution $\hat{F}_n$ converges uniformly to the true population distribution $F$.

2. The Plug-In Principle

The Plug-In Principle dictates that any unknown population characteristic $\theta = t(F)$ is estimated by the corresponding characteristic of the empirical distribution:

Similarly, the unknown standard error $\text{se}_F(\hat{\theta})$ is estimated by the plug-in estimate:

3. The Non-Parametric Bootstrap Algorithm

Because calculating $\text{se}_{\hat{F}_n}(\hat{\theta}^*)$ analytically is intractable for complex estimators, Efron approximated it via Monte Carlo simulation:

1. For $b = 1, 2, \dots, B$ (where $B \ge 1,000$ bootstrap replicates):

• Draw a bootstrap sample $\mathbf{X}^{*b} = (X_1^{*b}, X_2^{*b}, \dots, X_n^{*b})$ by sampling $n$ observations uniformly at random with replacement from the original dataset $\mathbf{X}$.
• Compute the bootstrap replicate of the statistic:

1. The Bootstrap Standard Error is the sample standard deviation of the $B$ bootstrap replications:

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

1. The Bootstrap Bias Estimate is:

The bias-corrected point estimate is $\hat{\theta}_{\text{bc}} = \hat{\theta} - \widehat{\text{bias}}_B = 2\hat{\theta} - \bar{\theta}^*$.

Constructing confidence intervals for parameter $\theta$ via bootstrapping admits several formulations:

1. The Standard Normal Bootstrap Interval

Assuming the sampling distribution of $\hat{\theta}$ is approximately Gaussian:

Limitation: Fails when the sampling distribution is skewed, heavy-tailed, or bounded.

2. The Percentile Bootstrap Interval

Sort the $B$ bootstrap estimates in ascending order:

The $(1 - \alpha) \times 100\%$ Percentile Interval uses the empirical quantiles of the bootstrap distribution:

For a $95\%$ confidence interval with $B = 2,000$:

Key Property (Transformation Invariance): If $g(\cdot)$ is any strictly monotonic transformation, the percentile interval for $g(\theta)$ is exactly $[g(\hat{\theta}_{\text{low}}), g(\hat{\theta}_{\text{high}})]$. Accuracy: First-order accurate: coverage error $= \mathcal{O}(n^{-1/2})$.

3. The Basic (Empirical) Bootstrap Interval

Let the estimation error be $\delta = \hat{\theta} - \theta$. The basic bootstrap approximates the distribution of $\delta$ by the bootstrap errors $\delta^* = \hat{\theta}^* - \hat{\theta}$:

4. Bias-Corrected and Accelerated (BCa) Intervals

Formulated by Efron in 1987, the BCa Interval corrects for both median bias and variance stabilization through two parameters:

• $z_0$ (Bias-Correction Factor): Measures the median bias of $\hat{\theta}^*$:

• $a$ (Acceleration Parameter): Measures the rate of change of the standard error of $\hat{\theta}$ with respect to the true parameter $\theta$, computed via jackknife influence values:

The adjusted quantile points $\alpha_1, \alpha_2$ are:

The BCa interval is second-order accurate: coverage error $= \mathcal{O}(n^{-1})$ and is transformation-invariant.

Developed by Maurice Quenouille (1949) for bias reduction and extended by John Tukey (1958) for variance estimation, the Jackknife is a deterministic leave-one-out resampling method that requires no random number generation.

1. Leave-One-Out Estimators and Pseudovalues

Let $\mathbf{X} = (X_1, \dots, X_n)$ be a sample of size $n$, with point estimate $\hat{\theta} = s(\mathbf{X})$. For each observation $i \in \{1, 2, \dots, n\}$, define the leave-one-out subsample $\mathbf{X}_{(i)}$ omitting the $i$-th data point:

Let $\hat{\theta}_{(i)} = s(\mathbf{X}_{(i)})$ be the statistic evaluated on this $(n-1)$-sized subsample.

The mean of the leave-one-out estimates is:

2. Jackknife Bias and Variance Formulas

• Jackknife Bias Estimator:

• Bias-Corrected Jackknife Estimator:

• Jackknife Variance Estimator:

• Jackknife Standard Error:

Theoretical Limitation: The jackknife is non-smooth: it fails completely for non-differentiable statistics such as the sample median or sample quantiles, where $\widehat{\text{var}}_{\text{jack}}$ is inconsistent. The bootstrap, by contrast, handles non-smooth statistics seamlessly.

When testing the null hypothesis $H_0$ that two samples come from identical distributions, Ronald Fisher formulated the Permutation Test (Randomization Test).

The Two-Sample Permutation Algorithm

Consider sample $\mathbf{X} = (X_1, \dots, X_{n_1})$ from group 1 and sample $\mathbf{Y} = (Y_1, \dots, Y_{n_2})$ from group 2.

• Null hypothesis: $H_0: F_X = F_Y$.
• Compute the observed test statistic (e.g., difference in sample means):

• Under $H_0$, group membership labels are completely exchangeable. Pool the $N = n_1 + n_2$ observations into a single dataset $\mathbf{Z} = (\mathbf{X}, \mathbf{Y})$.
• The total number of possible partitions of $N$ items into groups of size $n_1$ and $n_2$ is:

For each of the $M$ partitions $\pi$:

• Assign $n_1$ items to group 1 and $n_2$ items to group 2, and compute $T_\pi = \bar{X}_\pi - \bar{Y}_\pi$.

The Exact Permutation p-value (for a two-tailed test) is:

When $M$ is astronomically large (e.g., $n_1 = n_2 = 50 \implies \binom{100}{50} \approx 10^{29}$), exact enumeration is impossible. Statisticians draw a random sample of $B = 10,000$ random permutations to compute the Monte Carlo Permutation p-value:

Permutation tests are completely distribution-free, valid for any sample size, and require zero parametric assumptions.

In Section 8.4, we examined the Sequential Importance Resampling (SIR) particle filter. While SIR is simple, it suffers from a major defect known as the "blind proposal" problem: particles $\mathbf{x}_t^{(i)}$ are propagated forward according to the state transition density $f(\mathbf{x}_t \mid \mathbf{x}_{t-1}^{(i)})$ without taking into account the newly arrived observation $\mathbf{y}_t$. If the measurement noise is small, the likelihood function $g(\mathbf{y}_t \mid \mathbf{x}_t)$ forms a narrow peak in state space. Propagating particles blindly causes almost all particles to land in regions of near-zero likelihood, triggering catastrophic particle collapse.

Michael K. Pitt and Neil Shephard (1999) solved this deficiency by inventing the Auxiliary Particle Filter (APF).

1. The Auxiliary Joint Density

The key insight of the APF is to introduce an auxiliary index variable $k \in \{1, \dots, N\}$ that identifies which particle from time $t-1$ will be propagated forward to generate state $\mathbf{x}_t$. The joint target distribution of the pair $(k, \mathbf{x}_t)$ conditional on $\mathbf{y}_{1:t}$ is:

Instead of sampling from this intractable joint density directly, the APF samples from an auxiliary proposal distribution:

2. Pre-Weighting and Lookahead

To construct $q(k \mid \mathbf{y}_{1:t})$, the APF computes a deterministic point estimate $\boldsymbol{\mu}_t^{(k)}$ of the future state given $\mathbf{x}_{t-1}^{(k)}$ (such as the conditional mean $\boldsymbol{\mu}_t^{(k)} = \mathbb{E}[\mathbf{X}_t \mid \mathbf{x}_{t-1}^{(k)}]$ or the mode). The first-stage auxiliary weights are defined as:

The auxiliary weights $\lambda_k$ act as a predictive lookahead: particles whose predicted future states are consistent with the new observation $\mathbf{y}_t$ receive large auxiliary weights, while particles that would produce low likelihood are pruned before state propagation occurs.

3. Two-Stage APF Algorithm

1. Auxiliary Selection: Resample indices $k_1, \dots, k_N$ from $\{1, \dots, N\}$ with probabilities proportional to $\lambda_k$.

2. State Propagation: For each selected index $k_i$, draw the candidate state:

3. Weight Correction: Assign the second-stage importance weights:

When $q(\mathbf{x}_t \mid \mathbf{x}_{t-1}^{(k_i)}) = f(\mathbf{x}_t \mid \mathbf{x}_{t-1}^{(k_i)})$, the weight simplifies to:

Because $\mathbf{x}_t^{(i)}$ is close to $\boldsymbol{\mu}_t^{(k_i)}$, the ratio $\tilde{w}_t^{(i)}$ is nearly constant, keeping the Effective Sample Size close to $N$ and eliminating weight variance even in highly informative measurement regimes.

In evolutionary genomics, complex queueing networks, ecology, and agent-based macroeconomics, the underlying scientific forward simulation model can easily generate synthetic data $\mathbf{x} \sim p(\mathbf{x} \mid \boldsymbol{\theta})$, but the analytical likelihood function:

involves integrating over millions of unobserved latent variables $\mathbf{z}$, rendering the likelihood completely intractable.

Approximate Bayesian Computation (ABC) enables exact or approximate Bayesian inference using exclusively forward stochastic simulations.

1. The Rejection ABC Algorithm

Let $\mathbf{y}_{\text{obs}}$ be the observed dataset. Rejection ABC proceeds as follows:

1. Sample candidate parameter vector $\boldsymbol{\theta}^* \sim p(\boldsymbol{\theta})$ from the prior.
2. Simulate a synthetic dataset $\mathbf{x}^* \sim p(\mathbf{x} \mid \boldsymbol{\theta}^*)$ using the forward simulator.
3. Compute a vector of summary statistics $S(\mathbf{x}^*)$ and compare it to the observed summary statistics $S(\mathbf{y}_{\text{obs}})$ using a distance metric $\rho(\cdot, \cdot)$:

1. Accept $\boldsymbol{\theta}^*$ if $d \le \epsilon$; otherwise, reject.

Theorem (Rubin 1984, Tavare et al. 1997): If $S(\cdot)$ is a vector of sufficient statistics for $\boldsymbol{\theta}$ and the tolerance threshold $\epsilon \to 0$, the accepted draws sample exactly from the true posterior distribution $p(\boldsymbol{\theta} \mid \mathbf{y}_{\text{obs}})$.

2. Sequential Monte Carlo ABC (ABC-SMC)

In practice, setting $\epsilon$ small results in an acceptance rate near zero under standard rejection ABC. Tina Toni, David Welch, Natalja Strelkowa, Andreas Ipsen, and Michael P. H. Stumpf (2009) formulated ABC-SMC, which propagates a population of particles through a decaying sequence of tolerance thresholds $\epsilon_1 > \epsilon_2 > \dots > \epsilon_T$:

1. At population iteration $t = 1$, obtain particles using standard rejection ABC with broad tolerance $\epsilon_1$.
2. For subsequent stages $t = 2, \dots, T$:

• Sample a particle $\boldsymbol{\theta}^*$ from the previous population $t-1$ with weights $w_{t-1}^{(i)}$.
• Perturb the particle using a transition kernel: $\boldsymbol{\theta}^{**} \sim K_t(\boldsymbol{\theta} \mid \boldsymbol{\theta}^*)$.
• Simulate synthetic data $\mathbf{x} \sim p(\mathbf{x} \mid \boldsymbol{\theta}^{**})$.
• If $\rho(S(\mathbf{x}), S(\mathbf{y}_{\text{obs}})) \le \epsilon_t$, accept $\boldsymbol{\theta}^{**}$ and compute its updated importance weight:

ABC-SMC maintains high acceptance rates across orders-of-magnitude reductions in $\epsilon$, providing an industrial-strength engine for simulator calibration.

Over the past decade, the interface between machine learning and computational statistics has generated Simulation-Based Inference (SBI), replacing heuristic summary statistic distances with deep neural density estimators.

Prominent modern SBI paradigms include:

1. Neural Posterior Estimation (NPE): Trains a continuous Normalizing Flow $q_\phi(\boldsymbol{\theta} \mid \mathbf{y})$ (such as Masked Autoregressive Flow or RealNVP) by minimizing the expected Kullback-Leibler divergence over synthetic parameter-data pairs $(\boldsymbol{\theta}_i, \mathbf{x}_i) \sim p(\boldsymbol{\theta}) p(\mathbf{x} \mid \boldsymbol{\theta})$:

Once trained, the normalizing flow amortizes posterior inference: evaluating $q_\phi(\boldsymbol{\theta} \mid \mathbf{y}_{\text{obs}})$ generates millions of independent posterior draws in milliseconds without running any further simulations.

2. Neural Likelihood Estimation (NLE): Trains a neural density estimator $q_\psi(\mathbf{x} \mid \boldsymbol{\theta})$ to approximate the intractable simulator likelihood directly, which is subsequently embedded into standard MCMC or HMC samplers to obtain exact Bayesian posteriors.

When analyzing non-linear state space models, the joint parameter-state inference problem involves two coupled estimation layers:

1. The static parameter vector $\boldsymbol{\theta}$ governing system dynamics (e.g., volatility parameters, transmission rates).
2. The dynamic time-varying latent state trajectory $\mathbf{x}_{1:T}$.

Standard MCMC can alternate between updating $\boldsymbol{\theta} \mid \mathbf{x}_{1:T}$ and updating $\mathbf{x}_{1:T} \mid \boldsymbol{\theta}$. However, because the dimension of $\mathbf{x}_{1:T}$ grows linearly with time $T$, the parameters and states become tightly coupled, causing Gibbs samplers to mix with glacial slowness.

Christophe Andrieu, Arnaud Doucet, and Roman Holenstein (2010) formulated Particle Markov Chain Monte Carlo (PMCMC), embedding Sequential Monte Carlo inside MCMC to sample directly from the exact marginal posterior $p(\boldsymbol{\theta} \mid \mathbf{y}_{1:T})$.

1. The Particle Marginal Metropolis-Hastings (PMMH) Algorithm

Suppose we propose a new parameter vector $\boldsymbol{\theta}^* \sim q(\boldsymbol{\theta}^* \mid \boldsymbol{\theta})$. To evaluate the ideal marginal Metropolis-Hastings ratio:

the marginal likelihood $p(\mathbf{y}_{1:T} \mid \boldsymbol{\theta}) = \prod_{t=1}^T p(\mathbf{y}_t \mid \mathbf{y}_{1:t-1}, \boldsymbol{\theta})$ must be evaluated.

Because this likelihood is analytically intractable, PMMH runs a bootstrap particle filter with $N$ particles targeted at $\boldsymbol{\theta}^*$. The particle filter provides an unbiased Monte Carlo estimate of the marginal likelihood:

PMMH replaces the exact intractable likelihood with this particle-filter estimate in the acceptance ratio:

2. The Pseudo-Marginal Theorem

A naive statistician might assume that substituting a stochastic estimate $\hat{p}(\mathbf{y}_{1:T} \mid \boldsymbol{\theta})$ introduces approximation error into the resulting Markov chain. Remarkably, the Pseudo-Marginal Theorem (Beaumont 2003, Andrieu and Roberts 2009) proves that this is not the case:

Theorem (Pseudo-Marginal Validity): If the likelihood estimator $\hat{p}(\mathbf{y}_{1:T} \mid \boldsymbol{\theta})$ is non-negative and strictly unbiased:

for all $\boldsymbol{\theta}$, then the stationary distribution of the Markov chain generated by PMMH is identically and exactly the true posterior distribution $p(\boldsymbol{\theta} \mid \mathbf{y}_{1:T})$, regardless of the number of particles $N \ge 1$.

Proof Intuition: The particle filter's internal random numbers $\mathbf{U}$ can be viewed as auxiliary variables on an extended state space $(\boldsymbol{\theta}, \mathbf{U})$. Because the joint target distribution on this extended space integrates out to the exact marginal target, the Metropolis-Hastings acceptance ratio on the extended space matches $\alpha_{\text{PMMH}}$ exactly.

The number of particles $N$ affects only the variance of the likelihood estimator: if $N$ is too small, $\hat{p}$ occasionally overestimates the likelihood, trapping the chain at an unrepresentative parameter value. Doucet et al. (2015) proved that the optimal computational efficiency is achieved when $N$ is tuned so that the standard deviation of $\log \hat{p}(\mathbf{y}_{1:T} \mid \boldsymbol{\theta})$ is approximately $1.2$ to $1.8$.

In numerical weather prediction, oceanography, and atmospheric monitoring, dynamic state vectors represent temperature, atmospheric pressure, and wind velocities across three-dimensional spatial grids, with state dimensions reaching $d \approx 10^8$ to $10^9$. In this ultra-high-dimensional regime, standard particle filters collapse instantly due to particle degeneracy, while classical Kalman filters cannot store or invert the $10^9 \times 10^9$ error covariance matrix.

Geir Evensen (1994) introduced the Ensemble Kalman Filter (EnKF), combining Monte Carlo ensemble forecasting with Kalman update geometry to make global data assimilation computationally feasible.

1. The Stochastic Forecast and Analysis Cycle

The EnKF represents the probability distribution of the high-dimensional state $\mathbf{x}_t \in \mathbb{R}^d$ through an ensemble of $N$ model states (typically $N \approx 50$ to $100$):

1. Forecast Step: Each ensemble member is integrated forward in time through the non-linear forward geophysical simulation model $\mathcal{M}$:

The ensemble mean $\bar{\mathbf{x}}$ and ensemble perturbation matrix $\mathbf{A} = \mathbf{X} - \bar{\mathbf{x}} \mathbf{1}^T$ define an implicit low-rank sample covariance matrix:

Crucially, $\mathbf{P}^f$ is never explicitly formed in memory.

2. Perturbed Observation Analysis Step: Let $\mathbf{y}_t \in \mathbb{R}^p$ be the newly arrived satellite observation vector with measurement error covariance $\mathbf{R}$. To preserve proper ensemble spread, independent random perturbations are added to the observation vector for each member (Burgers et al. 1998):

The analysis update for member $i$ is:

where the Kalman gain matrix $\mathbf{K}$ is computed directly in the low-dimensional observation space:

2. Covariance Localization and Inflation

Because the ensemble size $N \approx 100$ is astronomically smaller than the state dimension $d \approx 10^8$, the sample covariance $\mathbf{P}^f$ suffers from two critical sampling artifacts:

1. Spurious Long-Range Correlations: Small random correlations appear between physically distant geographic locations (e.g., surface temperature in Paris and rainfall in the South Pacific). Covariance Localization applies a Schur (element-wise) product with a compact smooth spatial correlation function $\boldsymbol{\rho}$:

forcing correlations between distant grid points to zero.

2. Filter Divergence (Inbreeding): Because the ensemble systematically underestimates forecast error variance, the filter progressively ignores new incoming observations. Covariance Inflation multiplies the perturbation matrix by an inflation factor $\lambda > 1$:

restoring realistic ensemble spread and ensuring stable, continuous tracking of planetary weather dynamics.

While filtering estimates the marginal filtering distribution $p(\mathbf{x}_t \mid \mathbf{y}_{1:t})$ online as observations arrive, retrospective historical reconstruction requires estimating the full joint smoothing distribution:

Naively extracting the ancestral path trajectories of surviving particles from a forward particle filter results in severe Path Degeneracy: due to repeated resampling across $T$ steps, all surviving particles at time $T$ trace back to a single common ancestor at time $t \ll T$, collapsing historical trajectory diversity.

Forward-Filtering Backward-Smoothing (FFBS)

Doucet, Godsill, and Andrieu (2000) resolved path degeneracy with the Forward-Filtering Backward-Smoothing (FFBS) algorithm:

1. Run a standard forward particle filter, storing all particles $\{\mathbf{x}_t^{(i)}\}$ and normalized weights $\{w_t^{(i)}\}$ for all timesteps $t = 1, \dots, T$.
2. At the final timestep $T$, sample an index $k_T$ with probability $w_T^{(k_T)}$ and set state $\tilde{\mathbf{x}}_T = \mathbf{x}_T^{(k_T)}$.
3. For $t = T-1, T-2, \dots, 0$, compute backward transition weights:

1. Sample state index $k_t$ with probability $w_{t|t+1}^{(k_t)}$ and set $\tilde{\mathbf{x}}_t = \mathbf{x}_t^{(k_t)}$.

By decoupling historical smoothing from the forward resampling bottleneck, FFBS produces independent, non-degenerate historical state trajectories across the entire time horizon $[0, T]$.

In robotics and autonomous navigation, a vehicle must estimate its two-dimensional position $(x, y)$ and heading angle $\theta$ inside an urban canyon where GPS signals suffer from severe multipath interference. The vehicle relies on wheel odometry sensors (noisy kinematic drift) and range-bearing LIDAR landmark scans.

Comparing localization filters across a 5-kilometer urban route:

• Extended Kalman Filter (EKF): Linearizes measurement equations via Taylor expansions. When the vehicle executes sharp turns with high bearing non-linearities, the EKF covariance matrix diverges, causing the estimated path to drift by over 45 meters.
• Sequential Importance Resampling (SIR) Filter ($N = 2,000$ particles): Accurately tracks multimodal position hypotheses during GPS dropouts, maintaining a localization error within $\pm 0.85$ meters.
• Auxiliary Particle Filter (APF): Eliminates particle depletion during high-speed maneuvers, sustaining continuous sub-meter tracking with only $N = 500$ particles and reducing CPU utilization by $68\%$.

Unit 8 synthesized sequential simulation, state space tracking, and modern simulation-based inference. We formulated state-space hidden Markov models and derived the Bayesian recursive filtering equations. We detailed Sequential Importance Sampling (SIS) and demonstrated how weight degeneration inevitably collapses particle diversity. We explored Sequential Importance Resampling (SIR) under multinomial, stratified, residual, and systematic resampling schemes, evaluating the Effective Sample Size $N_{\text{eff}} = 1 / \sum (w_t^{(i)})^2$. We examined Pitt and Shephard's Auxiliary Particle Filter (APF), showing how predictive lookahead auxiliary weights prevent particle starvation during informative observation updates. We analyzed Approximate Bayesian Computation (ABC) rejection and ABC-SMC algorithms for likelihood-free scientific models, and surveyed the Pseudo-Marginal theorem underpinning Particle MCMC (PMMH). Finally, we examined the Ensemble Kalman Filter (EnKF) for planetary-scale data assimilation and Neural Posterior Estimation (NPE) via normalizing flows.

In contemporary computational statistics, Normalizing Flows provide an expressive class of invertible generative models that bridge parametric density estimation and Monte Carlo simulation.

A normalizing flow transforms a simple base distribution (such as a standard isotropic Gaussian $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$) into a complex, multimodal target density $\mathbf{x} \in \mathbb{R}^d$ through a sequence of $K$ invertible, differentiable transformations:

By the multivariate change-of-variables theorem, the exact target log-density is computable in closed form:

By designing transformations with triangular Jacobian matrices (such as RealNVP or Masked Autoregressive Flows), evaluating the determinant requires only $O(d)$ operations. Once trained on simulation trajectories, normalizing flows generate millions of exact synthetic data draws in milliseconds, amortizing inference across complex physical systems.