Particle Filter

Particle filter, formally known in statistics and signal processing as a Sequential Monte Carlo (SMC) method, is a simulation-based recursive Bayesian estimation algorithm in computation and algorithms that tracks the hidden state x_t of a dynamical system over time from noisy observations y_1:t = {y₁, y₂, …, y_t}. Whereas the linear Kalman filter computes exact closed-form Gaussian posteriors when both the state transition and observation models are linear with additive Gaussian noise, real-world physical, biological, and robotic systems frequently exhibit strong nonlinearities, multimodal likelihood landscapes, or non-Gaussian heavy-tailed noise where analytical integrals are intractable.

Introduced in its modern bootstrap form by Neil J. Gordon, David J. Salmond, and Adrian F. M. Smith in 1993 by combining sequential importance sampling with a Darwinian resampling step, the particle filter approximates the continuous filtering posterior density p(x_t | y_1:t) by a discrete empirical cloud of N weighted random samples—called particles—{(x_t⁽ⁱ⁾, w_t⁽ⁱ⁾)}_i=1^N. As a dynamic extension of the Monte Carlo method, expectations of any state function φ(x_t) converge to the true Bayesian posterior expectation at the dimension-independent Monte Carlo rate O(N^−1/2).

A discrete-time Markov state-space model (hidden Markov model with continuous state space ℝ^d) is specified by an initial state prior p(x₀) and two conditional probability distributions:

1. State transition (process) model: x_t ∼ p(x_t | x_t−1), governing the stochastic dynamics of the hidden state.
2. Observation (measurement) likelihood: y_t ∼ p(y_t | x_t), relating the noisy sensor reading y_t to the current hidden state x_t.

Recursive Bayesian filtering alternates between two conceptual steps at each time step t:

