1 Why this page exists

This RPubs page is designed for students and researchers in basic sciences who want a visual and computational introduction to filtering, especially particle filtering.

The scientific setting is common:

Typical examples:

Domain Hidden state Measurement Why filtering is needed
Astronomy true position, brightness, orbit, source state noisy telescope image or flux sequential tracking under measurement error
Geoscience river flow regime, seismic stress, subsurface contamination sensor readings, catalogues, sparse surveys hidden process evolves through time
Soft matter / materials latent microstructural state microscopy, spectroscopy, scattering observation is indirect and noisy
Biophysics cell state, voltage, molecular conformation fluorescence, patch-clamp, imaging signal is stochastic and partially observed
Epidemiology true infection process reported cases reporting is noisy and delayed

The main lesson is:

Filtering is not merely smoothing a curve. It is sequential probabilistic inference about a hidden state.

2 State-space model: the common language

Let

\[ X_t = \text{hidden state at time }t, \qquad Z_t = \text{observed measurement at time }t. \]

A state-space model has two parts:

\[ \textbf{State transition:}\qquad X_t \sim p(x_t\mid x_{t-1}), \]

\[ \textbf{Observation model:}\qquad Z_t \sim p(z_t\mid x_t). \]

The filtering target is the posterior distribution

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

where

\[ z_{1:t}=(z_1,z_2,\ldots,z_t). \]

In words: after seeing all measurements up to time \(t\), what do we believe about the current hidden state?

3 Kalman filter vs particle filter

The Kalman filter is elegant and exact when the system is linear and Gaussian:

\[ X_t = aX_{t-1}+\eta_t, \qquad Z_t = bX_t+\epsilon_t, \]

where \(\eta_t\) and \(\epsilon_t\) are Gaussian noises.

The particle filter is more general. It represents the filtering posterior by many simulated particles:

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

where:

The effective sample size is

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

If ESS is small, only a few particles are carrying most of the probability mass; resampling is needed.

alg <- tibble::tribble(
  ~Step, ~Question, ~Mathematical_object, ~Scientific_meaning,
  "Predict", "Where could the system move next?", "simulate x_t^(i) from p(x_t | x_{t-1}^(i))", "propagate physics/dynamics",
  "Update", "Which particles agree with the measurement?", "w_t^(i) proportional to p(z_t | x_t^(i))", "score particles by sensor likelihood",
  "Resample", "How do we remove bad guesses?", "sample particles with probability w_t^(i)", "keep plausible hypotheses alive"
)
knitr::kable(alg)
Step Question Mathematical_object Scientific_meaning
Predict Where could the system move next? simulate x_t^(i) from p(x_t | x_{t-1}^(i)) propagate physics/dynamics
Update Which particles agree with the measurement? w_t^(i) proportional to p(z_t | x_t^(i)) score particles by sensor likelihood
Resample How do we remove bad guesses? sample particles with probability w_t^(i) keep plausible hypotheses alive

4 Why the previous raw HTML animation may fail on RPubs

Some hosted platforms sanitize or escape raw JavaScript. If raw <script>, <canvas>, <button>, or custom HTML is escaped, then the page shows HTML tags as text rather than running the animation.

For RPubs, the safer route is to use an R htmlwidget, such as plotly, because it is knitted and published through the standard R Markdown widget pipeline.

So below we create an animation completely from R.

5 Interactive particle-filter animation in R using plotly

This animation shows a target moving in two dimensions. The true state is hidden. A noisy sensor gives measurements. Particles represent possible locations of the target.

