bayesSSM: Bayesian Inference for State-Space Models

Description

The bayesSSM package provides implementations of particle filtering, Particle MCMC, and related methods for Bayesian inference in state-space models. It includes tools for simulation, posterior inference, and diagnostics.

Model Specification

Particle filter implementations in this package assume a discrete-time state-space model defined by:

• A sequence of latent states x_0, x_1, \ldots, x_T evolving according to a Markov process.

• Observations y_1, \ldots, y_T that are conditionally independent given the corresponding latent states.

The model is specified as:

x_0 \sim \mu_\theta

x_t \sim f_\theta(x_t \mid x_{t-1}), \quad t = 1, \ldots, T

y_t \sim g_\theta(y_t \mid x_t), \quad t = 1, \ldots, T

where \theta denotes model parameters passed via ....

The user provides the following functions:

• init_fn: draws from the initial distribution \mu_\theta.

• transition_fn: generates or evaluates the transition density f_\theta.

• weight_fn: evaluates the observation likelihood g_\theta.

Author(s)

Maintainer: Bjarke Hautop bjarke.hautop@gmail.com [copyright holder]

See Also

Useful links:

• https://github.com/BjarkeHautop/bayesSSM

• https://bjarkehautop.github.io/bayesSSM/

• Report bugs at https://github.com/BjarkeHautop/bayesSSM/issues

Internal function to back-transform parameters

Description

Internal function to back-transform parameters

Usage

.back_transform_params(theta_trans, transform)

Arguments

Value

original parameter vector

Helper function to validate input of user-defined functions and priors

Description

Helper function to validate input of user-defined functions and priors

Usage

.check_params_match(
  init_fn,
  transition_fn,
  log_likelihood_fn,
  pilot_init_params,
  log_priors
)

Arguments

Internal function to compute the Jacobian of the transformation

Description

Internal function to compute the Jacobian of the transformation

Usage

.compute_log_jacobian(theta, transform)

Arguments

Value

log-Jacobian of the transformation

Ensure that a function has a '...' argument

Description

Ensure that a function has a '...' argument

Usage

.ensure_dots(fun)

Arguments

Value

The modified function with '...' added to its formals if it was not already present.

Core Particle Filter Function

Description

This function implements the underlying logic used for particle filters in a state space model using sequential Monte Carlo methods.

Usage

.particle_filter_core(
  y,
  num_particles,
  init_fn,
  transition_fn,
  weight_fn,
  aux_weight_fn = NULL,
  move_fn = NULL,
  obs_times = NULL,
  algorithm = c("BPF", "APF", "RMPF"),
  resample_algorithm = c("SIS", "SISR", "SISAR"),
  resample_fn = c("stratified", "systematic", "multinomial"),
  threshold = NULL,
  return_particles = TRUE,
  ...
)

Arguments

Value

A list with components:

state_est

Estimated states over time (weighted mean of particles).

ess

Effective sample size at each time step.

loglike

Total log-likelihood.

loglike_history

Log-likelihood at each time step.

algorithm

The filtering algorithm used.

particles_history

Matrix of particle states over time (if return_particles = TRUE).

weights_history

Matrix of particle weights over time (if return_particles = TRUE).

Model Specification

Particle filter implementations in this package assume a discrete-time state-space model defined by:

• A sequence of latent states x_0, x_1, \ldots, x_T evolving according to a Markov process.

• Observations y_1, \ldots, y_T that are conditionally independent given the corresponding latent states.

The model is specified as:

x_0 \sim \mu_\theta

x_t \sim f_\theta(x_t \mid x_{t-1}), \quad t = 1, \ldots, T

y_t \sim g_\theta(y_t \mid x_t), \quad t = 1, \ldots, T

where \theta denotes model parameters passed via ....

The user provides the following functions:

• init_fn: draws from the initial distribution \mu_\theta.

• transition_fn: generates or evaluates the transition density f_\theta.

• weight_fn: evaluates the observation likelihood g_\theta.

Pilot Run for Particle Filter Tuning

Description

This internal function repeatedly evaluates the particle filter in order to estimate the variance of the log-likelihoods and to compute a recommended target number of particles for the Particle Marginal Metropolis Hastings (PMMH) algorithm.

Usage