• Prediction (Chapman–Kolmogorov equation): p(x_t | y_1:t−1) = ∫ p(x_t | x_t−1) p(x_t−1 | y_1:t−1) dx_t−1
• Update (Bayes' rule): p(x_t | y_1:t) = [p(y_t | x_t) p(x_t | y_1:t−1)] / p(y_t | y_1:t−1).

When p(x_t | x_t−1) or p(y_t | x_t) is nonlinear or bimodal, the integral in the Chapman–Kolmogorov equation cannot be evaluated in closed form. A particle filter replaces the continuous prior density with the empirical measure:

p̂(x_t | y_1:t) = ∑_i=1^N w_t⁽ⁱ⁾ δ(x_t − x_t⁽ⁱ⁾), where ∑_i=1^N w_t⁽ⁱ⁾ = 1 and w_t⁽ⁱ⁾ ≥ 0,

and δ(·) denotes the Dirac delta measure.

Sequential Importance Sampling (SIS)

In Sequential Importance Sampling (SIS), each particle trajectory x_0:t⁽ⁱ⁾ is extended at step t by drawing a new state x_t⁽ⁱ⁾ ∼ q(x_t | x_t−1⁽ⁱ⁾, y_t) from a user-chosen proposal (importance) distribution q. The unnormalized importance weight is updated recursively by:

w̃_t⁽ⁱ⁾ = w_t−1⁽ⁱ⁾ · [p(y_t | x_t⁽ⁱ⁾) p(x_t⁽ⁱ⁾ | x_t−1⁽ⁱ⁾)] / q(x_t⁽ⁱ⁾ | x_t−1⁽ⁱ⁾, y_t),

followed by normalization w_t⁽ⁱ⁾ = w̃_t⁽ⁱ⁾ / ∑_j=1^N w̃_t^(j). In the widely used Bootstrap Particle Filter (Gordon et al., 1993), the proposal distribution is chosen directly as the transition prior q(x_t | x_t−1⁽ⁱ⁾, y_t) = p(x_t | x_t−1⁽ⁱ⁾), simplifying the weight update to the observation likelihood alone: w̃_t⁽ⁱ⁾ = w_t−1⁽ⁱ⁾ p(y_t | x_t⁽ⁱ⁾).

Weight Degeneracy and the SIR Resampling Step

Without intervention, pure SIS suffers from an unavoidable flaw: the variance of the importance weights increases stochastically with time t. Within a few dozen time steps, one particle accumulates normalized weight w_t^(k) ≈ 1 while the remaining N − 1 particles decay to w_t⁽ⁱ⁾ ≈ 0—a phenomenon known as weight degeneracy. The degree of degeneracy is quantified by Kong's Effective Sample Size (ESS):

N_eff ≈ 1 / ∑_i=1^N (w_t⁽ⁱ⁾)², where 1 ≤ N_eff ≤ N.

Whenever N_eff drops below a threshold (typically N/2), the Sequential Importance Resampling (SIR) algorithm draws N new particles with replacement from the discrete distribution {(x_t⁽ⁱ⁾, w_t⁽ⁱ⁾)} and resets all weights to uniform w_t⁽ⁱ⁾ = 1/N. Low-variance systematic resampling and stratified resampling perform this selection in O(N) time using a single uniform random number u₀ ∼ Uniform(0, 1/N], eliminating particles in low-probability regions and duplicating particles in high-likelihood regions.

Repeated resampling introduces a secondary challenge called sample impoverishment (or lineage collapse): if process noise is small, resampling duplicates identical particles so that after k steps most particles share a single common ancestor at time t − k. Advanced SMC variants mitigate impoverishment and high-dimensional weight collapse through three strategies:

1. Auxiliary Particle Filter (Pitt & Shephard, 1999): Looks one step ahead using the current observation y_t to weight particles before propagating and resampling them, steering particles toward high-likelihood regions.
2. Resample-Move and Particle Markov Chain Monte Carlo (PMCMC): Applies a Markov chain Monte Carlo (Metropolis–Hastings) jitter kernel after resampling to restore particle diversity, or uses the unbiased marginal likelihood estimator p̂(y_1:T) = ∏_t=1^T [(1/N) ∑_i w̃_t⁽ⁱ⁾] inside Particle Marginal Metropolis–Hastings (PMMH) (Andrieu, Doucet, & Holenstein, 2010) to perform exact Bayesian parameter inference in nonlinear stochastic differential equations and epidemiological SIR models.
3. Rao–Blackwellized Particle Filtering (FastSLAM): In mobile robotics Simultaneous Localization and Mapping (SLAM), the robot's 3-DoF trajectory x_1:t is tracked by N particles while thousands of landmark coordinates conditioned on each trajectory are updated analytically by independent 2 × 2 Kalman filters, dramatically reducing the variance of the Monte Carlo estimator.

Why is systematic resampling preferred over multinomial resampling?

Multinomial resampling draws N independent uniform random numbers and searches the cumulative weight array in O(N log N) time, introducing unnecessary Monte Carlo sampling variance even when all weights are equal (w_t⁽ⁱ⁾ = 1/N). Systematic resampling lays down an evenly spaced comb of N points with spacing 1/N in O(N) time, achieving zero sampling variance when weights are uniform and minimizing variance across all unbiased resampling schemes.

How does the particle filter estimate the marginal likelihood of a model?

At each time step t, prior to resampling, the average of the unnormalized weights Ẑ_t = (1/N) ∑_i=1^N w̃_t⁽ⁱ⁾ estimates the predictive likelihood p(y_t | y_1:t−1). By Del Moral's theorem, the product Ẑ_1:T = ∏_t=1^T Ẑ_t is an unbiased estimator of the total marginal likelihood p(y_1:T) for any finite particle count N ≥ 1.

1. Gordon, N. J., Salmond, D. J., & Smith, A. F. M. (1993). Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings F (Radar and Signal Processing), 140(2), 107–113. https://doi.org/10.1049/ip-f-2.1993.0015
2. Pitt, M. K., & Shephard, N. (1999). Filtering via simulation: Auxiliary particle filters. Journal of the American Statistical Association, 94(446), 590–599. https://doi.org/10.1080/01621459.1999.10474153
3. Doucet, A., Godsill, S., & Andrieu, C. (2000). On sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and Computing, 10(3), 197–208. https://doi.org/10.1023/A:1008935410038
4. Andrieu, C., Doucet, A., & Holenstein, R. (2010). Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3), 269–342. https://doi.org/10.1111/j.1467-9868.2009.00736.x