simulate_pf2d <- function(T = 35, N = 350, sensor_sd = 25, process_sd = 9) {
  target <- c(x = 80, y = 70, vx = 5.5, vy = 3.2)
  width <- 760
  height <- 430
  particles <- tibble::tibble(
    id = seq_len(N),
    x = runif(N, 0, width),
    y = runif(N, 0, height),
    w = rep(1 / N, N)
  )

  all_frames <- list()
  true_path <- tibble::tibble(x = numeric(), y = numeric())
  est_path <- tibble::tibble(x = numeric(), y = numeric())

  frame_no <- 1

  systematic_resample <- function(w) {
    # Robust systematic resampling.
    # This protects against tiny numerical errors in cumulative weights
    # and against accidental NA/Inf weights after likelihood underflow.
    N <- length(w)
    w[!is.finite(w)] <- 0
    sw <- sum(w)
    if (!is.finite(sw) || sw <= 0) {
      w <- rep(1 / N, N)
    } else {
      w <- w / sw
    }

    cw <- cumsum(w)
    cw <- cummax(cw)  # enforce non-decreasing cumulative probabilities
    cw[N] <- 1        # ensure the last interval closes exactly at 1

    u0 <- runif(1, 0, 1 / N)
    u <- u0 + (0:(N - 1)) / N
    # Avoid findInterval() failures on hosted knitting environments.
    idx <- integer(N)
    j <- 1L
    for (i in seq_len(N)) {
      while (u[i] > cw[j] && j < N) {
        j <- j + 1L
      }
      idx[i] <- j
    }
    idx
  }

  for (tt in seq_len(T)) {
    # true motion
    target["x"] <- target["x"] + target["vx"]
    target["y"] <- target["y"] + target["vy"]
    if (target["x"] < 0 || target["x"] > width) target["vx"] <- -target["vx"]
    if (target["y"] < 0 || target["y"] > height) target["vy"] <- -target["vy"]

    z <- c(
      x = target["x"] + rnorm(1, 0, sensor_sd),
      y = target["y"] + rnorm(1, 0, sensor_sd)
    )

    # PREDICT
    particles <- particles |>
      dplyr::mutate(
        x = pmin(pmax(x + target["vx"] + rnorm(N, 0, process_sd), 0), width),
        y = pmin(pmax(y + target["vy"] + rnorm(N, 0, process_sd), 0), height),
        phase = "1 Predict"
      )

    pred_est <- particles |>
      dplyr::summarise(x = mean(x), y = mean(y))

    true_path <- dplyr::bind_rows(true_path, tibble::tibble(x = target["x"], y = target["y"]))
    est_path <- dplyr::bind_rows(est_path, pred_est)

    all_frames[[frame_no]] <- dplyr::bind_rows(
      particles |> dplyr::transmute(frame = frame_no, phase, kind = "particles", x, y, size = 3, alpha = 0.35, text = paste0("particle ", id)),
      tibble::tibble(frame = frame_no, phase = "1 Predict", kind = "true target", x = target["x"], y = target["y"], size = 12, alpha = 1, text = "true target"),
      tibble::tibble(frame = frame_no, phase = "1 Predict", kind = "sensor", x = z["x"], y = z["y"], size = 10, alpha = 1, text = "noisy sensor"),
      tibble::tibble(frame = frame_no, phase = "1 Predict", kind = "estimate", x = pred_est$x, y = pred_est$y, size = 11, alpha = 1, text = "particle mean estimate")
    )
    frame_no <- frame_no + 1

    # UPDATE
    dist2 <- (particles$x - z["x"])^2 + (particles$y - z["y"])^2
    lw <- -dist2 / (2 * sensor_sd^2)
    lw <- lw - max(lw)
    raw_w <- exp(lw)
    sw <- sum(raw_w)
    if (!is.finite(sw) || sw <= 0) {
      particles$w <- rep(1 / N, N)
    } else {
      particles$w <- raw_w / sw
    }
    ess <- 1 / sum(particles$w^2)

    upd_est <- particles |>
      dplyr::summarise(x = sum(w * x), y = sum(w * y))

    all_frames[[frame_no]] <- dplyr::bind_rows(
      particles |>
        dplyr::mutate(weight_size = 2 + 12 * w / max(w)) |>
        dplyr::transmute(frame = frame_no, phase = paste0("2 Update: ESS = ", round(ess, 1)), kind = "particles", x, y, size = weight_size, alpha = 0.55, text = paste0("weight = ", signif(w, 3))),
      tibble::tibble(frame = frame_no, phase = paste0("2 Update: ESS = ", round(ess, 1)), kind = "true target", x = target["x"], y = target["y"], size = 12, alpha = 1, text = "true target"),
      tibble::tibble(frame = frame_no, phase = paste0("2 Update: ESS = ", round(ess, 1)), kind = "sensor", x = z["x"], y = z["y"], size = 10, alpha = 1, text = "noisy sensor"),
      tibble::tibble(frame = frame_no, phase = paste0("2 Update: ESS = ", round(ess, 1)), kind = "estimate", x = upd_est$x, y = upd_est$y, size = 11, alpha = 1, text = "weighted estimate")
    )
    frame_no <- frame_no + 1

    # RESAMPLE
    idx <- systematic_resample(particles$w)
    particles <- particles[idx, ] |>
      dplyr::mutate(id = seq_len(N), w = 1 / N)

    res_est <- particles |>
      dplyr::summarise(x = mean(x), y = mean(y))

    all_frames[[frame_no]] <- dplyr::bind_rows(
      particles |> dplyr::transmute(frame = frame_no, phase = "3 Resample", kind = "particles", x, y, size = 3, alpha = 0.35, text = paste0("resampled particle ", id)),
      tibble::tibble(frame = frame_no, phase = "3 Resample", kind = "true target", x = target["x"], y = target["y"], size = 12, alpha = 1, text = "true target"),
      tibble::tibble(frame = frame_no, phase = "3 Resample", kind = "sensor", x = z["x"], y = z["y"], size = 10, alpha = 1, text = "noisy sensor"),
      tibble::tibble(frame = frame_no, phase = "3 Resample", kind = "estimate", x = res_est$x, y = res_est$y, size = 11, alpha = 1, text = "resampled estimate")
    )
    frame_no <- frame_no + 1
  }

  dplyr::bind_rows(all_frames) |>
    dplyr::mutate(
      kind = factor(kind, levels = c("particles", "sensor", "estimate", "true target")),
      frame_label = paste0("Frame ", frame, ": ", phase)
    )
}

