run_mcmc: Bayesian Inference of State Space Models

In bssm: Bayesian Inference of Non-Linear and Non-Gaussian State Space Models

Bayesian Inference of State Space Models

Description

Adaptive Markov chain Monte Carlo simulation for SSMs using Robust Adaptive Metropolis algorithm by Vihola (2012). Several different MCMC sampling schemes are implemented, see parameter arguments, package vignette, Vihola, Helske, Franks (2020) and Helske and Vihola (2021) for details.

Usage

run_mcmc(model, ...)

## S3 method for class 'lineargaussian'
run_mcmc(
  model,
  iter,
  output_type = "full",
  burnin = floor(iter/2),
  thin = 1,
  gamma = 2/3,
  target_acceptance = 0.234,
  S,
  end_adaptive_phase = FALSE,
  threads = 1,
  seed = sample(.Machine$integer.max, size = 1),
  verbose,
  ...
)

## S3 method for class 'nongaussian'
run_mcmc(
  model,
  iter,
  particles,
  output_type = "full",
  mcmc_type = "is2",
  sampling_method = "psi",
  burnin = floor(iter/2),
  thin = 1,
  gamma = 2/3,
  target_acceptance = 0.234,
  S,
  end_adaptive_phase = FALSE,
  local_approx = TRUE,
  threads = 1,
  seed = sample(.Machine$integer.max, size = 1),
  max_iter = 100,
  conv_tol = 1e-08,
  verbose,
  ...
)

## S3 method for class 'ssm_nlg'
run_mcmc(
  model,
  iter,
  particles,
  output_type = "full",
  mcmc_type = "is2",
  sampling_method = "bsf",
  burnin = floor(iter/2),
  thin = 1,
  gamma = 2/3,
  target_acceptance = 0.234,
  S,
  end_adaptive_phase = FALSE,
  threads = 1,
  seed = sample(.Machine$integer.max, size = 1),
  max_iter = 100,
  conv_tol = 1e-08,
  iekf_iter = 0,
  verbose,
  ...
)

## S3 method for class 'ssm_sde'
run_mcmc(
  model,
  iter,
  particles,
  output_type = "full",
  mcmc_type = "is2",
  L_c,
  L_f,
  burnin = floor(iter/2),
  thin = 1,
  gamma = 2/3,
  target_acceptance = 0.234,
  S,
  end_adaptive_phase = FALSE,
  threads = 1,
  seed = sample(.Machine$integer.max, size = 1),
  verbose,
  ...
)

Arguments

