Particle Filter — particle_filter • bayesSSM

This function implements a bootstrap particle filter for estimating the hidden states in a state space model using sequential Monte Carlo methods. Three filtering variants are supported:

SIS: Sequential Importance Sampling (without resampling).
SISR: Sequential Importance Sampling with resampling at every time step.
SISAR: SIS with adaptive resampling based on the Effective Sample Size (ESS). Resampling is triggered when the ESS falls below a given threshold (default particles / 2).

It is recommended to use either SISR or SISAR to avoid weight degeneracy.

Usage

particle_filter(
  y,
  num_particles,
  init_fn,
  transition_fn,
  log_likelihood_fn,
  obs_times = NULL,
  algorithm = c("SISAR", "SISR", "SIS"),
  resample_fn = c("stratified", "systematic", "multinomial"),
  threshold = NULL,
  return_particles = TRUE,
  ...
)

Arguments

y: A numeric vector or matrix of observations. Each row represents an observation at a time step. If observations not equally spaced, use the obs_times argument to specify the time points at which observations are available.
num_particles: A positive integer specifying the number of particles.
init_fn: A function that initializes the particle states. It should take `num_particles` as an argument for initializing the particles and return a vector or matrix of initial particle states. It can include any model-specific parameters as named arguments.
transition_fn: A function describing the state transition model. It should take `particles` as an argument and return the propagated particles. The function can optionally depend on time by including a time step argument `t`. It can include any model-specific parameters as named arguments.
log_likelihood_fn: A function that computes the log-likelihoods for the particles. It should take a `y` argument for the observations, the current particles, and return a numeric vector of log-likelihood values. The function can optionally depend on time by including a time step argument `t`. It can include any model-specific parameters as named arguments.
obs_times: A numeric vector indicating the time points at which observations in y are available. Must be of the same length as the number of rows in y. If not specified, it is assumed that observations are available at consecutive time steps, i.e., obs_times = 1:nrow(y).
algorithm: A character string specifying the particle filtering algorithm to use. Must be one of "SISAR", "SISR", or "SIS". Defaults to "SISAR".
resample_fn: A character string specifying the resampling method. Must be one of "stratified", "systematic", or "multinomial". Defaults to "stratified".
threshold: A numeric value specifying the ESS threshold for triggering resampling in the "SISAR" algorithm. If not provided, it defaults to num_particles / 2.
return_particles: A logical value indicating whether to return the full particle history. Defaults to TRUE.
...: Additional arguments passed to init_fn, transition_fn, and log_likelihood_fn. I.e., parameter values if the functions requires them.

Value

A list containing:

state_est: A numeric vector of estimated states over time, computed as the weighted average of particles.
ess: A numeric vector of the Effective Sample Size (ESS) at each time step.
loglike: The accumulated log-likelihood of the observations given the model.
loglike_history: A numeric vector of the log-likelihood at each time step.
algorithm: A character string indicating the filtering algorithm used.
particles_history: (Optional) A matrix of particle states over time, with dimension (num_obs + 1) x num_particles. Returned if return_particles is TRUE.
weights_history: (Optional) A matrix of particle weights over time, with dimension (num_obs + 1) x num_particles. Returned if return_particles is TRUE.

Details

The particle filter is a sequential Monte Carlo method that approximates the posterior distribution of the state in a state space model. The three supported algorithms differ in their approach to resampling:

SIS: Particles are propagated and weighted without any resampling, which may lead to weight degeneracy over time.
SISR: Resampling is performed at every time step to combat weight degeneracy.
SISAR: Resampling is performed adaptively; particles are resampled only when the Effective Sample Size (ESS) falls below a specified threshold (defaulting to particles / 2).

The Effective Sample Size (ESS) in context of particle filters is defined as $$ESS = \left(\sum_{i=1}^{\text{n}} w_i^2\right)^{-1},$$ where $n$ is the number of particles and $w_i$ are the normalized weights of the particles.

The default resampling method is stratified resampling, as Douc et al., 2005 showed that it always gives a lower variance compared to multinomial resampling.

References

Douc, R., Cappé, O., & Moulines, E. (2005). Comparison of Resampling Schemes for Particle Filtering. Accessible at: https://arxiv.org/abs/cs/0507025

Examples

init_fn <- function(num_particles) rnorm(num_particles, 0, 1)
transition_fn <- function(particles) particles + rnorm(length(particles))
log_likelihood_fn <- function(y, particles) {
  dnorm(y, mean = particles, sd = 1, log = TRUE)
}

y <- cumsum(rnorm(50)) # dummy data
num_particles <- 100

# Run the particle filter using default settings.
result <- particle_filter(
  y = y,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn
)
plot(result$state_est, type = "l", col = "blue", main = "State Estimates",
  ylim = range(c(result$state_est, y)))
points(y, col = "red", pch = 20)


# With parameters
init_fn <- function(num_particles) rnorm(num_particles, 0, 1)
transition_fn <- function(particles, mu) {
  particles + rnorm(length(particles), mean = mu)
}
log_likelihood_fn <- function(y, particles, sigma) {
  dnorm(y, mean = particles, sd = sigma, log = TRUE)
}

y <- cumsum(rnorm(50)) # dummy data
num_particles <- 100

# Run the particle filter using default settings.
result <- particle_filter(
  y = y,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  mu = 1,
  sigma = 1
)
plot(result$state_est, type = "l", col = "blue", main = "State Estimates",
  ylim = range(c(result$state_est, y)))
points(y, col = "red", pch = 20)


# With observations gaps
init_fn <- function(num_particles) rnorm(num_particles, 0, 1)
transition_fn <- function(particles, mu) {
  particles + rnorm(length(particles), mean = mu)
}
log_likelihood_fn <- function(y, particles, sigma) {
  dnorm(y, mean = particles, sd = sigma, log = TRUE)
}

# Generate data using DGP
simulate_ssm <- function(num_steps, mu, sigma) {
  x <- numeric(num_steps)
  y <- numeric(num_steps)
  x[1] <- rnorm(1, mean = 0, sd = sigma)
  y[1] <- rnorm(1, mean = x[1], sd = sigma)
  for (t in 2:num_steps) {
    x[t] <- mu * x[t - 1] + sin(x[t - 1]) + rnorm(1, mean = 0, sd = sigma)
    y[t] <- x[t] + rnorm(1, mean = 0, sd = sigma)
  }
  y
}

data <- simulate_ssm(10, mu = 1, sigma = 1)
# Suppose we have data for t=1,2,3,5,6,7,8,9,10 (i.e., missing at t=4)

obs_times <- c(1, 2, 3, 5, 6, 7, 8, 9, 10)
data_obs <- data[obs_times]

num_particles <- 100
# Run the particle filter
# Specify observation times in the particle filter using obs_times
result <- particle_filter(
  y = data_obs,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  obs_times = obs_times,
  mu = 1,
  sigma = 1,
)
plot(result$state_est, type = "l", col = "blue", main = "State Estimates",
  ylim = range(c(result$state_est, data)))
points(data_obs, col = "red", pch = 20)