pf_anim <- simulate_pf2d(T = 32, N = 320, sensor_sd = 28, process_sd = 9)
plotly::plot_ly(
  data = pf_anim,
  x = ~x,
  y = ~y,
  frame = ~frame_label,
  type = "scatter",
  mode = "markers",
  color = ~kind,
  colors = c("particles" = "#777777",
             "sensor" = "blue",
             "estimate" = "green",
             "true target" = "red"),
  size = ~size,
  sizes = c(3, 16),
  text = ~text,
  hoverinfo = "text",
  marker = list(line = list(width = 0.5, color = "white"))
) |>
  plotly::layout(
    title = "Particle Filter: Predict → Update → Resample",
    xaxis = list(title = "x position", range = c(0, 760)),
    yaxis = list(title = "y position", range = c(0, 430), scaleanchor = "x"),
    legend = list(orientation = "h", x = 0.05, y = -0.12),
    margin = list(l = 50, r = 20, b = 80, t = 60)
  ) |>
  plotly::animation_opts(frame = 650, transition = 250, redraw = FALSE) |>
  plotly::animation_slider(currentvalue = list(prefix = "Step: ")) |>
  plotly::animation_button(x = 1, xanchor = "right", y = 0, yanchor = "bottom")

5.1 What to observe

  1. During Predict, the particles spread according to the motion model.
  2. During Update, particles near the sensor reading become larger because their likelihood is higher.
  3. During Resample, good hypotheses are duplicated and poor hypotheses disappear.
  4. With too few particles, the filter may collapse onto the wrong region.
  5. With too large sensor noise, measurements become weak; with too small sensor noise, the filter may overtrust a bad sensor.