.pilot_run(
  pf_wrapper,
  y,
  pilot_n,
  pilot_reps,
  init_fn,
  transition_fn,
  log_likelihood_fn,
  obs_times = NULL,
  resample_fn = NULL,
  ...
)

Arguments

Details

The function performs pilot_reps evaluations of the particle filter using the provided parameter vector theta. It then estimates the variance of the log-likelihoods and scales the initial particle number by this variance. The final number of particles is taken as the ceiling of the scaled value with a minimum of 50 and a maximum of 1000.

Value

A list containing:

variance_estimate

The estimated variance of the log-likelihoods from the pilot run.

target_N

The number of particles used in PMMH algorithm.

pilot_loglikes

A numeric vector of log-likelihood values computed during the run.

Run Pilot Chain for Posterior Estimation

Description

Run Pilot Chain for Posterior Estimation

Usage

.run_pilot_chain(
  pf_wrapper,
  y,
  pilot_m,
  pilot_n,
  pilot_reps,
  init_fn,
  transition_fn,
  log_likelihood_fn,
  log_priors,
  proposal_sd,
  obs_times = NULL,
  param_transform = NULL,
  pilot_init_params = NULL,
  verbose = FALSE,
  ...
)

Arguments

Details

This function runs a pilot chain to estimate the posterior mean and covariance of the model parameters using a particle filter. The chain is run for 'pilot_m' iterations, with each iteration proposing new parameters and evaluating their likelihood and prior. The chain is then used to estimate the posterior mean and covariance, which are used to tune the number of particles for the Particle Marginal Metropolis Hastings (PMMH) algorithm.

Value

A list containing:

Internal function to transform parameters

Description

Internal function to transform parameters

Usage

.transform_params(theta, transform)

Arguments

Value

transformed parameter vector

Auxiliary Particle Filter (APF)

Description

The Auxiliary Particle Filter differs from the bootstrap filter by incorporating a look-ahead step: particles are reweighted using an approximation of the likelihood of the next observation prior to resampling. This adjustment can help reduce particle degeneracy and, improve filtering efficiency compared to the bootstrap approach.

Usage

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

Arguments

Value

A list with components:

state_est

Estimated states over time (weighted mean of particles).

ess

Effective sample size at each time step.

loglike

Total log-likelihood.

loglike_history

Log-likelihood at each time step.

algorithm

The filtering algorithm used.

particles_history

Matrix of particle states over time (if return_particles = TRUE).

weights_history

Matrix of particle weights over time (if return_particles = TRUE).

The Auxiliary Particle Filter (APF)

The Auxiliary Particle Filter (APF) was introduced by Pitt and Shephard (1999) to improve upon the standard bootstrap filter by incorporating a look ahead step. Before resampling at time t, particles are weighted by an auxiliary weight proportional to an estimate of the likelihood of the next observation, guiding resampling to favor particles likely to contribute to future predictions.

Specifically, if w_{t-1}^i are the normalized weights and x_{t-1}^i are the particles at time t-1, then auxiliary weights are computed as

\tilde{w}_t^i \propto w_{t-1}^i \, p(y_t | \mu_t^i),

where \mu_t^i is a predictive summary (e.g., the expected next state) of the particle x_{t-1}^i. Resampling is performed using \tilde{w}_t^i instead of w_{t-1}^i. This can reduce the variance of the importance weights at time t and help mitigate particle degeneracy, especially if the auxiliary weights are chosen well.

Default resampling method is stratified resampling, which has lower variance than multinomial resampling (Douc et al., 2005).

Model Specification

Particle filter implementations in this package assume a discrete-time state-space model defined by:

• A sequence of latent states x_0, x_1, \ldots, x_T evolving according to a Markov process.

• Observations y_1, \ldots, y_T that are conditionally independent given the corresponding latent states.

The model is specified as:

x_0 \sim \mu_\theta

x_t \sim f_\theta(x_t \mid x_{t-1}), \quad t = 1, \ldots, T

y_t \sim g_\theta(y_t \mid x_t), \quad t = 1, \ldots, T

where \theta denotes model parameters passed via ....

The user provides the following functions:

• init_fn: draws from the initial distribution \mu_\theta.

• transition_fn: generates or evaluates the transition density f_\theta.