model	Model of class bssm_model.
...	Ignored.
iter	A positive integer defining the total number of MCMC iterations. Suitable value depends on the model, data, and the choice of specific algorithms (mcmc_type and sampling_method). As increasing iter also increases run time, it is is generally good idea to first test the performance with a small values, e.g., less than 10000.
output_type	Either "full" (default, returns posterior samples from the posterior p(\alpha, \theta | y)), "theta" (for marginal posterior of theta), or "summary" (return the mean and variance estimates of the states and posterior samples of theta). See details.
burnin	A positive integer defining the length of the burn-in period which is disregarded from the results. Defaults to iter / 2. Note that all MCMC algorithms of bssm use adaptive MCMC during the burn-in period in order to find good proposal distribution.
thin	A positive integer defining the thinning rate. All the MCMC algorithms in bssm use the jump chain representation (see refs), and the thinning is applied to these blocks. Defaults to 1. For IS-corrected methods, larger value can also be statistically more effective. Note: With output_type = "summary", the thinning does not affect the computations of the summary statistics in case of pseudo-marginal methods.
gamma	Tuning parameter for the adaptation of RAM algorithm. Must be between 0 and 1.
target_acceptance	Target acceptance rate for MCMC. Defaults to 0.234. Must be between 0 and 1.
S	Matrix defining the initial value for the lower triangular matrix of the RAM algorithm, so that the covariance matrix of the Gaussian proposal distribution is SS'. Note that for some parameters (currently the standard deviation, dispersion, and autoregressive parameters of the BSM and AR(1) models) the sampling is done in unconstrained parameter space, i.e. internal_theta = log(theta) (and logit(rho) or AR coefficient).
end_adaptive_phase	Logical, if TRUE, S is held fixed after the burnin period. Default is FALSE.
threads	Number of threads for state simulation. Positive integer (default is 1). Note that parallel computing is only used in the post-correction phase of IS-MCMC and when sampling the states in case of (approximate) Gaussian models.
seed	Seed for the C++ RNG (positive integer).
verbose	If TRUE, prints a progress bar to the console. If missing, defined by rlang::is_interactive. Set to FALSE if number of iterations is less than 50.
particles	A positive integer defining the number of state samples per MCMC iteration for models other than linear-Gaussian models. Ignored if mcmc_type is "approx" or "ekf". Suitable values depend on the model, the data, mcmc_type and sampling_method. While large values provide more accurate estimates, the run time also increases with respect to to the number of particles, so it is generally a good idea to test the run time firstwith a small number of particles, e.g., less than 100.
mcmc_type	What type of MCMC algorithm should be used for models other than linear-Gaussian models? Possible choices are "pm" for pseudo-marginal MCMC, "da" for delayed acceptance version of PMCMC , "approx" for approximate inference based on the Gaussian approximation of the model, "ekf" for approximate inference using extended Kalman filter (for ssm_nlg), or one of the three importance sampling type weighting schemes: "is3" for simple importance sampling (weight is computed for each MCMC iteration independently), "is2" for jump chain importance sampling type weighting (default), or "is1" for importance sampling type weighting where the number of particles used for weight computations is proportional to the length of the jump chain block.
sampling_method	Method for state sampling when for models other than linear-Gaussian models. If "psi", \psi-APF is used (default). If "spdk", non-sequential importance sampling based on Gaussian approximation is used. If "bsf", bootstrap filter is used. If "ekf", particle filter based on EKF-proposals are used (only for ssm_nlg models).
local_approx	If TRUE (default), Gaussian approximation needed for some of the methods is performed at each iteration. If FALSE, approximation is updated only once at the start of the MCMC using the initial model.
max_iter	Maximum number of iterations used in Gaussian approximation, as a positive integer. Default is 100 (although typically only few iterations are needed).
conv_tol	Positive tolerance parameter used in Gaussian approximation.
iekf_iter	Non-negative integer. The default zero corresponds to normal EKF, whereas iekf_iter > 0 corresponds to iterated EKF with iekf_iter iterations. Used only for models of class ssm_nlg.
L_c, L_f	For ssm_sde models, Positive integer values defining the discretization levels for first and second stages (defined as 2^L). For pseudo-marginal methods ("pm"), maximum of these is used.

Details

For linear-Gaussian models, option "summary" does not simulate states directly but computes the posterior means and variances of states using fast Kalman smoothing. This is slightly faster, more memory efficient and more accurate than calculations based on simulation smoother. In other cases, the means and covariances are computed using the full output of particle filter instead of subsampling one of these as in case of output_type = "full". The states are sampled up to the time point n+1 where n is the length of the input time series i.e. the last values are one-step-ahead predictions. (for predicting further, see ?predict.mcmc_output).

Initial values for the sampling are taken from the model object (model$theta). If you want to continue from previous run, you can reconstruct your original model by plugging in the previously obtained parameters to model$theta, providing the S matrix for the RAM algorithm and setting burnin = 0. See example. Note however, that this is not identical as running all the iterations once, due to the RNG "discontinuity" and because even without burnin bssm does include "theta_0" i.e. the initial theta in the final chain (even with burnin=0).

Value

An object of class mcmc_output.

References

Vihola M (2012). Robust adaptive Metropolis algorithm with coerced acceptance rate. Statistics and Computing, 22(5), p 997-1008. https://doi.org/10.1007/s11222-011-9269-5

Vihola, M, Helske, J, Franks, J (2020). Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo. Scand J Statist. 1-38. https://doi.org/10.1111/sjos.12492

Helske J, Vihola M (2021). bssm: Bayesian Inference of Non-linear and Non-Gaussian State Space Models in R. The R Journal (2021) 13:2, 578-589. https://doi.org/10.32614/RJ-2021-103

Vihola, M, Helske, J, Franks, J. Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo. Scand J Statist. 2020; 1-38. https://doi.org/10.1111/sjos.12492

Examples

model <- ar1_lg(LakeHuron, rho = uniform(0.5,-1,1),
  sigma = halfnormal(1, 10), mu = normal(500, 500, 500),
  sd_y = halfnormal(1, 10))

mcmc_results <- run_mcmc(model, iter = 2e4)
summary(mcmc_results, return_se = TRUE)

sumr <- summary(mcmc_results, variable = "states")
library("ggplot2")
ggplot(sumr, aes(time, Mean)) +
  geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), alpha = 0.25) +
  geom_line() + theme_bw() +
  geom_point(data = data.frame(Mean = LakeHuron, time = time(LakeHuron)),
    col = 2)