6 A rigorous one-dimensional nonlinear example

A famous nonlinear filtering example uses

\[ X_t = 0.5X_{t-1}+\frac{25X_{t-1}}{1+X_{t-1}^2}+8\cos(1.2t)+\eta_t, \]

\[ Z_t = \frac{X_t^2}{20}+\epsilon_t. \]

This is difficult because the observation depends on \(X_t^2\). Therefore a positive and a negative state may produce similar observations. A Gaussian linear filter is not naturally suited for this.

simulate_nonlinear_ssm <- function(T = 80, process_sd = 1.2, obs_sd = 1.8) {
  x <- numeric(T)
  z <- numeric(T)
  x[1] <- rnorm(1, 0, 2)
  z[1] <- x[1]^2 / 20 + rnorm(1, 0, obs_sd)
  for (t in 2:T) {
    x[t] <- 0.5 * x[t - 1] + 25 * x[t - 1] / (1 + x[t - 1]^2) + 8 * cos(1.2 * t) + rnorm(1, 0, process_sd)
    z[t] <- x[t]^2 / 20 + rnorm(1, 0, obs_sd)
  }
  tibble::tibble(t = seq_len(T), x = x, z = z)
}

particle_filter_1d <- function(z, N = 2000, process_sd = 1.2, obs_sd = 1.8) {
  T <- length(z)
  particles <- rnorm(N, 0, 5)
  weights <- rep(1 / N, N)

  out <- tibble::tibble(
    t = integer(), mean = double(), q05 = double(), q50 = double(), q95 = double(), ess = double()
  )

  systematic_resample <- function(w) {
    # Robust systematic resampling.
    # This protects against tiny numerical errors in cumulative weights
    # and against accidental NA/Inf weights after likelihood underflow.
    N <- length(w)
    w[!is.finite(w)] <- 0
    sw <- sum(w)
    if (!is.finite(sw) || sw <= 0) {
      w <- rep(1 / N, N)
    } else {
      w <- w / sw
    }

    cw <- cumsum(w)
    cw <- cummax(cw)  # enforce non-decreasing cumulative probabilities
    cw[N] <- 1        # ensure the last interval closes exactly at 1

    u0 <- runif(1, 0, 1 / N)
    u <- u0 + (0:(N - 1)) / N
    # Avoid findInterval() failures on hosted knitting environments.
    idx <- integer(N)
    j <- 1L
    for (i in seq_len(N)) {
      while (u[i] > cw[j] && j < N) {
        j <- j + 1L
      }
      idx[i] <- j
    }
    idx
  }

  for (tt in seq_len(T)) {
    if (tt > 1) {
      particles <- 0.5 * particles + 25 * particles / (1 + particles^2) + 8 * cos(1.2 * tt) + rnorm(N, 0, process_sd)
    }

    logw <- dnorm(z[tt], mean = particles^2 / 20, sd = obs_sd, log = TRUE)
    logw <- logw - max(logw)
    weights <- exp(logw)
    sw <- sum(weights)
    if (!is.finite(sw) || sw <= 0) {
      weights <- rep(1 / N, N)
    } else {
      weights <- weights / sw
    }

    ess <- 1 / sum(weights^2)

    # weighted quantiles by resampling an auxiliary representation
    idx_tmp <- sample(seq_len(N), size = N, replace = TRUE, prob = weights)
    px <- particles[idx_tmp]

    out <- dplyr::bind_rows(out, tibble::tibble(
      t = tt,
      mean = sum(weights * particles),
      q05 = as.numeric(quantile(px, 0.05)),
      q50 = as.numeric(quantile(px, 0.50)),
      q95 = as.numeric(quantile(px, 0.95)),
      ess = ess
    ))

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

  out
}

ssm <- simulate_nonlinear_ssm(T = 90)
pf1 <- particle_filter_1d(ssm$z, N = 2500)
ssm_pf <- dplyr::left_join(ssm, pf1, by = "t")
ggplot(ssm_pf, aes(x = t)) +
  geom_ribbon(aes(ymin = q05, ymax = q95), fill = "skyblue", alpha = 0.28) +
  geom_line(aes(y = x, color = "Hidden true state"), linewidth = 0.9) +
  geom_line(aes(y = mean, color = "Particle-filter mean"), linewidth = 0.9) +
  labs(
    title = "Nonlinear/non-Gaussian filtering by particle filter",
    subtitle = "The shaded band is a 90% particle-filter uncertainty band",
    x = "time", y = "state"
  ) +
  scale_color_manual(values = c("Hidden true state" = "black", "Particle-filter mean" = "red")) +
  theme_minimal(base_size = 13)

ggplot(pf1, aes(x = t, y = ess)) +
  geom_line(color = "purple", linewidth = 0.9) +
  geom_hline(yintercept = 2500 / 2, linetype = 2, color = "gray40") +
  labs(
    title = "Effective sample size over time",
    subtitle = "Low ESS means particle degeneracy; resampling restores diversity",
    x = "time", y = "ESS"
  ) +
  theme_minimal(base_size = 13)

7 Counterexample: what if the filter is misspecified?

Now we deliberately use the wrong observation noise. The true observation noise is larger, but the filter acts as if the sensor is very precise. It overtrusts noisy observations.

pf_good <- particle_filter_1d(ssm$z, N = 2000, process_sd = 1.2, obs_sd = 1.8) |>
  dplyr::mutate(model = "Correct observation noise")

pf_bad <- particle_filter_1d(ssm$z, N = 2000, process_sd = 1.2, obs_sd = 0.45) |>
  dplyr::mutate(model = "Misspecified: overtrust sensor")

compare_pf <- dplyr::bind_rows(pf_good, pf_bad) |>
  dplyr::left_join(ssm |> dplyr::select(t, x), by = "t") |>
  dplyr::mutate(abs_error = abs(mean - x))

compare_pf |>
  dplyr::group_by(model) |>
  dplyr::summarise(
    RMSE = sqrt(mean((mean - x)^2)),
    mean_abs_error = mean(abs_error),
    mean_ESS = mean(ess),
    .groups = "drop"
  ) |>
  knitr::kable(digits = 3)
model RMSE mean_abs_error mean_ESS
Correct observation noise 4.422 2.229 1281.202
Misspecified: overtrust sensor 5.636 3.042 472.224 
ggplot(compare_pf, aes(x = t)) +
  geom_line(aes(y = x), color = "black", linewidth = 0.9) +
  geom_line(aes(y = mean, color = model), linewidth = 0.8) +
  labs(
    title = "Filtering failure under misspecification",
    subtitle = "A filter can be mathematically correct for the wrong model and scientifically misleading",
    x = "time", y = "state"
  ) +
  theme_minimal(base_size = 13)

7.1 Lesson

This is where statistics and machine learning meet.

  • ML often supplies flexible prediction machinery.
  • Statistics asks: what is the data-generating mechanism, what is uncertainty, and what happens under misspecification?
  • Filtering combines both: an algorithmic engine plus a probabilistic model.

8 Real data example 1: Nile river flow as a hidden hydrological regime

The built-in Nile data contain annual measurements of the flow of the Nile at Aswan from 1871 to 1970.

A simple scientific story is:

\[ Y_t = \text{observed river flow}, \qquad \mu_t = \text{hidden hydrological level}. \]

A local-level state-space model is

\[ \mu_t = \mu_{t-1}+\eta_t, \qquad Y_t = \mu_t+\epsilon_t. \]

Here \(\mu_t\) evolves slowly, while \(Y_t\) has measurement and year-to-year noise.

kalman_local_level <- function(y, q, r, m0 = y[1], C0 = 10000) {
  n <- length(y)
  m <- C <- a <- R <- innov <- Q <- numeric(n)
  loglik <- 0

  for (tt in seq_len(n)) {
    if (tt == 1) {
      a[tt] <- m0
      R[tt] <- C0 + q
    } else {
      a[tt] <- m[tt - 1]
      R[tt] <- C[tt - 1] + q
    }

    innov[tt] <- y[tt] - a[tt]
    Q[tt] <- R[tt] + r
    A <- R[tt] / Q[tt]
    m[tt] <- a[tt] + A * innov[tt]
    C[tt] <- (1 - A) * R[tt]
    loglik <- loglik + dnorm(innov[tt], 0, sqrt(Q[tt]), log = TRUE)
  }

  tibble::tibble(
    t = seq_len(n),
    filtered_mean = m,
    filtered_sd = sqrt(C),
    lwr = m - 1.96 * sqrt(C),
    upr = m + 1.96 * sqrt(C),
    loglik = loglik
  )
}

nile_df <- tibble::tibble(
  year = as.numeric(time(Nile)),
  flow = as.numeric(Nile)
)

# Simple empirical choices. For teaching we avoid optimization and choose interpretable values.
r_obs <- var(diff(nile_df$flow), na.rm = TRUE) / 2
q_state <- r_obs * 0.08

nile_fit <- kalman_local_level(nile_df$flow, q = q_state, r = r_obs) |>
  dplyr::bind_cols(nile_df)

knitr::kable(head(nile_fit, 8), digits = 2)
t filtered_mean filtered_sd lwr upr loglik year flow
1 1120.00 78.91 965.33 1274.67 -638.51 1871 1120
2 1133.69 69.56 997.35 1270.03 -638.51 1872 1160
3 1083.01 64.78 956.03 1209.99 -638.51 1873 963
4 1117.77 62.20 995.85 1239.69 -638.51 1874 1210
5 1128.81 60.77 1009.69 1247.92 -638.51 1875 1160
6 1136.74 59.97 1019.20 1254.29 -638.51 1876 1160
7 1055.60 59.52 938.95 1172.26 -638.51 1877 813
8 1098.94 59.26 982.78 1215.09 -638.51 1878 1230 
ggplot(nile_fit, aes(x = year)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "steelblue", alpha = 0.22) +
  geom_line(aes(y = flow, color = "Observed flow"), linewidth = 0.75) +
  geom_line(aes(y = filtered_mean, color = "Filtered hidden level"), linewidth = 1.0) +
  labs(
    title = "Nile river flow: filtering a hidden hydrological level",
    subtitle = "The band shows filtering uncertainty, not just a smoothed curve",
    x = "year", y = "flow"
  ) +
  scale_color_manual(values = c("Observed flow" = "gray35", "Filtered hidden level" = "red")) +
  theme_minimal(base_size = 13)

8.1 Scientific interpretation

The filtered curve is not merely a moving average. It is the estimated hidden state \(\mu_t\) under an explicit stochastic model. The uncertainty band tells us when the latent hydrological level is well or poorly determined.

9 Real data example 2: atmospheric CO2 and model misspecification

The built-in co2 series records monthly atmospheric CO2 concentration at Mauna Loa.

A naive local-level filter is misspecified because the data contain:

  1. long-term trend;
  2. seasonal cycle;
  3. short-term noise.

A simple filter without seasonality will confuse seasonal oscillation with state noise.

co2_df <- tibble::tibble(
  t = seq_along(co2),
  year = as.numeric(time(co2)),
  month = factor(cycle(co2)),
  co2 = as.numeric(co2)
)

p_raw <- ggplot(co2_df, aes(x = year, y = co2)) +
  geom_line(color = "darkgreen", linewidth = 0.7) +
  labs(title = "Monthly atmospheric CO2 at Mauna Loa", x = "year", y = "CO2") +
  theme_minimal(base_size = 13)
p_raw

9.1 Naive filter: ignores seasonality

r_co2 <- var(diff(co2_df$co2), na.rm = TRUE) / 2
q_co2 <- r_co2 * 0.05
co2_naive <- kalman_local_level(co2_df$co2, q = q_co2, r = r_co2) |>
  dplyr::bind_cols(co2_df) |>
  dplyr::mutate(model = "Naive local level")

ggplot(co2_naive, aes(x = year)) +
  geom_line(aes(y = co2), color = "gray40", linewidth = 0.6) +
  geom_line(aes(y = filtered_mean), color = "red", linewidth = 1) +
  labs(
    title = "Misspecified filter: seasonal signal is mistaken for random movement",
    x = "year", y = "CO2"
  ) +
  theme_minimal(base_size = 13)

9.2 Better simple model: remove seasonal component, then filter trend

seasonal_model <- lm(co2 ~ month, data = co2_df)
co2_df$season_adj <- residuals(seasonal_model) + mean(co2_df$co2)

r_adj <- var(diff(co2_df$season_adj), na.rm = TRUE) / 2
q_adj <- r_adj * 0.08
co2_adj_fit <- kalman_local_level(co2_df$season_adj, q = q_adj, r = r_adj) |>
  dplyr::bind_cols(co2_df)

co2_summary <- tibble::tibble(
  Model = c("Naive local-level filter", "Seasonally adjusted local-level filter"),
  Residual_SD = c(
    sd(co2_naive$co2 - co2_naive$filtered_mean),
    sd(co2_adj_fit$season_adj - co2_adj_fit$filtered_mean)
  )
)
knitr::kable(co2_summary, digits = 3)
Model Residual_SD
Naive local-level filter 1.711
Seasonally adjusted local-level filter 0.387 
ggplot(co2_adj_fit, aes(x = year)) +
  geom_line(aes(y = season_adj), color = "gray55", linewidth = 0.55) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "orange", alpha = 0.22) +
  geom_line(aes(y = filtered_mean), color = "darkred", linewidth = 1.0) +
  labs(
    title = "After seasonal adjustment: filtering the underlying CO2 trend",
    subtitle = "A better statistical model separates trend from seasonal structure",
    x = "year", y = "seasonally adjusted CO2"
  ) +
  theme_minimal(base_size = 13)