• weight_fn: evaluates the observation likelihood g_\theta.

References

Pitt, M. K., & Shephard, N. (1999). Filtering via simulation: Auxiliary particle filters. Journal of the American Statistical Association, 94(446), 590–599. doi:10.1080/01621459.1999.10474153

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)
}
aux_log_likelihood_fn <- function(y, particles) {
  # Predict next state (mean stays same) and compute log p(y | x)
  mean_forecast <- particles # since E[x'] = x in this model
  dnorm(y, mean = mean_forecast, sd = 1, log = TRUE)
}

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

result <- auxiliary_filter(
  y = y,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  aux_log_likelihood_fn = aux_log_likelihood_fn
)
plot(result$state_est,
  type = "l", col = "blue", main = "APF: 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)
}
aux_log_likelihood_fn <- function(y, particles, mu, sigma) {
  # Forecast mean of x' given x, then evaluate p(y | forecast)
  forecast <- particles + mu
  dnorm(y, mean = forecast, sd = sigma, log = TRUE)
}

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

result <- auxiliary_filter(
  y = y,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  aux_log_likelihood_fn = aux_log_likelihood_fn,
  mu = 1,
  sigma = 1
)
plot(result$state_est,
  type = "l", col = "blue", main = "APF with Parameters",
  ylim = range(c(result$state_est, y))
)
points(y, col = "red", pch = 20)

# ---- With observation gaps ----
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)
obs_times <- c(1, 2, 3, 5, 6, 7, 8, 9, 10) # Missing at t = 4
data_obs <- data[obs_times]

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)
}
aux_log_likelihood_fn <- function(y, particles, mu, sigma) {
  forecast <- particles + mu
  dnorm(y, mean = forecast, sd = sigma, log = TRUE)
}

num_particles <- 100
result <- auxiliary_filter(
  y = data_obs,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  aux_log_likelihood_fn = aux_log_likelihood_fn,
  obs_times = obs_times,
  mu = 1,
  sigma = 1
)
plot(result$state_est,
  type = "l", col = "blue", main = "APF with Observation Gaps",
  ylim = range(c(result$state_est, data))
)
points(data_obs, col = "red", pch = 20)

Bootstrap Particle Filter (BPF)

Description

Implements a bootstrap particle filter for sequential Bayesian inference in state space models using sequential Monte Carlo methods.

Usage

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

Arguments

Value

A list with components:

state_est

Estimated states over time (weighted mean of particles).

ess

Effective sample size at each time step.

loglike

Total log-likelihood.

loglike_history

Log-likelihood at each time step.

algorithm

The filtering algorithm used.

particles_history

Matrix of particle states over time (if return_particles = TRUE).

weights_history

Matrix of particle weights over time (if return_particles = TRUE).

The Effective Sample Size (ESS) is defined as

ESS = \left(\sum_{i=1}^{n} w_i^2\right)^{-1},

where w_i are the normalized weights of the particles.

Default resampling method is stratified resampling, which has lower variance than multinomial resampling (Douc et al., 2005).

Model Specification

Particle filter implementations in this package assume a discrete-time state-space model defined by:

• A sequence of latent states x_0, x_1, \ldots, x_T evolving according to a Markov process.

• Observations y_1, \ldots, y_T that are conditionally independent given the corresponding latent states.

The model is specified as:

x_0 \sim \mu_\theta

x_t \sim f_\theta(x_t \mid x_{t-1}), \quad t = 1, \ldots, T

y_t \sim g_\theta(y_t \mid x_t), \quad t = 1, \ldots, T

where \theta denotes model parameters passed via ....

The user provides the following functions:

• init_fn: draws from the initial distribution \mu_\theta.

• transition_fn: generates or evaluates the transition density f_\theta.

• weight_fn: evaluates the observation likelihood g_\theta.

References

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. doi:10.1049/ip-f-2.1993.0015

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 <- bootstrap_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 bootstrap particle filter using default settings.
result <- bootstrap_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
# Specify observation times in the bootstrap particle filter using obs_times
result <- bootstrap_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)

Create Tuning Control Parameters

Description

This function creates a list of tuning parameters used by the pmmh function. The tuning choices are inspired by Pitt et al. [2012] and Dahlin and Schön [2019].

Usage

default_tune_control(
  pilot_proposal_sd = 0.5,
  pilot_n = 100,
  pilot_m = 2000,
  pilot_target_var = 1,
  pilot_burn_in = 500,
  pilot_reps = 100,
  pilot_resample_algorithm = c("SISAR", "SISR", "SIS"),
  pilot_resample_fn = c("stratified", "systematic", "multinomial")
)

Arguments

Value

A list of tuning control parameters.

References

M. K. Pitt, R. d. S. Silva, P. Giordani, and R. Kohn. On some properties of Markov chain Monte Carlo simulation methods based on the particle filter. Journal of Econometrics, 171(2):134–151, 2012. doi: https://doi.org/10.1016/j.jeconom.2012.06.004

J. Dahlin and T. B. Schön. Getting started with particle Metropolis-Hastings for inference in nonlinear dynamical models. Journal of Statistical Software, 88(2):1–41, 2019. doi: 10.18637/jss.v088.c02

Estimate effective sample size (ESS) of MCMC chains.

Description

Estimate effective sample size (ESS) of MCMC chains.

Usage

ess(chains)

Arguments

Details

Uses the formula for ESS proposed by Vehtari et al. (2021).

Value

The estimated effective sample size (ess) if given a matrix, or a named vector of ESS values if given a data frame.

References

Vehtari et al. (2021). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Available at: https://doi.org/10.1214/20-BA1221

Examples

# With a matrix:
chains <- matrix(rnorm(3000), nrow = 1000, ncol = 3)
ess(chains)

# With a data frame:
chains_df <- data.frame(
  chain = rep(1:3, each = 1000),
  param1 = rnorm(3000),
  param2 = rnorm(3000)
)
ess(chains_df)

Particle filter functions

Description

The package provides several particle filter implementations for state-space models for estimating the intractable marginal likelihood p(y_{1:T}\mid \theta):

• auxiliary_filter

• bootstrap_filter

• resample_move_filter

The simplest one is the bootstrap_filter, and is thus recommended as a starting point.

Common Parameters for Particle Filters

Description

These parameters are shared by particle filter implementations such as the bootstrap filter, auxiliary particle filter, and resample-move particle filter.

Arguments

Model Specification for Particle Filters

Description

Model Specification for Particle Filters

Model Specification

Particle filter implementations in this package assume a discrete-time state-space model defined by:

• A sequence of latent states x_0, x_1, \ldots, x_T evolving according to a Markov process.

• Observations y_1, \ldots, y_T that are conditionally independent given the corresponding latent states.

The model is specified as:

x_0 \sim \mu_\theta

x_t \sim f_\theta(x_t \mid x_{t-1}), \quad t = 1, \ldots, T

y_t \sim g_\theta(y_t \mid x_t), \quad t = 1, \ldots, T

where \theta denotes model parameters passed via ....

The user provides the following functions:

• init_fn: draws from the initial distribution \mu_\theta.

• transition_fn: generates or evaluates the transition density f_\theta.

• weight_fn: evaluates the observation likelihood g_\theta.

Shared Return Values for Particle Filters

Description

This block documents the common return value for particle filtering functions.

Value

A list with components:

state_est

Estimated states over time (weighted mean of particles).

ess

Effective sample size at each time step.

loglike

Total log-likelihood.

loglike_history

Log-likelihood at each time step.

algorithm

The filtering algorithm used.

particles_history

Matrix of particle states over time (if return_particles = TRUE).

weights_history

Matrix of particle weights over time (if return_particles = TRUE).

Particle Marginal Metropolis-Hastings (PMMH) for State-Space Models

Description

This function implements a Particle Marginal Metropolis-Hastings (PMMH) resample_algorithm to perform Bayesian inference in state-space models. It first runs a pilot chain to tune the proposal distribution and the number of particles for the particle filter, and then runs the main PMMH chain.

Usage

pmmh(
  pf_wrapper,
  y,
  m,
  init_fn,
  transition_fn,
  log_likelihood_fn,
  log_priors,
  pilot_init_params,
  burn_in,
  num_chains = 4,
  obs_times = NULL,
  resample_algorithm = c("SISAR", "SISR", "SIS"),
  resample_fn = c("stratified", "systematic", "multinomial"),
  param_transform = NULL,
  tune_control = default_tune_control(),
  verbose = FALSE,
  return_latent_state_est = FALSE,
  seed = NULL,
  num_cores = 1,
  ...
)

Arguments

Details

The PMMH resample_algorithm is essentially a Metropolis Hastings algorithm, where instead of using the intractable marginal likelihood p(y_{1:T}\mid \theta) it instead uses the estimated likelihood using a particle filter (see particle_filter for available particle filters). Values are proposed using a multivariate normal distribution in the transformed space (specified using 'param_transform'). The proposal covariance and the number of particles is chosen based on a pilot run. The number of particles is chosen such that the variance of the log-likelihood estimate at the estimated posterior mean is approximately 1 (with a minimum of 50 particles and a maximum of 1000).

Value

A list containing:

theta_chain

A dataframe of post burn-in parameter samples.

latent_state_chain

If return_latent_state_est is TRUE, a list of matrices containing the latent state estimates for each time step.

diagnostics

Diagnostics containing ESS and Rhat for each parameter (see ess and rhat for documentation).

References

Andrieu et al. (2010). Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342. doi: 10.1111/j.1467-9868.2009.00736.x

Examples

init_fn <- function(num_particles) {
  rnorm(num_particles, mean = 0, sd = 1)
}
transition_fn <- function(particles, phi, sigma_x) {
  phi * particles + sin(particles) +
    rnorm(length(particles), mean = 0, sd = sigma_x)
}
log_likelihood_fn <- function(y, particles, sigma_y) {
  dnorm(y, mean = cos(particles), sd = sigma_y, log = TRUE)
}
log_prior_phi <- function(phi) {
  dnorm(phi, mean = 0, sd = 1, log = TRUE)
}
log_prior_sigma_x <- function(sigma) {
  dexp(sigma, rate = 1, log = TRUE)
}
log_prior_sigma_y <- function(sigma) {
  dexp(sigma, rate = 1, log = TRUE)
}
log_priors <- list(
  phi = log_prior_phi,
  sigma_x = log_prior_sigma_x,
  sigma_y = log_prior_sigma_y
)
# Generate data
t_val <- 10
x <- numeric(t_val)
y <- numeric(t_val)
phi <- 0.8
sigma_x <- 1
sigma_y <- 0.5

init_state <- rnorm(1, mean = 0, sd = 1)
x[1] <- phi * init_state + sin(init_state) + rnorm(1, mean = 0, sd = sigma_x)
y[1] <- x[1] + rnorm(1, mean = 0, sd = sigma_y)
for (t in 2:t_val) {
  x[t] <- phi * x[t - 1] + sin(x[t - 1]) + rnorm(1, mean = 0, sd = sigma_x)
  y[t] <- cos(x[t]) + rnorm(1, mean = 0, sd = sigma_y)
}
x <- c(init_state, x)

# Should use higher MCMC iterations in practice (m)
pmmh_result <- pmmh(
  pf_wrapper = bootstrap_filter,
  y = y,
  m = 1000,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  log_priors = log_priors,
  pilot_init_params = list(
    c(phi = 0.8, sigma_x = 1, sigma_y = 0.5),
    c(phi = 1, sigma_x = 0.5, sigma_y = 1)
  ),
  burn_in = 100,
  num_chains = 2,
  param_transform = list(
    phi = "identity",
    sigma_x = "log",
    sigma_y = "log"
  ),
  tune_control = default_tune_control(pilot_m = 500, pilot_burn_in = 100)
)
# Convergence warning is expected with such low MCMC iterations.

# 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)
y <- y[obs_times]

# Specify observation times in the pmmh using obs_times
pmmh_result <- pmmh(
  pf_wrapper = bootstrap_filter,
  y = y,
  m = 1000,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  log_priors = log_priors,
  pilot_init_params = list(
    c(phi = 0.8, sigma_x = 1, sigma_y = 0.5),
    c(phi = 1, sigma_x = 0.5, sigma_y = 1)
  ),
  burn_in = 100,
  num_chains = 2,
  obs_times = obs_times,
  param_transform = list(
    phi = "identity",
    sigma_x = "log",
    sigma_y = "log"
  ),
  tune_control = default_tune_control(pilot_m = 500, pilot_burn_in = 100)
)

Print method for PMMH output

Description

Displays a concise summary of parameter estimates from a PMMH output object, including means, standard deviations, medians, 95% credible intervals, effective sample sizes (ESS), and Rhat. This provides a quick overview of the posterior distribution and convergence diagnostics.

Usage

## S3 method for class 'pmmh_output'
print(x, ...)

Arguments

Value

The object 'x' invisibly.

Examples

# Create dummy chains for two parameters across two chains
chain1 <- data.frame(param1 = rnorm(100), param2 = rnorm(100), chain = 1)
chain2 <- data.frame(param1 = rnorm(100), param2 = rnorm(100), chain = 2)
dummy_output <- list(
  theta_chain = rbind(chain1, chain2),
  diagnostics = list(
    ess = c(param1 = 200, param2 = 190),
    rhat = c(param1 = 1.01, param2 = 1.00)
  )
)
class(dummy_output) <- "pmmh_output"
print(dummy_output)

Resample-Move Particle Filter (RMPF)

Description

The Resample-Move Particle Filter differs from standard resampling methods by including a Metropolis–Hastings move step after resampling. This additional step can increase particle diversity and, in some contexts, help mitigate sample impoverishment.

Usage

resample_move_filter(
  y,
  num_particles,
  init_fn,
  transition_fn,
  log_likelihood_fn,
  move_fn,
  obs_times = NULL,
  resample_fn = c("stratified", "systematic", "multinomial"),
  return_particles = TRUE,
  ...
)

Arguments

Value

A list with components:

state_est

Estimated states over time (weighted mean of particles).

ess

Effective sample size at each time step.

loglike

Total log-likelihood.

loglike_history

Log-likelihood at each time step.

algorithm

The filtering algorithm used.

particles_history

Matrix of particle states over time (if return_particles = TRUE).

weights_history

Matrix of particle weights over time (if return_particles = TRUE).

The Resample-Move Particle Filter (RMPF)

The Resample-Move Particle Filter enhances the standard particle filtering framework by introducing a move step after resampling. After resampling at time t, particles \{x_t^{(i)}\}_{i=1}^N are propagated via a Markov kernel K_t(x' \mid x) that leaves the target posterior p(x_t \mid y_{1:t}) invariant:

x_t^{(i)} \sim K_t(\cdot \mid x_t^{(i)}).

This move step often uses a Metropolis-Hastings update that preserves the posterior distribution as the invariant distribution of K_t.

The goal of the move step is to mitigate particle impoverishment — the collapse of diversity caused by resampling selecting only a few unique particles — by rejuvenating particles and exploring the state space more thoroughly. This leads to improved approximation of the filtering distribution and reduces Monte Carlo error.

The move_fn argument represents this transition kernel and should take the current particle set as input and return the updated particles. Additional model-specific parameters may be passed via ....

Default resampling method is stratified resampling, which has lower variance than multinomial resampling (Douc et al., 2005).

In this implementation, resampling is performed at every time step using the specified method (default: stratified), followed immediately by the move step. This follows the standard Resample-Move framework as described by Gilks and Berzuini (2001). Unlike other particle filtering variants that may use an ESS threshold to decide whether to resample, RMPF requires resampling at every step to ensure the effectiveness of the subsequent rejuvenation step.

Model Specification

Particle filter implementations in this package assume a discrete-time state-space model defined by:

• A sequence of latent states x_0, x_1, \ldots, x_T evolving according to a Markov process.

• Observations y_1, \ldots, y_T that are conditionally independent given the corresponding latent states.

The model is specified as:

x_0 \sim \mu_\theta

x_t \sim f_\theta(x_t \mid x_{t-1}), \quad t = 1, \ldots, T

y_t \sim g_\theta(y_t \mid x_t), \quad t = 1, \ldots, T

where \theta denotes model parameters passed via ....

The user provides the following functions:

• init_fn: draws from the initial distribution \mu_\theta.

• transition_fn: generates or evaluates the transition density f_\theta.

• weight_fn: evaluates the observation likelihood g_\theta.

References

Gilks, W. R., & Berzuini, C. (2001). Following a moving target—Monte Carlo inference for dynamic Bayesian models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(1), 127–146. doi:10.2307/2670179

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)
}

