R/run_mcmc.R In bssm: Bayesian Inference of Non-Linear and Non-Gaussian State Space Models

Documented in run_mcmcrun_mcmc.gaussianrun_mcmc.nongaussianrun_mcmc.ssm_nlgrun_mcmc.ssm_sde

#' Bayesian Inference of State Space Models
#'
#' 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.
#'
#' @details
#'
#' For linear-Gaussian models, option \code{"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
#' \code{output_type = "full"}.
#'
#' @importFrom stats tsp
#' @param model Model of class \code{bssm_model}.
#' @param iter Positive integer defining the total number of MCMC iterations.
#' @param output_type Either \code{"full"}
#' (default, returns posterior samples from the posterior
#' \eqn{p(\alpha, \theta | y)}), \code{"theta"} (for marginal posterior of
#' theta), or \code{"summary"} (return the mean and variance estimates of the
#' states and posterior samples of theta). See details.
#' @param burnin Positive integer defining the length of the burn-in period
#' which is disregarded from the results. Defaults to \code{iter / 2}.
#' Note that all MCMC algorithms of \code{bssm} use adaptive MCMC during the
#' burn-in period in order to find good proposal distribution.
#' @param thin Positive integer defining the thinning rate. All MCMC algorithms
#' in \code{bssm} use the jump chain representation (see ref [2]), 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 \code{output_type = "summary"},
#' the thinning does not affect the computations of the summary statistics in
#' case of pseudo-marginal methods.
#' @param gamma Tuning parameter for the adaptation of RAM algorithm. Must be
#' between 0 and 1.
#' @param target_acceptance Target acceptance rate for MCMC. Defaults to 0.234.
#' Must be between 0 and 1.
#' @param 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 \eqn{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).
#' @param end_adaptive_phase Logical, if \code{TRUE}, S is held fixed after the
#' burnin period. Default is \code{FALSE}.
#' (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.
#' @param seed Seed for the random number generator (positive integer).
#' @param local_approx If \code{TRUE} (default), Gaussian approximation
#' needed for some of the methods is performed at each iteration.
#' If \code{FALSE}, approximation is updated only once at the start of the
#' MCMC using the initial model.
#' @param max_iter Maximum number of iterations used in Gaussian approximation,
#' as a positive integer.
#' Default is 100 (although typically only few iterations are needed).
#' @param conv_tol Positive tolerance parameter used in Gaussian approximation.
#' @param particles Number of state samples per MCMC iteration for models other
#' than linear-Gaussian models.
#' Ignored if \code{mcmc_type} is \code{"approx"} or \code{"ekf"}.
#' @param mcmc_type What type of MCMC algorithm should be used for models other
#' than linear-Gaussian models? Possible choices are
#' \code{"pm"} for pseudo-marginal MCMC,
#' \code{"da"} for delayed acceptance version of PMCMC ,
#' \code{"approx"} for approximate inference based on the Gaussian
#' approximation of the model,
#' \code{"ekf"} for approximate inference using extended Kalman filter
#' (for \code{ssm_nlg}),
#' or one of the three importance sampling type weighting schemes:
#' \code{"is3"} for simple importance sampling (weight is computed for each
#' MCMC iteration independently),
#' \code{"is2"} for jump chain importance sampling type weighting (default), or
#' \code{"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.
#' @param sampling_method Method for state sampling when for models other than
#' linear-Gaussian models. If \code{"psi"}, \eqn{\psi}-APF is used (default).
#' If \code{"spdk"}, non-sequential importance sampling
#' based on Gaussian approximation is used. If \code{"bsf"}, bootstrap filter
#' is used. If \code{"ekf"}, particle filter based on EKF-proposals are used
#' (only for \code{ssm_nlg} models).
#' @param iekf_iter Non-negative integer. The default zero corresponds to
#' normal EKF, whereas \code{iekf_iter > 0} corresponds to iterated EKF
#' with \code{iekf_iter} iterations. Used only for models of class
#' \code{ssm_nlg}.
#' @param L_c,L_f For \code{ssm_sde} models, Positive integer values defining
#' the discretization levels for first and second stages (defined as 2^L).
#' For pseudo-marginal methods (\code{"pm"}), maximum of these is used.
#' @param ... Ignored.
#' @export
#' @rdname run_mcmc
#' @references
#' [1] Vihola M (2012). Robust adaptive Metropolis algorithm with
#' coerced acceptance rate. Statistics and Computing, 22(5), p 997-1008.
#'
#' [2] 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
#'
#' [3] Helske, J, Vihola, M (2019). bssm: Bayesian Inference of Non-linear and
#' Non-Gaussian State Space Models in R. Arxiv preprint 2101.08492.
#'
run_mcmc <- function(model, ...) {
UseMethod("run_mcmc", model)
}
#' @method run_mcmc gaussian
#' @rdname run_mcmc
#' @export
#' @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)
#'
#' library("dplyr")
#' sumr <- as.data.frame(mcmc_results, variable = "states") %>%
#'   group_by(time) %>%
#'   summarise(mean = mean(value),
#'     lwr = quantile(value, 0.025),
#'     upr = quantile(value, 0.975))
#' library("ggplot2")
#' sumr %>% ggplot(aes(time, mean)) +
#'   geom_ribbon(aes(ymin = lwr, ymax = upr),alpha=0.25) +
#'   geom_line() + theme_bw() +
#'   geom_point(data = data.frame(mean = LakeHuron, time = time(LakeHuron)),
#'     col = 2)
run_mcmc.gaussian <- function(model, iter, output_type = "full",
burnin = floor(iter / 2), thin = 1, gamma = 2 / 3,
seed = sample(.Machine$integer.max, size = 1), ...) { if (!test_flag(end_adaptive_phase)) stop("Argument 'end_adaptive_phase' should be TRUE or FALSE. ") seed <- check_integer(seed, "seed", FALSE, max = .Machine$integer.max)

thin <- check_integer(thin, "thin", max = 100)
iter <- check_integer(iter, "iter", positive = FALSE, max = 1e12)
burnin <- check_integer(burnin, "burnin", max = 1e12)

if (length(model$theta) == 0) stop("No unknown parameters ('model$theta' has length of zero).")
a <- proc.time()

check_prop(target_acceptance)
check_prop(gamma, "gamma")
output_type <- pmatch(output_type, c("full", "summary", "theta"))

if (inherits(model, "bsm_lg")) {
names_ind <- !model$fixed & c(TRUE, TRUE, model$slope, model$seasonal) transformed <- c("sd_y", "sd_level", "sd_slope", "sd_seasonal")[names_ind] model$theta[transformed] <- log(pmax(0.001, model$theta[transformed])) } else { if (inherits(model, "ar1_lg")) { transformed <- c("sigma", "sd_y") model$theta[transformed] <- log(pmax(0.001, model$theta[transformed])) } } if (missing(S)) { S <- diag(0.1 * pmax(0.1, abs(model$theta)), length(model$theta)) } if(output_type == "full") { nsamples <- ifelse(!is.null(nrow(model$y)), nrow(model$y), length(model$y)) *
length(model$a1) * (iter - burnin) / thin * target_acceptance if (nsamples > 1e12) { warning(paste("Number of state samples to be stored is approximately", nsamples, "you might run out of memory.")) } } out <- gaussian_mcmc(model, output_type, iter, burnin, thin, gamma, target_acceptance, S, seed, end_adaptive_phase, threads, model_type(model)) if (output_type == 1) { colnames(out$alpha) <- names(model$a1) } else { if (output_type == 2) { colnames(out$alphahat) <- colnames(out$Vt) <- rownames(out$Vt) <-
names(model$a1) out$alphahat <- ts(out$alphahat, start = start(model$y),
frequency = frequency(model$y)) } } colnames(out$theta) <- rownames(out$S) <- colnames(out$S) <-
names(model$theta) if (inherits(model, "bsm_lg")) { out$theta[, transformed] <- exp(out$theta[, transformed]) } else { if (inherits(model, "ar1_lg")) { out$theta[, transformed] <- exp(out$theta[, transformed]) } } out$call <- match.call()
out$seed <- seed out$iter <- iter
out$burnin <- burnin out$thin <- thin
out$mcmc_type <- "gaussian_mcmc" out$output_type <- output_type
out$time <- proc.time() - a class(out) <- "mcmc_output" attr(out, "model_type") <- class(model)[1] attr(out, "ts") <- list(start = start(model$y), end = end(model$y), frequency = frequency(model$y))
out
}
#' @method run_mcmc nongaussian
#' @rdname run_mcmc
#' @export
#' @examples
#' 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 = 5000,
#'   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 = Hmisc::wtd.mean(value, weight, normwt = TRUE),
#'     lwr = Hmisc::wtd.quantile(value, weight,
#'       0.025, normwt = TRUE),
#'     upr = Hmisc::wtd.quantile(value, weight,
#'       0.975, normwt = TRUE))
#'
#' # 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 = 5000, 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()
#'
run_mcmc.nongaussian <- function(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,
seed = sample(.Machine$integer.max, size = 1), max_iter = 100, conv_tol = 1e-8, ...) { if (!test_flag(end_adaptive_phase)) stop("Argument 'end_adaptive_phase' should be TRUE or FALSE. ") seed <- check_integer(seed, "seed", FALSE, max = .Machine$integer.max)

if (missing(particles)) {
nsim <- eval(match.call(expand.dots = TRUE)$nsim) if (!is.null(nsim)) { warning(paste("Argument nsim is deprecated. Use argument particles", "instead.", sep = " ")) particles <- nsim particles <- check_integer(particles, "particles") } } else { particles <- check_integer(particles, "particles") } threads <- check_integer(threads, "threads") model$max_iter <- check_integer(max_iter, "max_iter", positive = FALSE)
model$conv_tol <- check_positive_real(conv_tol, "conv_tol") thin <- check_integer(thin, "thin", max = 100) iter <- check_integer(iter, "iter", positive = FALSE, max = 1e12) burnin <- check_integer(burnin, "burnin", max = 1e12) if (!test_flag(local_approx)) { stop("Argument 'local_approx' should be TRUE or FALSE. ") } else model$local_approx <- local_approx

if (length(model$theta) == 0) stop("No unknown parameters ('model$theta' has length of zero).")
a <- proc.time()
check_prop(target_acceptance)
check_prop(gamma, "gamma")

output_type <- pmatch(output_type, c("full", "summary", "theta"))
mcmc_type <- match.arg(mcmc_type, c("pm", "da", paste0("is", 1:3), "approx"))
if (mcmc_type == "approx") particles <- 0
if (particles < 2 && mcmc_type != "approx")
stop(paste("Number of state samples less than 2, use 'mcmc_type' 'approx'",

sampling_method <-
pmatch(match.arg(sampling_method, c("psi", "bsf", "spdk")),
c("psi", "bsf", "spdk"))

dists <-
c("svm", "poisson", "binomial", "negative binomial", "gamma", "gaussian")
model$distribution <- pmatch(model$distribution, dists, duplicates.ok = TRUE) - 1

if(output_type == "full") {
nsamples <- ifelse(!is.null(nrow(model$y)), nrow(model$y), length(model$y)) * length(model$a1) * (iter - burnin) / thin * target_acceptance
if (nsamples > 1e12) {
warning(paste("Number of state samples to be stored is approximately",
nsamples, "you might run out of memory."))
}
}

if (inherits(model, "bsm_ng")) {

names_ind <-
c(!model$fixed & c(TRUE, model$slope, model$seasonal), model$noise)

transformed <- c(
c("sd_level", "sd_slope", "sd_seasonal", "sd_noise")[names_ind],
if (dists[model$distribution + 1] %in% dists[4:5]) "phi") model$theta[transformed] <- log(pmax(0.001, model$theta[transformed])) } else { if (inherits(model, "ar1_ng")) { transformed <- c("sigma", if (dists[model$distribution + 1] %in% dists[4:5]) "phi")

model$theta[transformed] <- log(pmax(0.001, model$theta[transformed]))
}
}

if (missing(S)) {
S <- diag(0.1 * pmax(0.1, abs(model$theta)), length(model$theta))
}

switch(mcmc_type,
"da" = {
out <- nongaussian_da_mcmc(model,
output_type, particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method, model_type(model))
},
"pm" = {
out <- nongaussian_pm_mcmc(model, output_type,
particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method, model_type(model))
},
"is1" =,
"is2" =,
"is3" = {
out <- nongaussian_is_mcmc(model, output_type,
particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method,
pmatch(mcmc_type, paste0("is", 1:3)), model_type(model), FALSE)
},
"approx" = {
out <- nongaussian_is_mcmc(model, output_type,
particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method, 2, model_type(model), TRUE)
})
if (output_type == 1) {
colnames(out$alpha) <- names(model$a1)
} else {
if (output_type == 2) {
colnames(out$alphahat) <- colnames(out$Vt) <- rownames(out$Vt) <- names(model$a1)
out$alphahat <- ts(out$alphahat, start = start(model$y), frequency = frequency(model$y))
}
}

colnames(out$theta) <- rownames(out$S) <- colnames(out$S) <- names(model$theta)
if (inherits(model, "bsm_ng")) {
out$theta[, transformed] <- exp(out$theta[, transformed])
} else {
if (inherits(model, "ar1_ng")) {
out$theta[, transformed] <- exp(out$theta[, transformed])
}
}
out$iter <- iter out$burnin <- burnin
out$thin <- thin out$mcmc_type <- mcmc_type
out$output_type <- output_type out$call <- match.call()
out$seed <- seed out$time <- proc.time() - a
class(out) <- "mcmc_output"
attr(out, "model_type") <- class(model)[1]
attr(out, "ts") <-
list(start = start(model$y), end = end(model$y),
frequency = frequency(model$y)) out } #' @method run_mcmc ssm_nlg #' @rdname run_mcmc #' @export run_mcmc.ssm_nlg <- function(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-8, iekf_iter = 0, ...) {

stop("Argument 'end_adaptive_phase' should be TRUE or FALSE. ")

seed <- check_integer(seed, "seed", FALSE, max = .Machine$integer.max) if (missing(particles)) { nsim <- eval(match.call(expand.dots = TRUE)$nsim)
if (!is.null(nsim)) {
warning(paste("Argument nsim is deprecated. Use argument particles",
particles <- nsim
particles <- check_integer(particles, "particles")
}
} else {
particles <- check_integer(particles, "particles")
}

max_iter <- check_integer(max_iter, "max_iter", positive = FALSE)
conv_tol <- check_positive_real(conv_tol, "conv_tol")
thin <- check_integer(thin, "thin", max = 100)
iter <- check_integer(iter, "iter", positive = FALSE, max = 1e12)
burnin <- check_integer(burnin, "burnin", max = 1e12)
iekf_iter <- check_integer(iekf_iter, "iekf_iter", positive = FALSE)

if (length(model$theta) == 0) stop("No unknown parameters ('model$theta' has length of zero).")
a <- proc.time()

check_prop(target_acceptance)
check_prop(gamma, "gamma")

output_type <- pmatch(output_type, c("full", "summary", "theta"))
mcmc_type <- match.arg(mcmc_type, c("pm", "da", paste0("is", 1:3),
"ekf", "approx"))
if (mcmc_type %in% c("ekf", "approx")) particles <- 0
sampling_method <- pmatch(match.arg(sampling_method, c("psi", "bsf", "ekf")),
c("psi", "bsf", NA, "ekf"))

if (missing(S)) {
S <- diag(0.1 * pmax(0.1, abs(model$theta)), length(model$theta))
}

if (particles < 2 && !(mcmc_type %in% c("ekf", "approx")))
stop(paste("Number of state samples less than 2, use 'mcmc_type'",
"'approx' or 'ekf' instead.", sep = " "))

if(output_type == "full") {
nsamples <- ifelse(!is.null(nrow(model$y)), nrow(model$y), length(model$y)) * model$n_states * (iter - burnin) / thin * target_acceptance
if (nsamples > 1e12) {
warning(paste("Number of state samples to be stored is approximately",
nsamples, "you might run out of memory."))
}
}
out <- switch(mcmc_type,
"da" = {
nonlinear_da_mcmc(t(model$y), model$Z, model$H, model$T,
model$R, model$Z_gn, model$T_gn, model$a1, model$P1, model$theta, model$log_prior_pdf, model$known_params,
model$known_tv_params, as.integer(model$time_varying),
model$n_states, model$n_etas, seed,
particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method, iekf_iter, output_type)
},
"pm" = {
nonlinear_pm_mcmc(t(model$y), model$Z, model$H, model$T,
model$R, model$Z_gn, model$T_gn, model$a1, model$P1, model$theta, model$log_prior_pdf, model$known_params,
model$known_tv_params, as.integer(model$time_varying),
model$n_states, model$n_etas, seed,
particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method, iekf_iter, output_type)
},
"is1" =,
"is2" =,
"is3" = {
if (sampling_method == 4)
stop("IS-MCMC with extended particle filter is (not yet) supported.")
nonlinear_is_mcmc(t(model$y), model$Z, model$H, model$T,
model$R, model$Z_gn, model$T_gn, model$a1, model$P1, model$theta, model$log_prior_pdf, model$known_params,
model$known_tv_params, as.integer(model$time_varying),
model$n_states, model$n_etas, seed,
particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method, max_iter, conv_tol, iekf_iter,
output_type, FALSE)
},
"ekf" = {
nonlinear_ekf_mcmc(t(model$y), model$Z, model$H, model$T,
model$R, model$Z_gn, model$T_gn, model$a1, model$P1, model$theta, model$log_prior_pdf, model$known_params,
model$known_tv_params, as.integer(model$time_varying),
model$n_states, model$n_etas, seed,
iter, burnin, thin, gamma, target_acceptance, S,
},
"approx" = {
nonlinear_is_mcmc(t(model$y), model$Z, model$H, model$T,
model$R, model$Z_gn, model$T_gn, model$a1, model$P1, model$theta, model$log_prior_pdf, model$known_params,
model$known_tv_params, as.integer(model$time_varying),
model$n_states, model$n_etas, seed,
particles, iter, burnin, thin, gamma, target_acceptance, S,
sampling_method, max_iter, conv_tol,
iekf_iter, output_type, TRUE)
}
)
if (output_type == 1) {
colnames(out$alpha) <- model$state_names
} else {
if (output_type == 2) {
colnames(out$alphahat) <- colnames(out$Vt) <- rownames(out$Vt) <- model$state_names
out$alphahat <- ts(out$alphahat, start = start(model$y), frequency = frequency(model$y))
}
}

colnames(out$theta) <- rownames(out$S) <- colnames(out$S) <- names(model$theta)

out$iter <- iter out$burnin <- burnin
out$thin <- thin out$mcmc_type <- mcmc_type
out$output_type <- output_type out$call <- match.call()
out$seed <- seed out$time <- proc.time() - a
class(out) <- "mcmc_output"
attr(out, "model_type") <- "ssm_nlg"
attr(out, "ts") <-
list(start = start(model$y), end = end(model$y),
frequency = frequency(model$y)) out } #' @method run_mcmc ssm_sde #' @rdname run_mcmc #' @export #' @references #' 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 run_mcmc.ssm_sde <- function(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), ...) {

if (any(c(model$drift, model$diffusion, model$ddiffusion, model$prior_pdf,
model$obs_pdf) %in% c("<pointer: (nil)>", "<pointer: 0x0>"))) { stop(paste("NULL pointer detected, please recompile the pointer file", "and reconstruct the model.", sep = " ")) } if (!test_flag(end_adaptive_phase)) stop("Argument 'end_adaptive_phase' should be TRUE or FALSE. ") seed <- check_integer(seed, "seed", FALSE, max = .Machine$integer.max)

if (missing(particles)) {
nsim <- eval(match.call(expand.dots = TRUE)$nsim) if (!is.null(nsim)) { warning(paste("Argument nsim is deprecated. Use argument particles", "instead.", sep = " ")) particles <- nsim } } particles <- check_integer(particles, "particles") threads <- check_integer(threads, "threads") thin <- check_integer(thin, "thin", max = 100) iter <- check_integer(iter, "iter", positive = FALSE, max = 1e12) burnin <- check_integer(burnin, "burnin", max = 1e12) if (length(model$theta) == 0)
stop("No unknown parameters ('model$theta' has length of zero).") a <- proc.time() check_prop(target_acceptance) check_prop(gamma, "gamma") output_type <- pmatch(output_type, c("full", "summary", "theta")) mcmc_type <- match.arg(mcmc_type, c("pm", "da", paste0("is", 1:3))) if (missing(S)) { S <- diag(0.1 * pmax(0.1, abs(model$theta)), length(model$theta)) } if (mcmc_type != "pm") { if (L_f <= L_c) stop("L_f should be larger than L_c.") if (L_c < 1) stop("L_c should be at least 1") } else { if (missing(L_c)) L_c <- 0 if (missing(L_f)) L_f <- 0 L <- max(L_c, L_f) if (L <= 0) stop("L should be positive.") } if(output_type == "full") { nsamples <- length(model$y) * (iter - burnin) / thin * target_acceptance
if (nsamples > 1e12) {
warning(paste("Number of state samples to be stored is approximately",
nsamples, "you might run out of memory."))
}
}
out <- switch(mcmc_type,
"da" = {
out <- sde_da_mcmc(model$y, model$x0, model$positive, model$drift, model$diffusion, model$ddiffusion,
model$prior_pdf, model$obs_pdf, model$theta, particles, L_c, L_f, seed, iter, burnin, thin, gamma, target_acceptance, S, end_adaptive_phase, output_type) }, "pm" = { out <- sde_pm_mcmc(model$y, model$x0, model$positive,
model$drift, model$diffusion, model$ddiffusion, model$prior_pdf, model$obs_pdf, model$theta,
particles, L, seed,
iter, burnin, thin, gamma, target_acceptance, S,
},
"is1" =,
"is2" =,
"is3" = {
out <- sde_is_mcmc(model$y, model$x0, model$positive, model$drift, model$diffusion, model$ddiffusion,
model$prior_pdf, model$obs_pdf, model$theta, particles, L_c, L_f, seed, iter, burnin, thin, gamma, target_acceptance, S, end_adaptive_phase, pmatch(mcmc_type, paste0("is", 1:3)), threads, output_type) }) colnames(out$alpha) <- model$state_names colnames(out$theta) <- rownames(out$S) <- colnames(out$S) <-
names(model$theta) out$iter <- iter
out$burnin <- burnin out$thin <- thin
out$mcmc_type <- mcmc_type out$output_type <- output_type
out$call <- match.call() out$seed <- seed
out$time <- proc.time() - a class(out) <- "mcmc_output" attr(out, "model_type") <- "ssm_sde" attr(out, "ts") <- list(start = start(model$y), end = end(model$y), frequency = frequency(model$y))
out
}


Try the bssm package in your browser

Any scripts or data that you put into this service are public.

bssm documentation built on Sept. 21, 2021, 5:09 p.m.