9.3 Lesson

A flexible ML curve might fit this dataset visually, but the scientific question is sharper:

  • What is trend?
  • What is seasonality?
  • What is noise?
  • What is uncertainty?

Statistical modelling forces these distinctions.

10 Real data example 3: sunspots and time-order validation

Solar activity is a dynamical process. If we randomly split a time series, the future leaks into the training distribution. For forecasting, chronological validation is more honest.

sun_df <- tibble::tibble(
  year = as.numeric(time(sunspot.year)),
  sunspot = as.numeric(sunspot.year)
) |>
  dplyr::mutate(
    lag1 = dplyr::lag(sunspot, 1),
    lag2 = dplyr::lag(sunspot, 2),
    lag3 = dplyr::lag(sunspot, 3)
  ) |>
  tidyr::drop_na()

train <- sun_df |> dplyr::filter(year <= 1950)
test  <- sun_df |> dplyr::filter(year > 1950)

lm_ar <- lm(sunspot ~ lag1 + lag2 + lag3, data = train)
tree_ar <- rpart::rpart(sunspot ~ lag1 + lag2 + lag3, data = train, control = rpart::rpart.control(cp = 0.01))

test_pred <- test |>
  dplyr::mutate(
    pred_linear_AR = as.numeric(predict(lm_ar, newdata = test)),
    pred_tree_ML = as.numeric(predict(tree_ar, newdata = test))
  )

