1 Motivation

A particle filter, also called a Sequential Monte Carlo method, is a simulation-based filtering algorithm for learning the hidden state of a dynamic system from noisy, partial, or nonlinear measurements.

Scientific picture. A planet, robot, animal, fluid tracer, molecule, cyclone, or epidemic state evolves over time. We do not observe the true state directly. We observe noisy sensor data. A filter updates our best uncertainty statement about the present state as data arrive sequentially.

Think of the algorithm as a probabilistic version of Marco Polo:

  • the hidden target is the true state;
  • the noisy voice is the measurement;
  • the laws of motion are the physical model;
  • the particles are many simulated guesses;
  • the posterior cloud shows what is still plausible after seeing the data.

2 The state-space model

Let

\[ X_t \in \mathcal X \]

denote the hidden state at time \(t\), and let

\[ Z_t \in \mathcal Z \]

denote the observation at time \(t\). A state-space model has two probability laws.

\[ \underbrace{p(x_t\mid x_{t-1})}_{\text{transition / physical model}}, \qquad \underbrace{p(z_t\mid x_t)}_{\text{measurement / sensor model}}. \]

The filtering goal is not merely to estimate one number. The goal is to compute or approximate the filtering distribution

\[ p(x_t\mid z_{1:t}), \qquad z_{1:t}=(z_1,\ldots,z_t). \]

This distribution expresses our uncertainty about the present hidden state after all observations up to time \(t\).

Important distinction. A point estimate answers: “Where is the target most likely?” A filtering distribution answers: “Which states are plausible, and with what uncertainty?”

3 Kalman filter versus particle filter

The classical Kalman filter is elegant and exact under the linear-Gaussian state-space model:

\[ X_t = A X_{t-1}+\varepsilon_t,\qquad \varepsilon_t\sim N(0,Q), \]

\[ Z_t = H X_t+\eta_t,\qquad \eta_t\sim N(0,R). \]

Under these assumptions, the filtering distribution remains Gaussian, and only its mean and covariance need to be updated.

A particle filter is more general. It approximates

\[ p(x_t\mid z_{1:t}) \]

by a weighted empirical distribution:

\[ \widehat p_N(dx_t\mid z_{1:t}) = \sum_{i=1}^N w_t^{(i)}\delta_{x_t^{(i)}}(dx_t), \]

where \(x_t^{(i)}\) is the \(i\)-th particle and \(w_t^{(i)}\) is its weight.

4 The three-step algorithm

Step Question Mathematical_object
Predict Where could the system move next under the physical model? Draw x_t^(i) from p(x_t | x_{t-1}^(i))
Update Which particles are more consistent with the new measurement? Set w_t^(i) proportional to p(z_t | x_t^(i))
Resample How do we remove bad guesses and duplicate good guesses? Sample new particles with probability proportional to w_t^(i)

The effective sample size is often monitored:

\[ \mathrm{ESS}_t=\frac{1}{\sum_{i=1}^N (w_t^{(i)})^2}. \]

A small ESS means that only a few particles dominate; resampling is then needed.

5 Interactive demonstration

Use the controls below. The red dot is the true target. The blue cross is the noisy sensor reading. The grey dots are particles. The green cross is the filter estimate.

Interactive Particle Filter Demonstration 

<span style="color: red;">● True Target</span>
<span style="color: blue;">✖ Sensor Reading</span>
<span style="color: gray;">• Particles</span>
<span style="color: green;">✚ Estimated State</span>
<span style="color: purple;">— True Path</span>
Status will appear here. 
<div style="text-align:center;">
  <button id="stepBtn">Start Demo: PREDICT</button>
  <button id="autoBtn">Auto Run</button>
  <button id="resetBtn">Reset</button>
</div>
<div class="slider-row">
  <label>Particles: <input type="range" id="pCount" min="50" max="1500" value="300" step="50"></label>
  <span id="pCountVal">300</span>
</div>
<div class="slider-row">
  <label>Sensor Noise: <input type="range" id="noiseLevel" min="5" max="90" value="30"></label>
  <span id="noiseVal">30</span>
</div>
<div class="slider-row">
  <label>Process Noise: <input type="range" id="processNoise" min="1" max="45" value="12"></label>
  <span id="processVal">12</span>
</div>

6 What the animation teaches

  • Predict: the particles move according to the physical model. Uncertainty spreads.
  • Update: the sensor reading arrives. Particles near the observation receive large weights.
  • Resample: low-weight guesses disappear; high-weight guesses are copied. The cloud concentrates around plausible states.

Why this is Bayesian. The algorithm recursively approximates the posterior distribution \(p(x_t\mid z_{1:t})\). It is not merely a tracking trick; it is sequential probabilistic inference.

7 A rigorous R simulation: nonlinear filtering

Now we implement a one-dimensional nonlinear state-space model in R. This example is famous because the observation only sees a nonlinear transformation of the state. Thus a naive method can be misled.

\[ X_t = 0.5X_{t-1}+\frac{25X_{t-1}}{1+X_{t-1}^2}+8\cos(1.2t)+\varepsilon_t, \qquad \varepsilon_t\sim N(0,q). \]

\[ Y_t=\frac{X_t^2}{20}+\eta_t, \qquad \eta_t\sim N(0,r). \]

The observation \(Y_t\) depends on \(X_t^2\), so it cannot easily distinguish positive and negative states. This creates a non-Gaussian and sometimes multi-modal filtering problem.