# Define a simple random-walk Metropolis move function
move_fn <- function(particle, y) {
  proposal <- particle + rnorm(1, 0, 0.1)
  log_p_current <- log_likelihood_fn(y = y, particles = particle)
  log_p_proposal <- log_likelihood_fn(y = y, particles = proposal)
  if (log(runif(1)) < (log_p_proposal - log_p_current)) {
    return(proposal)
  } else {
    return(particle)
  }
}

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

result <- resample_move_filter(
  y = y,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  move_fn = move_fn
)
plot(result$state_est,
  type = "l", col = "blue", main = "RMPF 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)
}
move_fn <- function(particle, y, sigma) {
  proposal <- particle + rnorm(1, 0, 0.1)
  log_p_curr <- log_likelihood_fn(y = y, particles = particle, sigma = sigma)
  log_p_prop <- log_likelihood_fn(y = y, particles = proposal, sigma = sigma)
  if (log(runif(1)) < (log_p_prop - log_p_curr)) {
    return(proposal)
  } else {
    return(particle)
  }
}

y <- cumsum(rnorm(50))
num_particles <- 100

result <- resample_move_filter(
  y = y,
  num_particles = num_particles,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  move_fn = move_fn,
  mu = 1,
  sigma = 1
)
plot(result$state_est,
  type = "l", col = "blue", main = "RMPF with Parameters",
  ylim = range(c(result$state_est, y))
)
points(y, col = "red", pch = 20)


# With observation gaps
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)
obs_times <- c(1, 2, 3, 5, 6, 7, 8, 9, 10) # skip t=4
data_obs <- data[obs_times]

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)
}
move_fn <- function(particle, y, sigma) {
  proposal <- particle + rnorm(1, 0, 0.1)
  log_p_cur <- log_likelihood_fn(y = y, particles = particle, sigma = sigma)
  log_p_prop <- log_likelihood_fn(y = y, particles = proposal, sigma = sigma)
  if (log(runif(1)) < (log_p_prop - log_p_cur)) {
    return(proposal)
  } else {
    return(particle)
  }
}

result <- resample_move_filter(
  y = data_obs,
  num_particles = 100,
  init_fn = init_fn,
  transition_fn = transition_fn,
  log_likelihood_fn = log_likelihood_fn,
  move_fn = move_fn,
  obs_times = obs_times,
  mu = 1,
  sigma = 1
)
plot(result$state_est,
  type = "l", col = "blue", main = "RMPF with Observation Gaps",
  ylim = range(c(result$state_est, data))
)
points(data_obs, col = "red", pch = 20)

Compute split Rhat statistic

Description

Compute split Rhat statistic

Usage

rhat(chains)

Arguments

Details

Uses the formula for split-Rhat proposed by Gelman et al. (2013).

Value

Rhat value (matrix input) or named vector of Rhat values.

References

Gelman et al. (2013). Bayesian Data Analysis, 3rd Edition.

Examples

# Example with matrix
chains <- matrix(rnorm(3000), nrow = 1000, ncol = 3)
rhat(chains)
#' # Example with data frame
chains_df <- data.frame(
  chain = rep(1:3, each = 1000),
  param1 = rnorm(3000),
  param2 = rnorm(3000)
)
rhat(chains_df)

Summary method for PMMH output

Description

This function returns summary statistics for PMMH output objects, including means, standard deviations, medians, credible intervals, and diagnostics.

Usage

## S3 method for class 'pmmh_output'
summary(object, ...)

Arguments

Value

A data frame containing summary statistics for each parameter.

Examples

# Create dummy chains for two parameters across two chains
chain1 <- data.frame(param1 = rnorm(100), param2 = rnorm(100), chain = 1)
chain2 <- data.frame(param1 = rnorm(100), param2 = rnorm(100), chain = 2)
dummy_output <- list(
  theta_chain = rbind(chain1, chain2),
  diagnostics = list(
    ess = c(param1 = 200, param2 = 190),
    rhat = c(param1 = 1.01, param2 = 1.00)
  )
)
class(dummy_output) <- "pmmh_output"
summary(dummy_output)