rmse <- function(a, b) sqrt(mean((a - b)^2))

perf <- tibble::tibble(
  Model = c("Statistical AR regression", "ML regression tree"),
  Test_RMSE = c(
    rmse(test_pred$sunspot, test_pred$pred_linear_AR),
    rmse(test_pred$sunspot, test_pred$pred_tree_ML)
  )
)
knitr::kable(perf, digits = 3)
Model Test_RMSE
Statistical AR regression 24.462
ML regression tree 34.648 
ggplot(test_pred, aes(x = year)) +
  geom_line(aes(y = sunspot, color = "Observed"), linewidth = 0.9) +
  geom_line(aes(y = pred_linear_AR, color = "Statistical AR"), linewidth = 0.85) +
  geom_line(aes(y = pred_tree_ML, color = "ML tree"), linewidth = 0.85) +
  labs(
    title = "Sunspot forecasting: statistics and ML under chronological validation",
    subtitle = "The task is not random interpolation; it is prediction of future dynamics",
    x = "year", y = "sunspot number"
  ) +
  theme_minimal(base_size = 13)

10.1 Interpretation

  • The statistical AR model is transparent and gives a simple dynamical baseline.
  • The ML tree can capture nonlinear threshold patterns, but may be unstable outside the training regime.
  • A serious scientific workflow should use both: interpretable stochastic baselines and flexible ML competitors.