Tn <- 70
q <- 10
r <- 1
x_true <- numeric(Tn)
y_obs <- numeric(Tn)
x_true[1] <- rnorm(1, 0, sqrt(5))
y_obs[1] <- x_true[1]^2 / 20 + rnorm(1, 0, sqrt(r))

for (t in 2:Tn) {
  x_true[t] <- 0.5 * x_true[t - 1] +
    25 * x_true[t - 1] / (1 + x_true[t - 1]^2) +
    8 * cos(1.2 * t) + rnorm(1, 0, sqrt(q))
  y_obs[t] <- x_true[t]^2 / 20 + rnorm(1, 0, sqrt(r))
}

df_sim <- data.frame(t = 1:Tn, x_true = x_true, y_obs = y_obs)
head(df_sim)
##   t     x_true      y_obs
## 1 1  -3.689294  0.3511910
## 2 2 -14.359539 11.3957252
## 3 3 -15.905507 11.2443604
## 4 4  -9.223936  2.9993445
## 5 5  -1.843543  0.7124554
## 6 6  -1.642937  0.2680033

The observed data are not the hidden state; they are noisy nonlinear measurements of it.

8 Particle filter from scratch in R

weighted_quantile <- function(x, w, probs = c(0.025, 0.5, 0.975)) {
  ord <- order(x)
  x <- x[ord]
  w <- w[ord] / sum(w)
  cw <- cumsum(w)
  sapply(probs, function(p) x[which(cw >= p)[1]])
}

systematic_resample <- function(w) {
  N <- length(w)
  u0 <- runif(1, 0, 1 / N)
  u <- u0 + (0:(N - 1)) / N
  cw <- cumsum(w)
  idx <- integer(N)
  j <- 1
  for (i in seq_len(N)) {
    while (u[i] > cw[j]) j <- j + 1
    idx[i] <- j
  }
  idx
}

pf_nonlinear <- function(y, N = 3000, q = 10, r = 1) {
  Tn <- length(y)
  particles <- rnorm(N, 0, sqrt(5))
  weights <- rep(1 / N, N)
  out <- data.frame(t = 1:Tn, mean = NA_real_, lo = NA_real_, med = NA_real_, hi = NA_real_, ess = NA_real_)

  for (t in seq_len(Tn)) {
    if (t > 1) {
      particles <- 0.5 * particles +
        25 * particles / (1 + particles^2) +
        8 * cos(1.2 * t) + rnorm(N, 0, sqrt(q))
    }

    logw <- dnorm(y[t], mean = particles^2 / 20, sd = sqrt(r), log = TRUE)
    logw <- logw - max(logw)
    weights <- exp(logw)
    weights <- weights / sum(weights)

    qs <- weighted_quantile(particles, weights)
    out$mean[t] <- sum(weights * particles)
    out$lo[t] <- qs[1]
    out$med[t] <- qs[2]
    out$hi[t] <- qs[3]
    out$ess[t] <- 1 / sum(weights^2)

    if (out$ess[t] < N / 2) {
      idx <- systematic_resample(weights)
      particles <- particles[idx]
      weights <- rep(1 / N, N)
    }
  }
  out
}

pf_out <- pf_nonlinear(y_obs, N = 3000, q = q, r = r)
pf_out$x_true <- x_true
pf_out$abs_sensor <- sqrt(pmax(20 * y_obs, 0))

Quantity Value
RMSE: particle filter mean 3.126
RMSE: naive sqrt(20Y) inversion 17.236
Empirical 95% interval coverage 0.957

9 What can go wrong?

Particle filters are powerful, but not magical.

Method RMSE Coverage
Well-tuned particle filter 3.126 0.957
Overconfident misspecified filter 6.311 0.486

Statistical lesson. Filtering is not just computation. It requires a scientifically meaningful transition model, a realistic measurement model, and uncertainty checks. A filter with wrong noise assumptions can be more dangerous than no filter, because it may look precise while being wrong.

10 Where basic-science researchers use this idea

Domain Hidden_state Observation Why_filtering
Astronomy object position/brightness telescope image / photon counts tracking under noise
Geoscience seismic or groundwater state sensor and field measurements sequential hazard updating
Soft matter particle/tracer state microscopy frames stochastic motion inference
Biophysics molecular/ion-channel state noisy electrophysiology hidden-state dynamics
Ecology animal location/behaviour GPS tags / acoustic detections movement ecology
Epidemiology latent infection level case reports / tests real-time nowcasting

11 Practical advice for students

  1. Start with a scientific state variable: what is hidden but important?
  2. Write a transition model: how should the state evolve?
  3. Write an observation model: how does the instrument see the state?
  4. Simulate before using real data.
  5. Check uncertainty, not only the point estimate.
  6. Test sensitivity to noise parameters and number of particles.

12 Concluding message

Particle filtering is a beautiful example of Bayesian data science in motion. It combines probability, simulation, scientific modelling, and computation. It is especially useful when the system is dynamic, partially observed, nonlinear, and noisy.

Final message. Machine learning may learn patterns from data, but a filtering model asks a scientific question: how does a hidden state evolve, and what do the measurements tell us about it now?

13 References for further reading

  • Doucet, A., de Freitas, N., and Gordon, N. (eds.). Sequential Monte Carlo Methods in Practice. Springer, 2001.
  • Gordon, N. J., Salmond, D. J., and Smith, A. F. M. (1993). Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings F.
  • Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering.
  • Särkkä, S. (2013). Bayesian Filtering and Smoothing. Cambridge University Press.
  • Arulampalam, M. S., Maskell, S., Gordon, N., and Clapp, T. (2002). A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing.