particle_filter {bayesSSM} | R Documentation |
Particle Filter
Description
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
|
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 |
algorithm |
A character string specifying the particle filtering
algorithm to use. Must be one of |
resample_fn |
A character string specifying the resampling method.
Must be one of |
threshold |
A numeric value specifying the ESS threshold for triggering
resampling in the |
return_particles |
A logical value indicating whether to return the full
particle history. Defaults to |
... |
Additional arguments passed to |
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.
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 ifreturn_particles
isTRUE
.- weights_history
(Optional) A matrix of particle weights over time, with dimension
(num_obs + 1) x num_particles
. Returned ifreturn_particles
isTRUE
.
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)