11 Where ML helps, where statistics helps

comp <- tibble::tribble(
  ~Question, ~Statistical_inference, ~Machine_learning, ~Best_scientific_practice,
  "What is the hidden state?", "Posterior/filtering distribution", "Latent representation or predictor", "Use a probabilistic state-space model with flexible components",
  "How uncertain is the answer?", "Confidence/credible/prediction intervals", "Often not automatic", "Calibrate ML using statistical uncertainty tools",
  "Can the model extrapolate?", "Model assumptions visible", "Can fail silently under shift", "Use physical constraints and honest validation",
  "Is the mechanism interpretable?", "Parameters tied to scientific quantities", "High predictive accuracy possible", "Combine interpretable structure with flexible residual learning",
  "What if the model is wrong?", "Model checking and sensitivity", "Stress testing and external validation", "Do both"
)
knitr::kable(comp)
Question Statistical_inference Machine_learning Best_scientific_practice
What is the hidden state? Posterior/filtering distribution Latent representation or predictor Use a probabilistic state-space model with flexible components
How uncertain is the answer? Confidence/credible/prediction intervals Often not automatic Calibrate ML using statistical uncertainty tools
Can the model extrapolate? Model assumptions visible Can fail silently under shift Use physical constraints and honest validation
Is the mechanism interpretable? Parameters tied to scientific quantities High predictive accuracy possible Combine interpretable structure with flexible residual learning
What if the model is wrong? Model checking and sensitivity Stress testing and external validation Do both