# Continue from the previous run
model$theta[] <- mcmc_results$theta[nrow(mcmc_results$theta), ]
run_more <- run_mcmc(model, S = mcmc_results$S, iter = 1000, burnin = 0)

set.seed(1)
n <- 50
slope <- cumsum(c(0, rnorm(n - 1, sd = 0.001)))
level <- cumsum(slope + c(0, rnorm(n - 1, sd = 0.2)))
y <- rpois(n, exp(level))
poisson_model <- bsm_ng(y,
  sd_level = halfnormal(0.01, 1),
  sd_slope = halfnormal(0.01, 0.1),
  P1 = diag(c(10, 0.1)), distribution = "poisson")

# Note small number of iterations for CRAN checks
mcmc_out <- run_mcmc(poisson_model, iter = 1000, particles = 10,
  mcmc_type = "da")
summary(mcmc_out, what = "theta", return_se = TRUE)

set.seed(123)
n <- 50
sd_level <- 0.1
drift <- 0.01
beta <- -0.9
phi <- 5

level <- cumsum(c(5, drift + rnorm(n - 1, sd = sd_level)))
x <- 3 + (1:n) * drift + sin(1:n + runif(n, -1, 1))
y <- rnbinom(n, size = phi, mu = exp(beta * x + level))

model <- bsm_ng(y, xreg = x,
  beta = normal(0, 0, 10),
  phi = halfnormal(1, 10),
  sd_level = halfnormal(0.1, 1),
  sd_slope = halfnormal(0.01, 0.1),
  a1 = c(0, 0), P1 = diag(c(10, 0.1)^2),
  distribution = "negative binomial")

# run IS-MCMC
# Note small number of iterations for CRAN checks
fit <- run_mcmc(model, iter = 4000,
  particles = 10, mcmc_type = "is2", seed = 1)

# extract states
d_states <- as.data.frame(fit, variable = "states", time = 1:n)

library("dplyr")
library("ggplot2")

 # compute summary statistics
level_sumr <- d_states |>
  filter(variable == "level") |>
  group_by(time) |>
  summarise(mean = diagis::weighted_mean(value, weight),
    lwr = diagis::weighted_quantile(value, weight,
      0.025),
    upr = diagis::weighted_quantile(value, weight,
      0.975))

# visualize
level_sumr |> ggplot(aes(x = time, y = mean)) +
  geom_line() +
  geom_line(aes(y = lwr), linetype = "dashed", na.rm = TRUE) +
  geom_line(aes(y = upr), linetype = "dashed", na.rm = TRUE) +
  theme_bw() +
  theme(legend.title = element_blank()) +
  xlab("Time") + ylab("Level")

# theta
d_theta <- as.data.frame(fit, variable = "theta")
ggplot(d_theta, aes(x = value)) +
 geom_density(aes(weight = weight), adjust = 2, fill = "#92f0a8") +
 facet_wrap(~ variable, scales = "free") +
 theme_bw()


# Bivariate Poisson model:

set.seed(1)
x <- cumsum(c(3, rnorm(19, sd = 0.5)))
y <- cbind(
  rpois(20, exp(x)),
  rpois(20, exp(x)))

prior_fn <- function(theta) {
  # half-normal prior using transformation
  dnorm(exp(theta), 0, 1, log = TRUE) + theta # plus jacobian term
}

update_fn <- function(theta) {
  list(R = array(exp(theta), c(1, 1, 1)))
}

model <- ssm_mng(y = y, Z = matrix(1,2,1), T = 1,
  R = 0.1, P1 = 1, distribution = "poisson",
  init_theta = log(0.1),
  prior_fn = prior_fn, update_fn = update_fn)

# Note small number of iterations for CRAN checks
out <- run_mcmc(model, iter = 4000, mcmc_type = "approx")

sumr <- as.data.frame(out, variable = "states") |>
  group_by(time) |> mutate(value = exp(value)) |>
  summarise(mean = mean(value),
    ymin = quantile(value, 0.05), ymax = quantile(value, 0.95))
ggplot(sumr, aes(time, mean)) +
geom_ribbon(aes(ymin = ymin, ymax = ymax),alpha = 0.25) +
geom_line() +
geom_line(data = data.frame(mean = y[, 1], time = 1:20),
  colour = "tomato") +
geom_line(data = data.frame(mean = y[, 2], time = 1:20),
  colour = "tomato") +
theme_bw()

bssm documentation built on Nov. 2, 2023, 6:25 p.m.