12 A hybrid philosophy: ML + probabilistic inference

A modern scientific data-science pipeline may look like this:

  1. Use scientific knowledge to specify state variables and measurement structure.
  2. Use statistics to define uncertainty, likelihood, validation, and inferential targets.
  3. Use ML to learn flexible nonlinear components.
  4. Use filtering or Bayesian/sequential inference to propagate uncertainty over time.
  5. Use model checking to avoid confident nonsense.

In formula form, one may replace a rigid transition model by a learned component:

\[ X_t = f_{\mathrm{ML}}(X_{t-1},u_t)+\eta_t, \]

but still keep an observation model

\[ Z_t\sim p(z_t\mid X_t), \]

and still infer

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

That is the useful synthesis:

ML learns complicated structure; statistics keeps uncertainty honest.

13 Final take-home messages

  1. A particle filter represents uncertainty by a cloud of weighted simulated states.
  2. It is useful when the system is nonlinear, non-Gaussian, or multi-modal.
  3. Filtering is a probability statement, not just a smoothing trick.
  4. Misspecification can make a filter confidently wrong.
  5. Real scientific data often require decomposition into trend, seasonality, noise, and hidden state.
  6. Statistics and ML are complementary: ML gives flexible prediction; statistics gives uncertainty, assumptions, estimands, and diagnostics.

14 Suggested further reading