Nothing
R/exdqlmMCMC.R
In exdqlm: Extended Dynamic Quantile Linear Models
Defines functions
.exdqlm_mcmc_multichain_diag
exdqlmMCMC
.exdqlm_mcmc_warmup_state
.exdqlm_theta_state_mcmc_controls
.exdqlm_dqlm_sigma_mcmc_controls
.exdqlm_latent_state_mcmc_controls
Documented in
exdqlmMCMC
#' exDQLM - MCMC algorithm
#'
#' The function applies a Markov chain Monte Carlo (MCMC) algorithm to sample the posterior of an exDQLM.
#'
#' @param y A univariate time-series.
#' @param p0 The quantile of interest, a value between 0 and 1.
#' @param model List of the state-space model including \code{GG}, \code{FF},
#' prior parameters \code{m0} and \code{C0}.
#' @param df Discount factors for each block.
#' @param dim.df Dimension of each block of discount factors.
#' @param fix.gamma Logical value indicating whether to fix gamma at
#' \code{gam.init}. Default is \code{FALSE}.
#' @param gam.init Initial value for gamma (skewness parameter), or value at
#' which gamma will be fixed if \code{fix.gamma = TRUE}.
#' @param fix.sigma Logical value indicating whether to fix sigma at `sig.init`.
#' Default is `FALSE`.
#' @param sig.init Initial value for sigma (scale parameter), or value at which
#' sigma will be fixed if \code{fix.sigma = TRUE}.
#' @param dqlm.ind Logical value indicating whether to fix gamma at \code{0},
#' reducing the exDQLM to the AL/DQLM special case. Default is \code{FALSE}.
#' @param Sig.mh Covariance matrix used in the random walk MH step to jointly sample sigma and gamma.
#' @param joint.sample Logical value indicating whether or not to recompute `Sig.mh` based off the initial burn-in samples of gamma and sigma. Default is `FALSE`.
#' @param n.burn Number of MCMC iterations to burn. Default is `n.burn = 2000`.
#' @param n.mcmc Number of MCMC iterations to sample. Default is `n.mcmc = 1500`.
#' @param init.from.isvb Logical value indicating whether to use the legacy ISVB
#' warm start when `init.from.vb = TRUE`. Default is `FALSE`, which favors
#' `LDVB` as the default VB warm start. This flag only chooses the
#' warm-start source; it does not change the subsequent MCMC proposal kernel.
#' @param init.from.vb Optional logical. If `TRUE`, run a VB pre-initialization
#' step (`LDVB` by default, or `ISVB` when `init.from.isvb = TRUE`) and
#' initialize MCMC from converged VB moments. Default is `TRUE`. If
#' explicitly set to `NULL`, it falls back to `init.from.isvb` behavior for
#' backward compatibility.
#' @param vb_init_controls Optional list controlling VB warm start. Supported keys:
#' `method` (`"isvb"` or `"ldvb"`), `tol`, `n.IS`, `n.samp`, `max_iter`, `verbose`.
#' @param vb_init_fit Optional precomputed VB fit object. If supplied, warm start
#' uses this object directly and does not rerun VB internally.
#' @param mcmc_control Optional normalized MCMC control list, usually from
#' [exal_make_mcmc_control()]. When supplied, the core MCMC arguments and
#' warmup blocks are read from `mcmc_control` first and then merged with the
#' explicit function arguments. When omitted, exAL-style dynamic MCMC uses the
#' package's conservative default `(sigma, gamma)` warmup profile
#' automatically; explicit controls remain the advanced override path.
#' @param sigmagam_controls Optional list controlling warmup/freeze for the
#' exDQLM sigma/gamma block during MCMC.
#' @param latent_state_controls Optional list controlling early latent-state
#' warmup/freeze in dynamic MCMC. Supported keys include
#' `freeze_burnin_iters`, `freeze_only_during_burn`, `force_after_warmup`,
#' and `mode` (`"u_only"` or `"u_st_pair"`).
#' @param theta_state_controls Optional list controlling early theta-state
#' warmup/freeze in dynamic MCMC. Supported keys include
#' `freeze_burnin_iters`, `freeze_only_during_burn`, and
#' `force_after_warmup`.
#' @param dqlm_sigma_controls Optional list controlling sigma-only
#' warmup/freeze in the DQLM branch. Supported keys mirror
#' `sigmagam_controls`.
#' @param mh.proposal Character; proposal kernel for the exDQLM scale/skew block.
#' `"slice"` (default) uses
#' an exact sigma GIG update plus a bounded univariate slice sampler directly
#' on `gamma`; `"laplace_rw"` uses a Laplace-informed covariance then RW;
#' and `"rw"` uses joint random-walk MH on `(log sigma, logit gamma)`.
#' This choice is separate from the VB warm-start method.
#' @param mh.adapt Logical; adapt MH proposal scale during burn-in.
#' @param mh.adapt.interval Integer; adaptation interval (iterations).
#' @param mh.target.accept Numeric length-2 vector with lower/upper target acceptance rates.
#' @param mh.scale.bounds Numeric length-2 vector with min/max global scaling for MH covariance.
#' @param mh.max_scale.step Numeric in (0,1); maximum fractional scale change per adaptation step.
#' @param mh.min_burn_adapt Minimum burn-in iterations required to enable adaptation.
#' @param slice.width Positive numeric width for the bounded slice sampler when
#' `mh.proposal = "slice"`. Default `0.1` for parity with `bqrgal`.
#' @param slice.max.steps Positive integer or `Inf`; maximum stepping-out
#' expansions for the slice sampler.
#' @param trace.diagnostics Logical; if `TRUE`, retain per-iteration
#' sigma/gamma/s/u diagnostics under `mh.diagnostics$trace`. Set `FALSE` for
#' lighter-weight runs.
#' @param trace.every Positive integer; when `trace.diagnostics = TRUE`, record
#' one diagnostics row every `trace.every` iterations.
#' @param verbose.every Positive integer controlling how often console progress
#' is printed when `verbose = TRUE`. Default `500`, independent of
#' `trace.every`.
#' @param progress_callback Optional callback invoked with a named list at MCMC
#' start, at each progress checkpoint, and on completion. Intended for
#' workflow-level progress logging.
#' @param PriorSigma List of parameters for inverse gamma prior on sigma; shape
#' \code{a_sig} and scale \code{b_sig}. Default is an inverse gamma with mean
#' 1, or \code{sig.init} when supplied, and variance 10.
#' @param PriorGamma List of parameters for truncated Student-t prior on gamma;
#' center \code{m_gam}, scale \code{s_gam}, and degrees of freedom
#' \code{df_gam}. Default is a standard Student-t with 1 degree of freedom,
#' truncated to the support of gamma.
#' @param verbose Logical value indicating whether progress should be displayed.
#'
#' @return An object with classes "\code{exdqlmMCMC}" and "\code{exdqlmFit}"
#' containing the following:
#' \itemize{
#' \item `y` - Time-series data used to fit the model.
#' \item `run.time` - Algorithm run time in seconds.
#' \item `dqlm.ind` - Logical value indicating whether gamma was fixed at `0`, reducing the exDQLM to the special case of the DQLM.
#' \item `model` - List of the state-space model including `GG`, `FF`, prior parameters `m0` and `C0`.
#' \item `p0` - The quantile which was estimated.
#' \item `df` - Discount factors used for each block.
#' \item `dim.df` - Dimension used for each block of discount factors.
#' \item `samp.theta` - Posterior sample of the state vector.
#' \item `samp.post.pred` - Sample of the posterior predictive distributions.
#' \item `map.standard.forecast.errors` - MAP standardized one-step-ahead forecast errors.
#' \item `samp.sigma` - Posterior sample of scale parameter sigma.
#' \item `samp.vts` - Posterior sample of latent parameters, v_t.
#' \item `theta.out` - List containing the distributions of the state vector including filtered distribution parameters (`fm` and `fC`) and smoothed distribution parameters (`sm` and `sC`).
#' \item `n.burn` Number of MCMC iterations that were burned.
#' \item `n.mcmc` Number of MCMC iterations that were sampled.
#' }
#' If `dqlm.ind=FALSE`, the object also contains the following:
#' \itemize{
#' \item `samp.gamma` - Posterior sample of skewness parameter gamma.
#' \item `samp.sts` - Posterior sample of latent parameters, s_t.
#' \item `init.log.sigma` - Burned samples of log sigma from the random walk MH joint sampling of sigma and gamma.
#' \item `init.logit.gamma` - Burned samples of logit gamma from the random walk MH joint sampling of sigma and gamma.
#' \item `accept.rate` - Acceptance rate of the MH step.
#' \item `accept.rate.burn` - MH acceptance rate during burn-in.
#' \item `accept.rate.keep` - MH acceptance rate in kept MCMC samples.
#' \item `Sig.mh` - Covariance matrix used in MH step to jointly sample sigma and gamma.
#' \item `mh.diagnostics` - MH tuning diagnostics (proposal mode, scaling path, adaptation summary).
#' \item `diagnostics` - ESS and chain-ready summaries for sigma/gamma.
#' }
#' @export
#'
#' @examples
#' \donttest{
#' data("scIVTmag", package = "exdqlm")
#' y = scIVTmag[1:80]
#' trend.comp = polytrendMod(order = 1, m0 = stats::quantile(y, 0.85), C0 = 10)
#' seas.comp = seasMod(p = 365, h = c(1,2), C0 = 10*diag(4))
#' model = trend.comp + seas.comp
#' M2 = exdqlmMCMC(y, p0=0.85, model, df = c(1,1), dim.df = c(1,4),
#' gam.init = -3.5, sig.init = 15,
#' n.burn = 40, n.mcmc = 40,
#' init.from.vb = FALSE, verbose = FALSE)
#'
#' M2_al = exdqlmMCMC(y, p0=0.85, model, df = c(1,1), dim.df = c(1,4),
#' dqlm.ind = TRUE, sig.init = 15,
#' n.burn = 30, n.mcmc = 30,
#' init.from.vb = FALSE, verbose = FALSE)
#' }
#'
#' @name exdqlmMCMC
NULL
.exdqlm_latent_state_mcmc_controls <- function(latent_cfg = NULL, default_mode = c("u_only", "u_st_pair")) {
default_mode <- match.arg(default_mode)
latent_cfg <- latent_cfg %||% list()
freeze_burnin_iters <- suppressWarnings(as.integer(
latent_cfg$freeze_burnin_iters %||%
latent_cfg$freeze_latent_burnin_iters %||%
latent_cfg$freeze_ust_burnin_iters %||%
0L
)[1L])
if (!is.finite(freeze_burnin_iters) || freeze_burnin_iters < 0L) freeze_burnin_iters <- 0L
mode <- tolower(as.character(
latent_cfg$mode %||%
latent_cfg$latent_mode %||%
latent_cfg$freeze_mode %||%
default_mode
)[1L])
if (!(mode %in% c("u_only", "u_st_pair"))) mode <- default_mode
min_postwarmup_updates <- suppressWarnings(as.integer(
latent_cfg$min_postwarmup_updates %||%
latent_cfg$latent_min_postwarmup_updates %||%
0L
)[1L])
if (!is.finite(min_postwarmup_updates) || min_postwarmup_updates < 0L) {
min_postwarmup_updates <- 0L
}
list(
mode = mode,
freeze_burnin_iters = freeze_burnin_iters,
freeze_only_during_burn = if (is.null(latent_cfg$freeze_only_during_burn)) TRUE else isTRUE(latent_cfg$freeze_only_during_burn),
force_after_warmup = if (is.null(latent_cfg$force_after_warmup)) TRUE else isTRUE(latent_cfg$force_after_warmup),
min_postwarmup_updates = min_postwarmup_updates,
trace = if (is.null(latent_cfg$trace)) TRUE else isTRUE(latent_cfg$trace)
)
}
.exdqlm_dqlm_sigma_mcmc_controls <- function(sigma_cfg = NULL) {
sigma_cfg <- sigma_cfg %||% list()
freeze_burnin_iters <- suppressWarnings(as.integer(
sigma_cfg$freeze_burnin_iters %||%
sigma_cfg$freeze_sigma_burnin_iters %||%
sigma_cfg$freeze_dqlm_sigma_burnin_iters %||%
0L
)[1L])
if (!is.finite(freeze_burnin_iters) || freeze_burnin_iters < 0L) freeze_burnin_iters <- 0L
list(
freeze_burnin_iters = freeze_burnin_iters,
freeze_only_during_burn = if (is.null(sigma_cfg$freeze_only_during_burn)) TRUE else isTRUE(sigma_cfg$freeze_only_during_burn),
force_after_warmup = if (is.null(sigma_cfg$force_after_warmup)) TRUE else isTRUE(sigma_cfg$force_after_warmup),
trace = if (is.null(sigma_cfg$trace)) TRUE else isTRUE(sigma_cfg$trace)
)
}
.exdqlm_theta_state_mcmc_controls <- function(theta_cfg = NULL) {
theta_cfg <- theta_cfg %||% list()
freeze_burnin_iters <- suppressWarnings(as.integer(
theta_cfg$freeze_burnin_iters %||%
theta_cfg$freeze_theta_burnin_iters %||%
theta_cfg$freeze_theta_state_burnin_iters %||%
0L
)[1L])
if (!is.finite(freeze_burnin_iters) || freeze_burnin_iters < 0L) freeze_burnin_iters <- 0L
list(
freeze_burnin_iters = freeze_burnin_iters,
freeze_only_during_burn = if (is.null(theta_cfg$freeze_only_during_burn)) TRUE else isTRUE(theta_cfg$freeze_only_during_burn),
force_after_warmup = if (is.null(theta_cfg$force_after_warmup)) TRUE else isTRUE(theta_cfg$force_after_warmup)
)
}
.exdqlm_mcmc_warmup_state <- function(iter, ctrl, postwarmup_update_count, n_burn) {
active <- isTRUE(
ctrl$freeze_burnin_iters > 0L &&
iter <= ctrl$freeze_burnin_iters &&
(!ctrl$freeze_only_during_burn || iter <= n_burn)
)
force_now <- isTRUE(
!active &&
ctrl$freeze_burnin_iters > 0L &&
ctrl$force_after_warmup &&
postwarmup_update_count <= 0L
)
reason <- if (isTRUE(active)) {
"warmup"
} else if (isTRUE(force_now)) {
"force_after_warmup"
} else {
"scheduled"
}
list(active = active, force_now = force_now, reason = reason)
}
exdqlmMCMC <- function(y,p0,model,df,dim.df,fix.gamma=FALSE,gam.init=NA,fix.sigma=FALSE,sig.init=NA,dqlm.ind=FALSE,
Sig.mh,joint.sample=FALSE,n.burn=2000,n.mcmc=1500,init.from.isvb=FALSE,PriorSigma=NULL,PriorGamma=NULL,verbose=TRUE,
init.from.vb=TRUE,vb_init_controls=NULL,vb_init_fit=NULL,mcmc_control=NULL,sigmagam_controls=NULL,latent_state_controls=NULL,theta_state_controls=NULL,dqlm_sigma_controls=NULL,
mh.proposal=c("slice","laplace_rw","rw"),mh.adapt=TRUE,mh.adapt.interval=50L,
mh.target.accept=c(0.20,0.45),mh.scale.bounds=c(0.1,10),
mh.max_scale.step=0.35,mh.min_burn_adapt=50L,
slice.width=0.1,slice.max.steps=Inf,
trace.diagnostics=TRUE,trace.every=1L,verbose.every=500L,
progress_callback=NULL){
if (!is.null(mcmc_control)) {
mcmc_control <- exal_make_mcmc_control(control = mcmc_control)
if (!is.null(mcmc_control$n_burn)) n.burn <- as.integer(mcmc_control$n_burn)[1L]
if (!is.null(mcmc_control$n_mcmc)) n.mcmc <- as.integer(mcmc_control$n_mcmc)[1L]
if (!is.null(mcmc_control$verbose)) verbose <- isTRUE(mcmc_control$verbose)
if (!is.null(mcmc_control$progress_every)) verbose.every <- as.integer(mcmc_control$progress_every)[1L]
if (!is.null(mcmc_control$init_from_vb)) init.from.vb <- isTRUE(mcmc_control$init_from_vb)
if (!is.null(mcmc_control$vb_warm_start_control)) {
vb_init_controls <- utils::modifyList(vb_init_controls %||% list(), mcmc_control$vb_warm_start_control)
}
if (!is.null(mcmc_control$sigmagam)) sigmagam_controls <- utils::modifyList(sigmagam_controls %||% list(), mcmc_control$sigmagam)
if (!is.null(mcmc_control$latent_state)) latent_state_controls <- utils::modifyList(latent_state_controls %||% list(), mcmc_control$latent_state)
if (!is.null(mcmc_control$theta)) theta_state_controls <- utils::modifyList(theta_state_controls %||% list(), mcmc_control$theta)
if (!is.null(mcmc_control$dqlm_sigma)) dqlm_sigma_controls <- utils::modifyList(dqlm_sigma_controls %||% list(), mcmc_control$dqlm_sigma)
}
# check inputs
y = check_ts(y)
model = check_mod(model)
rv = check_logics(gam.init,sig.init,fix.gamma,fix.sigma,dqlm.ind)
gam.init = rv$gam.init
dqlm.ind = rv$dqlm.ind
fix.gamma = rv$fix.gamma
### MCMC iterations
if(n.mcmc<=0){
stop("number of mcmc samples must be positive")
}
if(verbose & n.burn<=0){
warning("mcmc will be sampled without burn-in, a burn-in is recommended even if initializing using the isvb algorithm")
n.burn=0
}
I = n.mcmc + n.burn
mh.proposal <- match.arg(mh.proposal)
if (is.null(init.from.vb)) {
init.from.vb <- isTRUE(init.from.isvb)
}
if (!is.null(vb_init_fit)) {
init.from.vb <- TRUE
}
init.from.vb <- isTRUE(init.from.vb)
vb.ctrl.default <- list(
method = if (isTRUE(init.from.isvb)) "isvb" else "ldvb",
tol = 0.5,
n.IS = 200L,
n.samp = 200L,
max_iter = getOption("exdqlm.max_iter", 200L),
verbose = FALSE,
ld_controls = NULL
)
if (is.null(vb_init_controls)) vb_init_controls <- list()
vb.ctrl <- utils::modifyList(vb.ctrl.default, vb_init_controls)
vb.ctrl$method <- tolower(as.character(vb.ctrl$method)[1])
if (!(vb.ctrl$method %in% c("isvb", "ldvb"))) vb.ctrl$method <- "isvb"
vb.ctrl$tol <- as.numeric(vb.ctrl$tol)[1]
if (!is.finite(vb.ctrl$tol) || vb.ctrl$tol <= 0) vb.ctrl$tol <- 0.5
vb.ctrl$n.IS <- suppressWarnings(as.integer(vb.ctrl$n.IS)[1])
if (!is.finite(vb.ctrl$n.IS) || vb.ctrl$n.IS < 20L) vb.ctrl$n.IS <- 200L
vb.ctrl$n.samp <- suppressWarnings(as.integer(vb.ctrl$n.samp)[1])
if (!is.finite(vb.ctrl$n.samp) || vb.ctrl$n.samp < 20L) vb.ctrl$n.samp <- 200L
vb.ctrl$max_iter <- suppressWarnings(as.integer(vb.ctrl$max_iter)[1])
if (!is.finite(vb.ctrl$max_iter) || vb.ctrl$max_iter < 5L) vb.ctrl$max_iter <- 200L
vb.ctrl$verbose <- isTRUE(vb.ctrl$verbose)
if (!is.null(vb.ctrl$ld_controls) && !is.list(vb.ctrl$ld_controls)) {
stop("vb_init_controls$ld_controls must be a list or NULL")
}
mh.adapt <- isTRUE(mh.adapt)
mh.adapt.interval <- suppressWarnings(as.integer(mh.adapt.interval)[1])
if (!is.finite(mh.adapt.interval) || mh.adapt.interval < 5L) mh.adapt.interval <- 50L
mh.min_burn_adapt <- suppressWarnings(as.integer(mh.min_burn_adapt)[1])
if (!is.finite(mh.min_burn_adapt) || mh.min_burn_adapt < 20L) mh.min_burn_adapt <- 50L
if (length(mh.target.accept) != 2L) mh.target.accept <- c(0.20, 0.45)
mh.target.accept <- as.numeric(mh.target.accept)
mh.target.accept <- sort(pmin(pmax(mh.target.accept, 0.01), 0.99))
if (length(mh.scale.bounds) != 2L) mh.scale.bounds <- c(0.1, 10)
mh.scale.bounds <- sort(as.numeric(mh.scale.bounds))
if (!all(is.finite(mh.scale.bounds)) || mh.scale.bounds[1] <= 0 || mh.scale.bounds[2] <= mh.scale.bounds[1]) {
mh.scale.bounds <- c(0.1, 10)
}
mh.max_scale.step <- as.numeric(mh.max_scale.step)[1]
if (!is.finite(mh.max_scale.step) || mh.max_scale.step <= 0 || mh.max_scale.step >= 1) {
mh.max_scale.step <- 0.35
}
mh.laplace.refresh.interval <- suppressWarnings(as.integer(getOption("exdqlm.mcmc.laplace_refresh_interval", mh.adapt.interval))[1])
if (!is.finite(mh.laplace.refresh.interval) || mh.laplace.refresh.interval < 5L) {
mh.laplace.refresh.interval <- mh.adapt.interval
}
mh.laplace.refresh.start <- suppressWarnings(as.integer(getOption("exdqlm.mcmc.laplace_refresh_start", mh.min_burn_adapt))[1])
if (!is.finite(mh.laplace.refresh.start) || mh.laplace.refresh.start < 1L) {
mh.laplace.refresh.start <- mh.min_burn_adapt
}
mh.laplace.refresh.weight <- as.numeric(getOption("exdqlm.mcmc.laplace_refresh_weight", 0.60))[1]
if (!is.finite(mh.laplace.refresh.weight) || mh.laplace.refresh.weight <= 0 || mh.laplace.refresh.weight > 1) {
mh.laplace.refresh.weight <- 0.60
}
sigmagam_ctrl <- .exal_sigmagam_mcmc_controls(sigmagam_controls)
sigmagam_ctrl <- .exal_clamp_mcmc_sigmagam_control(sigmagam_ctrl, n_burn = n.burn)
latent_ctrl <- .exdqlm_latent_state_mcmc_controls(
latent_state_controls,
default_mode = if (isTRUE(dqlm.ind)) "u_only" else "u_st_pair"
)
theta_ctrl <- .exdqlm_theta_state_mcmc_controls(theta_state_controls)
dqlm_sigma_ctrl <- .exdqlm_dqlm_sigma_mcmc_controls(dqlm_sigma_controls)
slice.width <- as.numeric(slice.width)[1]
if (!is.finite(slice.width) || slice.width <= 0) slice.width <- 0.1
slice.max.steps <- as.numeric(slice.max.steps)[1]
if (!(is.infinite(slice.max.steps) || (is.finite(slice.max.steps) && slice.max.steps >= 1 && floor(slice.max.steps) == slice.max.steps))) {
slice.max.steps <- Inf
}
trace.diagnostics <- isTRUE(trace.diagnostics)
trace.every <- suppressWarnings(as.integer(trace.every)[1])
if (!is.finite(trace.every) || trace.every < 1L) trace.every <- 1L
verbose.every <- suppressWarnings(as.integer(verbose.every)[1])
if (!is.finite(verbose.every) || verbose.every < 1L) verbose.every <- 500L
verbose_every_env <- suppressWarnings(as.integer(Sys.getenv("EXDQLM_MCMC_PROGRESS_EVERY", NA_character_))[1])
if (is.finite(verbose_every_env) && !is.na(verbose_every_env) && verbose_every_env >= 1L) {
verbose.every <- verbose_every_env
}
safe_progress_callback <- function(info) {
if (!is.function(progress_callback)) return(invisible(NULL))
try(progress_callback(info), silent = TRUE)
invisible(NULL)
}
state_signal <- function(FF_local, theta_mat) {
drop(colSums(FF_local * theta_mat))
}
warmup_state <- function(iter, ctrl, postwarmup_update_count) {
.exdqlm_mcmc_warmup_state(
iter = iter,
ctrl = ctrl,
postwarmup_update_count = postwarmup_update_count,
n_burn = n.burn
)
}
max_abs_finite <- function(x) {
x <- as.numeric(x)
if (!length(x) || any(!is.finite(x))) return("NA")
format(max(abs(x), na.rm = TRUE), digits = 6)
}
count_nonfinite <- function(x) {
sum(!is.finite(as.numeric(x)))
}
count_nonpositive <- function(x) {
x <- as.numeric(x)
sum(is.finite(x) & x <= 0)
}
if (n.burn < mh.min_burn_adapt) mh.adapt <- FALSE
### Define L and U
bounds = .gamma_bounds(p0)
L = bounds["L"]; U = bounds["U"]
if(!is.na(gam.init)){
if(gam.init < L | gam.init > U){
stop(sprintf("gam.init must be between %s and %s for %s quantile",round(L,3),round(U,3),p0))
}
}
### sigma and gamma priors
# sigma ~ IG(a_sig,b_sig)
if(is.null(PriorSigma)){
m_sigma = 1
v_sigma = 10
PriorSigma$a_sig = (m_sigma^2)/(v_sigma) + 2
PriorSigma$b_sig = (m_sigma^3)/(v_sigma) + m_sigma
}else{
if(!is.list(PriorSigma) | any( is.na( match(c("a_sig", "b_sig"),names(PriorSigma)) ) )){
stop("`PriorSigma` must be a list containing `a_sig` and `b_sig`")
}
}
PriorSigmaDens<-function(sigma){ LaplacesDemon::dinvgamma(sigma,shape=PriorSigma$a_sig,scale=PriorSigma$b_sig) }
# gamma ~ truncated student t on L,U
PriorGamma <- .normalize_gamma_prior_trunc_t(PriorGamma)
PriorGammaDens <- function(gamma) {
.gamma_prior_density_trunc_t(gamma, bounds = c(L, U), PriorGamma = PriorGamma, log = FALSE)
}
### state-space model
## prior, theta ~ N(m0,C0)
m0 = model$m0
C0 = model$C0
#
TT = length(y)
p = length(m0)
if(!is.na(dim(model$GG)[3])){
if(dim(model$GG)[3] != TT){stop("time-varying dimension of GG does not match length of y")}
}
GG = array(model$GG,c(p,p,TT)); model$GG = GG
if(ncol(model$FF)>1){
if(ncol(model$FF) != TT){stop("time-varying dimension of FF does not match length of y")}
}
FF = matrix(model$FF,p,TT); model$FF = FF
## discount factor blocking
if(!methods::hasArg(dim.df)){
if(length(df)!=1){
stop("length of component discount factors does not match length of component dimensions")
}
dim.df = p
}
df.mat = make_df_mat(df,dim.df,p)
### backend controls (MCMC-specific)
use_cpp_mcmc_opt <- isTRUE(getOption("exdqlm.use_cpp_mcmc", FALSE))
cpp_mcmc_mode <- tolower(as.character(getOption("exdqlm.cpp_mcmc_mode", "strict")))
if (!(cpp_mcmc_mode %in% c("strict", "fast"))) {
warning("Invalid exdqlm.cpp_mcmc_mode; using 'strict'.")
cpp_mcmc_mode <- "strict"
}
has_cpp_mcmc <- exists("mcmc_ffbs_smooth_cpp", mode = "function") &&
exists("mcmc_ffbs_sample_cpp", mode = "function")
if (use_cpp_mcmc_opt && !has_cpp_mcmc) {
warning("exdqlm.use_cpp_mcmc=TRUE but C++ MCMC FFBS kernels not available; using R backend.")
}
# strict mode keeps R kernels to preserve exact legacy path; fast enables C++ FFBS.
use_cpp_mcmc <- isTRUE(use_cpp_mcmc_opt && has_cpp_mcmc && identical(cpp_mcmc_mode, "fast"))
mcmc_backend <- if (use_cpp_mcmc) "C++" else "R"
if (verbose) {
cat(sprintf("MCMC backend: %s (mode=%s)\n", mcmc_backend, cpp_mcmc_mode))
}
cpp_ffbs_smooth <- function(ex.f, ex.q) {
out <- mcmc_ffbs_smooth_cpp(
GG = GG,
m0 = as.numeric(m0),
C0 = C0,
FF = FF,
y = as.numeric(y),
ex_f = as.numeric(ex.f),
ex_q = as.numeric(ex.q),
df_mat = df.mat
)
out$standard.forecast.errors <- as.numeric(out$standard.forecast.errors)
out$sm <- as.matrix(out$sm)
out$fm <- as.matrix(out$fm)
out$sC <- array(out$sC, dim = c(p, p, TT))
out$fC <- array(out$fC, dim = c(p, p, TT))
out
}
cpp_ffbs_sample <- function(ex.f, ex.q) {
out <- mcmc_ffbs_sample_cpp(
GG = GG,
m0 = as.numeric(m0),
C0 = C0,
FF = FF,
y = as.numeric(y),
ex_f = as.numeric(ex.f),
ex_q = as.numeric(ex.q),
df_mat = df.mat
)
out$standard.forecast.errors <- as.numeric(out$standard.forecast.errors)
out$sam.theta <- as.matrix(out$sam.theta)
out$fm <- as.matrix(out$fm)
out$fC <- array(out$fC, dim = c(p, p, TT))
out
}
# function to produce smoothed estimates for return value
smoothed_theta<-function(ex.f,ex.q){
# initialize ffbs
m <- sm <- matrix(NA,p,TT)
C <- sC <- array(NA,c(p,p,TT))
standard.forecast.errors <- rep(NA,TT)
## forward filter
# first iteration
a = as.vector(GG[,,1]%*%m0)
P = .exdqlm_regularize_cov(GG[,,1]%*%C0%*%t(GG[,,1]), context = "mcmc_smooth_P_t1")
R = .exdqlm_regularize_cov(P + df.mat*P, context = "mcmc_smooth_R_t1")
f = t(FF[,1])%*%a + ex.f[1]
q = .exdqlm_regularize_var(t(FF[,1])%*%R%*%FF[,1] + ex.q[1], context = "mcmc_smooth_q_t1")
m[,1] = a + t(R)%*%FF[,1]%*%(y[1]-f)/q[1]
C[,,1] = .exdqlm_regularize_cov(
R - t(R)%*%FF[,1]%*%t(FF[,1])%*%R/q[1],
context = "mcmc_smooth_C_t1"
)
standard.forecast.errors[1] = (y[1]-f)/sqrt(q)
# t = 2:TT
for(t in 2:TT){
a = as.vector(GG[,,t]%*%m[,(t-1)])
P = .exdqlm_regularize_cov(GG[,,t]%*%C[,,(t-1)]%*%t(GG[,,t]), context = sprintf("mcmc_smooth_P_t%d", t))
R = .exdqlm_regularize_cov(P + df.mat*P, context = sprintf("mcmc_smooth_R_t%d", t))
f = t(FF[,t])%*%a + ex.f[t]
fB = t(FF[,t])%*%R
q = .exdqlm_regularize_var(fB%*%FF[,t] + ex.q[t], context = sprintf("mcmc_smooth_q_t%d", t))
m[,t] = a + t(fB)%*%(y[t]-f)/q[1]
C[,,t] = .exdqlm_regularize_cov(
R - t(fB)%*%fB/q[1],
context = sprintf("mcmc_smooth_C_t%d", t)
)
standard.forecast.errors[t] = (y[t]-f)/sqrt(q)
}
## backwards smoothing
sC[,,TT] = C[,,TT]
sm[,TT] = m[,TT]
for(t in (TT-1):1){
P = .exdqlm_regularize_cov(GG[,,(t+1)]%*%C[,,(t)]%*%t(GG[,,(t+1)]), context = sprintf("mcmc_smooth_back_P_t%d", t + 1L))
R.info = .exdqlm_cov_inverse(P + df.mat*P, context = sprintf("mcmc_smooth_back_R_t%d", t + 1L))
sB = C[,,t]%*%t(GG[,,(t+1)])%*%R.info$inverse
sm[,t] = m[,t] + sB%*%(sm[,(t+1)]-as.vector(GG[,,(t+1)]%*%m[,(t)]))
sC[,,t] = .exdqlm_regularize_cov(
C[,,t] + sB%*%(sC[,,(t+1)]-R.info$Sigma)%*%t(sB),
context = sprintf("mcmc_smooth_back_C_t%d", t)
)
}
return(list(standard.forecast.errors=standard.forecast.errors,sm=sm,sC=sC,fm=m,fC=C))
}
if (use_cpp_mcmc) {
smoothed_theta <- function(ex.f, ex.q) cpp_ffbs_smooth(ex.f, ex.q)
}
### Initialize MCMC
init.log.sigma <- init.logit.gamma <- rep(NA,n.burn)
save.sigma <- save.gamma <- rep(NA,n.mcmc)
save.Ut <- save.st <- matrix(NA,TT,n.mcmc)
save.theta <- array(NA,c(p,TT,n.mcmc))
save.post.pred <- matrix(NA,TT,n.mcmc)
vb.out <- NULL
gig_backend <- "cpp_devroye_required"
gig_eps <- 1e-12
gig_b_floor <- .gig_b_floor()
current_iter <- NA_integer_
sample_gig_cpp_required <- function(chi, psi, lambda = 0.5, context = "gig") {
if (!exists("sample_gig_devroye_vector", mode = "function")) {
stop(sprintf("%s requires sample_gig_devroye_vector(), but it is not available", context))
}
chi <- as.numeric(chi)
psi <- as.numeric(psi)[1]
lambda <- as.numeric(lambda)[1]
iter_suffix <- if (is.finite(current_iter)) sprintf(" (iter=%d)", current_iter) else ""
bad <- which(!is.finite(chi))
if (length(bad)) {
stop(sprintf("%s%s chi has %d non-finite values (first index=%d)", context, iter_suffix, length(bad), bad[1]))
}
badneg <- which(chi < 0)
if (length(badneg)) {
stop(sprintf("%s%s chi has %d negative values (first index=%d, value=%.6g)", context, iter_suffix, length(badneg), badneg[1], chi[badneg[1]]))
}
if (!is.finite(psi)) {
stop(sprintf("%s%s psi must be finite; got %.6g", context, iter_suffix, psi))
}
if (!is.finite(lambda)) {
stop(sprintf("%s%s lambda must be finite; got %.6g", context, iter_suffix, lambda))
}
chi <- pmax(chi, gig_b_floor)
psi <- max(psi, gig_eps)
draws <- as.numeric(sample_gig_devroye_vector(
1L, p = lambda, a = psi, b_vec = chi
)[1, ])
bad_draws <- which(!is.finite(draws) | draws <= 0)
if (length(bad_draws)) {
stop(sprintf("%s%s sample_gig_devroye_vector returned %d invalid draws (first index=%d, value=%.6g)",
context, iter_suffix, length(bad_draws), bad_draws[1], draws[bad_draws[1]]))
}
pmax(draws, gig_eps)
}
run_vb_init <- function() {
old_opt <- options(exdqlm.max_iter = vb.ctrl$max_iter)
on.exit(options(old_opt), add = TRUE)
if (vb.ctrl$method == "ldvb") {
if (is.list(vb.ctrl$ld_controls) && length(vb.ctrl$ld_controls)) {
ld_opt_names <- paste0("exdqlm.dynamic.ldvb.", names(vb.ctrl$ld_controls))
old_ld_opt <- options(stats::setNames(vb.ctrl$ld_controls, ld_opt_names))
on.exit(options(old_ld_opt), add = TRUE)
}
exdqlmLDVB(
y = y, p0 = p0, model = model, df = df, dim.df = dim.df,
fix.gamma = fix.gamma, gam.init = gam.init,
fix.sigma = fix.sigma, sig.init = sig.init,
dqlm.ind = dqlm.ind,
tol = vb.ctrl$tol, n.samp = vb.ctrl$n.samp,
PriorSigma = PriorSigma, PriorGamma = PriorGamma,
verbose = vb.ctrl$verbose
)
} else {
exdqlmISVB(
y = y, p0 = p0, model = model, df = df, dim.df = dim.df,
fix.gamma = fix.gamma, gam.init = gam.init,
fix.sigma = fix.sigma, sig.init = sig.init,
dqlm.ind = dqlm.ind,
tol = vb.ctrl$tol, n.IS = vb.ctrl$n.IS, n.samp = vb.ctrl$n.samp,
PriorSigma = PriorSigma, PriorGamma = PriorGamma,
verbose = vb.ctrl$verbose
)
}
}
# Set initial values
if(init.from.vb){
if(verbose){
cat(sprintf("running %s algorithm to initialize mcmc\n", toupper(vb.ctrl$method)))
}
if (!is.null(vb_init_fit)) {
vb.out <- vb_init_fit
if (verbose) {
cat("using provided vb_init_fit object for MCMC initialization\n")
}
} else {
vb.out <- run_vb_init()
}
cursam.sigma <- ifelse(fix.sigma,sig.init,ifelse(dqlm.ind,vb.out$sig.out$E.sigma,vb.out$gammasig.out$E.sigma))
cursam.Ut <- vb.out$vts.out$E.uts
cursam.theta <- vb.out$theta.out$sm
}else{
cursam.sigma <- m_sigma
cursam.Ut <- rep(1/m_sigma,TT)
cursam.theta <- matrix(m0,p,TT)
}
if (verbose) {
cat("GIG backend: C++ Devroye (required)\n")
}
######## exDQLM
if(!dqlm.ind){
### Define logit and inverse logit functions
logit = function(x){log((x-L)/(U-x))}
inv.logit = function(x){(U*exp(x)+L)/(exp(x)+1)}
log_prior_gamma <- function(gamma) {
.gamma_log_prior_trunc_t(gamma, bounds = c(L, U), PriorGamma = PriorGamma)
}
### Additional initial values
if(!is.null(vb.out)){
cursam.st <- vb.out$sts.out$E.sts
cursam.gamma <- ifelse(fix.gamma,gam.init,vb.out$gammasig.out$E.gam)
cursam.logit.gamma <- logit(cursam.gamma)
cursam.log.sigma <- log(cursam.sigma)
}else{
cursam.st <- truncnorm::rtruncnorm(TT,a=0,b=Inf,mean=0,sd=1)
cursam.gamma <- ifelse(!is.na(gam.init),gam.init,(L+U)/2)
cursam.logit.gamma <- logit(cursam.gamma)
cursam.log.sigma <- log(cursam.sigma)
}
### Initialize MH
n.accept = 0
n.accept.burn = 0L
n.accept.keep = 0L
n.trial.burn = 0L
n.trial.keep = 0L
adapt.history <- data.frame(
iter = integer(0),
window_accept = numeric(0),
mh_scale = numeric(0),
sig11 = numeric(0),
sig22 = numeric(0),
laplace_refreshed = logical(0),
stringsAsFactors = FALSE
)
trace_rows <- if (trace.diagnostics) vector("list", ceiling(I / trace.every)) else NULL
trace_idx <- 0L
mh.scale <- 1
window.accept <- 0L
window.total <- 0L
laplace_refresh_attempts <- 0L
laplace_refresh_success <- 0L
prep_Sig_mh <- function(S) {
S <- suppressWarnings(as.matrix(S))
if (!all(dim(S) == c(2L, 2L))) {
S <- diag(c(ifelse(fix.sigma, 0, 0.005), ifelse(fix.gamma, 0, 0.005)))
}
S[!is.finite(S)] <- 0
S <- (S + t(S)) / 2
if (fix.sigma) {
S[1, ] <- 0
S[, 1] <- 0
}
if (fix.gamma) {
S[2, ] <- 0
S[, 2] <- 0
}
for (j in 1:2) {
if (!is.finite(S[j, j]) || S[j, j] < 0) S[j, j] <- 0
}
if (!fix.sigma && S[1, 1] <= 0) S[1, 1] <- 0.005
if (!fix.gamma && S[2, 2] <= 0) S[2, 2] <- 0.005
S
}
build_chol <- function(S) {
S <- prep_Sig_mh(S)
if (fix.gamma || fix.sigma) {
sqrt(S)
} else {
out <- tryCatch(t(chol(S)), error = function(e) NULL)
if (is.null(out)) {
eig <- eigen(S, symmetric = TRUE)
vals <- pmax(eig$values, 1e-8)
out <- eig$vectors %*% diag(sqrt(vals), 2, 2) %*% t(eig$vectors)
}
out
}
}
if(!methods::hasArg(Sig.mh)){
if(!is.null(vb.out)){
sig.samples <- NULL
gam.samples <- NULL
if (!is.null(vb.out$gammasig.out$sigma.samples) && !is.null(vb.out$gammasig.out$gamma.samples)) {
sig.samples <- as.numeric(vb.out$gammasig.out$sigma.samples)
gam.samples <- as.numeric(vb.out$gammasig.out$gamma.samples)
} else if (!is.null(vb.out$samp.sigma) && !is.null(vb.out$samp.gamma)) {
sig.samples <- as.numeric(vb.out$samp.sigma)
gam.samples <- as.numeric(vb.out$samp.gamma)
}
if (!is.null(sig.samples) && !is.null(gam.samples) &&
all(is.finite(sig.samples)) && all(sig.samples > 0) &&
all(is.finite(gam.samples)) && all(gam.samples > L) && all(gam.samples < U)) {
Sig.mh <- stats::cov(cbind(log(sig.samples), logit(gam.samples)))
} else {
Sig.mh <- diag(c(ifelse(fix.sigma,0,0.005),ifelse(fix.gamma,0,0.005)))
}
}else{
Sig.mh = diag(c(ifelse(fix.sigma,0,0.005),ifelse(fix.gamma,0,0.005)))
}
}
Sig.mh <- prep_Sig_mh(Sig.mh)
Sig.mh.initial <- Sig.mh
chol_Sig.mh <- build_chol(Sig.mh)
# exdqlm function sample theta ffbs
ex_samp_theta<-function(ex.f,ex.q,gamma,sigma,sts,tau,c_tau){
# initialize ffbs
m <- sam.theta <- matrix(NA,p,TT)
C <- array(NA,c(p,p,TT))
standard.forecast.errors <- post.pred <- rep(NA,TT)
## forward filter
# first iteration
a = as.vector(GG[,,1]%*%m0)
P = .exdqlm_regularize_cov(GG[,,1]%*%C0%*%t(GG[,,1]), context = "mcmc_dqlm_sample_P_t1")
R = .exdqlm_regularize_cov(P + df.mat*P, context = "mcmc_dqlm_sample_R_t1")
f = t(FF[,1])%*%a + ex.f[1]
q = .exdqlm_regularize_var(t(FF[,1])%*%R%*%FF[,1] + ex.q[1], context = "mcmc_dqlm_sample_q_t1")
m[,1] = a + t(R)%*%FF[,1]%*%(y[1]-f)/q[1]
C[,,1] = .exdqlm_regularize_cov(
R - t(R)%*%FF[,1]%*%t(FF[,1])%*%R/q[1],
context = "mcmc_dqlm_sample_C_t1"
)
standard.forecast.errors[1] = (y[1]-f)/sqrt(q)
# t = 2:TT
for(t in 2:TT){
a = as.vector(GG[,,t]%*%m[,(t-1)])
P = .exdqlm_regularize_cov(GG[,,t]%*%C[,,(t-1)]%*%t(GG[,,t]), context = sprintf("mcmc_dqlm_sample_P_t%d", t))
R = .exdqlm_regularize_cov(P + df.mat*P, context = sprintf("mcmc_dqlm_sample_R_t%d", t))
f = t(FF[,t])%*%a + ex.f[t]
fB = t(FF[,t])%*%R
q = .exdqlm_regularize_var(fB%*%FF[,t] + ex.q[t], context = sprintf("mcmc_dqlm_sample_q_t%d", t))
m[,t] = a + t(fB)%*%(y[t]-f)/q[1]
C[,,t] = .exdqlm_regularize_cov(
R - t(fB)%*%fB/q[1],
context = sprintf("mcmc_dqlm_sample_C_t%d", t)
)
standard.forecast.errors[t] = (y[t]-f)/sqrt(q)
}
## backwards sample
sC_TT = .exdqlm_regularize_cov(C[,,TT], context = "mcmc_dqlm_sample_sC_TT")
svd.sC = svd(sC_TT)
sam.theta[,TT] = m[,TT] + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta <- numeric(TT)
reg_theta[TT] <- drop(crossprod(FF[,TT], sam.theta[,TT]))
post.pred[TT] = rexal(1,tau,reg_theta[TT]+c_tau*sigma*abs(gamma)*sts[TT],sigma,0)
for(t in (TT-1):1){
P = .exdqlm_regularize_cov(GG[,,(t+1)]%*%C[,,(t)]%*%t(GG[,,(t+1)]), context = sprintf("mcmc_dqlm_back_P_t%d", t + 1L))
R.info = .exdqlm_cov_inverse(P + df.mat*P, context = sprintf("mcmc_dqlm_back_R_t%d", t + 1L))
sB = C[,,t]%*%t(GG[,,(t+1)])%*%R.info$inverse
sm = m[,t] + sB%*%(sam.theta[,(t+1)]-as.vector(GG[,,(t+1)]%*%m[,(t)]))
sC = .exdqlm_regularize_cov(
C[,,t] - sB%*%GG[,,(t+1)]%*%C[,,t],
context = sprintf("mcmc_dqlm_back_sC_t%d", t)
)
svd.sC = svd(sC)
sam.theta[,t] = sm + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta[t] <- drop(crossprod(FF[,t], sam.theta[,t]))
post.pred[t] = rexal(1,tau,reg_theta[t]+c_tau*sigma*abs(gamma)*sts[t],sigma,0)
}
return(list(standard.forecast.errors=standard.forecast.errors,post.pred=post.pred,sam.theta=sam.theta,fm=m,fC=C))
}
if (use_cpp_mcmc) {
ex_samp_theta <- function(ex.f, ex.q, gamma, sigma, sts, tau, c_tau) {
out <- cpp_ffbs_sample(ex.f, ex.q)
sam.theta <- out$sam.theta
reg_theta <- state_signal(FF, sam.theta)
post.pred <- vapply(seq_len(TT), function(t) {
rexal(1, tau,
reg_theta[t] + c_tau * sigma * abs(gamma) * sts[t],
sigma, 0)
}, numeric(1))
list(
standard.forecast.errors = out$standard.forecast.errors,
post.pred = post.pred,
sam.theta = sam.theta,
fm = out$fm,
fC = out$fC
)
}
}
# exdqlm function sample uts
ex_samp_uts<-function(reg1,gamma,sigma,sts,a_tau,b_tau,c_tau){
chi <- as.numeric(((y-reg1-sigma*c_tau*abs(gamma)*sts)^2)/(b_tau*sigma))
psi <- (a_tau^2)/(b_tau*sigma) + (2/sigma)
sample_gig_cpp_required(chi = chi, psi = psi, lambda = 0.5, context = "exdqlm_mcmc_uts")
}
# exdqlm function sample sts
ex_samp_sts<-function(reg1,gamma,sigma,uts,a_tau,b_tau,c_tau){
s.sig2<-1/(1+c_tau^2*abs(gamma)^2*sigma/(b_tau*uts))
s.sig2<-pmax(s.sig2, 1e-12)
s.mu<-s.sig2*c_tau*abs(gamma)*(y-(reg1+a_tau*uts))/(b_tau*uts)
truncnorm::rtruncnorm(TT,rep(0,TT),rep(Inf,TT),s.mu,sqrt(s.sig2))
}
# exdqlm function sample sigma and gamma
logL<-function(reg1,log.sigma,logit.gamma,sts,uts){
sigma=exp(log.sigma); gamma=inv.logit(logit.gamma)
temp.p = p.fn(p0,gamma)
a = (1-2*temp.p)/(temp.p*(1-temp.p))
b = (2)/(temp.p*(1-temp.p))
c = (as.numeric(gamma>0)-temp.p)^(-1)
logJ<-logit.gamma-2*log(1+exp(logit.gamma))+log.sigma
PriorGamma<- PriorGammaDens(gamma)
PriorSigma<- PriorSigmaDens(sigma)
sum(stats::dnorm(y,reg1+sigma*c*abs(gamma)*sts+a*uts,sqrt(sigma*b*uts),log = TRUE)) +
sum(stats::dexp(uts,rate = 1/sigma,log=TRUE)) +
log(PriorSigma) + log(PriorGamma) + logJ
}
make_logpost_gamma <- function(reg1, sigma, sts, uts) {
sigma <- as.numeric(sigma)[1]
reg1 <- as.numeric(reg1)
sts <- as.numeric(sts)
uts <- as.numeric(uts)
valid_inputs <- is.finite(sigma) && sigma > 0 &&
all(is.finite(reg1)) && all(is.finite(sts)) && all(is.finite(uts)) &&
all(uts > 0)
if (!valid_inputs) {
return(function(gamma) -Inf)
}
y_center <- y - reg1
sigma_sts <- sigma * sts
sqrt_sigma_uts <- sqrt(sigma * uts)
function(gamma) {
gamma <- as.numeric(gamma)[1]
if (!is.finite(gamma) || gamma <= L || gamma >= U) return(-Inf)
temp.p <- p.fn(p0, gamma)
a <- (1 - 2 * temp.p) / (temp.p * (1 - temp.p))
b <- 2 / (temp.p * (1 - temp.p))
c <- (as.numeric(gamma > 0) - temp.p)^(-1)
if (!all(is.finite(c(a, b, c))) || b <= 0) return(-Inf)
mu_shift <- c * abs(gamma) * sigma_sts + a * uts
ll <- sum(stats::dnorm(
y_center,
mean = mu_shift,
sd = sqrt(b) * sqrt_sigma_uts,
log = TRUE
))
lp <- log_prior_gamma(gamma)
if (!is.finite(ll) || !is.finite(lp)) return(-Inf)
ll + lp
}
}
samp_sigma_exact <- function(reg1, sigma, gamma, sts, uts) {
temp.p <- p.fn(p0, gamma)
a <- (1 - 2 * temp.p) / (temp.p * (1 - temp.p))
b <- 2 / (temp.p * (1 - temp.p))
c <- (as.numeric(gamma > 0) - temp.p)^(-1)
r <- y - reg1 - a * uts
chi_sigma <- sum((r * r) / (b * uts)) + 2 * sum(uts) + 2 * PriorSigma$b_sig
psi_sigma <- ((c * abs(gamma))^2 / b) * sum((sts * sts) / uts)
k_sigma <- -(PriorSigma$a_sig + 1.5 * TT)
sigma_new <- sample_gig_cpp_required(
chi = chi_sigma,
psi = psi_sigma,
lambda = k_sigma,
context = "exdqlm_mcmc_sigma"
)[1]
if (is.finite(sigma_new) && sigma_new > 0) sigma_new else sigma
}
laplace_cov_init <- function(reg1, log.sigma, logit.gamma, sts, uts) {
fn <- function(z) {
val <- logL(reg1, z[1], z[2], sts, uts)
if (is.finite(val)) -val else 1e12
}
opt <- tryCatch(
stats::optim(c(log.sigma, logit.gamma), fn = fn, method = "BFGS",
control = list(maxit = 100), hessian = TRUE),
error = function(e) NULL
)
H <- if (!is.null(opt) && !is.null(opt$hessian)) opt$hessian else NULL
if (is.null(H) || any(!is.finite(H))) {
H <- tryCatch(numDeriv::hessian(fn, x = c(log.sigma, logit.gamma)), error = function(e) NULL)
}
if (is.null(H) || any(!is.finite(H))) return(NULL)
H <- (H + t(H)) / 2
eig <- eigen(H, symmetric = TRUE)
vals <- pmax(eig$values, 1e-6)
cov <- eig$vectors %*% diag(1 / vals, 2, 2) %*% t(eig$vectors)
cov
}
if (identical(mh.proposal, "laplace_rw") && !fix.gamma && !fix.sigma) {
reg1.init <- state_signal(FF, cursam.theta)
cov.lap <- laplace_cov_init(reg1.init, cursam.log.sigma, cursam.logit.gamma, cursam.st, cursam.Ut)
if (!is.null(cov.lap)) {
Sig.mh <- prep_Sig_mh(cov.lap)
}
}
chol_Sig.mh <- build_chol(Sig.mh * (mh.scale^2))
ex_samp_lsiglgam<-function(reg1,log.sigma,logit.gamma,sts,uts,chol_Sig){
prop<-c(log.sigma,logit.gamma)+chol_Sig%*%stats::rnorm(2)
if(inv.logit(prop[2]) < U && inv.logit(prop[2]) > L){
logr<-logL(reg1,prop[1],prop[2],sts,uts)-logL(reg1,log.sigma,logit.gamma,sts,uts)
accept=(log(stats::runif(1))<logr)
}else{
accept=FALSE
}
log.sigma.new<-accept*prop[1]+(1-accept)*log.sigma
logit.gamma.new<-accept*prop[2]+(1-accept)*logit.gamma
return(list(log.sigma=log.sigma.new,logit.gamma=logit.gamma.new,accept=accept))
}
callback.every <- if (trace.diagnostics) {
trace.every
} else {
100L
}
callback.every <- max(1L, as.integer(callback.every)[1])
sigmagam_frozen_trace <- rep(FALSE, I)
sigmagam_forced_postwarmup_trace <- rep(FALSE, I)
sigmagam_update_performed_trace <- rep(FALSE, I)
sigmagam_update_reason_trace <- rep(NA_character_, I)
sigmagam_update_count_trace <- integer(I)
sigmagam_first_active_iter <- NA_integer_
sigmagam_update_count <- 0L
sigmagam_postwarmup_update_count <- 0L
sigmagam_updates_burn <- 0L
sigmagam_updates_keep <- 0L
latent_frozen_trace <- rep(FALSE, I)
latent_forced_postwarmup_trace <- rep(FALSE, I)
latent_update_performed_trace <- rep(FALSE, I)
latent_update_reason_trace <- rep(NA_character_, I)
latent_update_count_trace <- integer(I)
latent_first_active_iter <- NA_integer_
latent_update_count <- 0L
latent_postwarmup_update_count <- 0L
latent_updates_burn <- 0L
latent_updates_keep <- 0L
theta_frozen_trace <- rep(FALSE, I)
theta_forced_postwarmup_trace <- rep(FALSE, I)
theta_update_performed_trace <- rep(FALSE, I)
theta_update_reason_trace <- rep(NA_character_, I)
theta_update_count_trace <- integer(I)
theta_first_active_iter <- NA_integer_
theta_update_count <- 0L
theta_postwarmup_update_count <- 0L
theta_updates_burn <- 0L
theta_updates_keep <- 0L
# Sample from exdqlm posterior
tictoc::tic()
.exdqlm_progress(
"MCMC start",
burn = n.burn,
keep = n.mcmc,
kernel = mh.proposal,
warm_start = if (isTRUE(init.from.vb)) vb.ctrl$method else "none",
.verbose = verbose
)
safe_progress_callback(list(
event = "start",
iter = 0L,
total_iter = as.integer(I),
phase = "burn",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = cursam.gamma,
kernel = mh.proposal,
accept = if (identical(mh.proposal, "slice")) NA_real_ else 0
))
for (i in 1:I){
current_iter <- as.integer(i)
# counter
if (i %% verbose.every == 0L) {
kept_now <- if (i <= n.burn) 0L else (i - n.burn)
.exdqlm_progress(
"MCMC progress",
model = "exDQLM",
phase = if (i <= n.burn) "burn" else "keep",
iter = sprintf("%d/%d", i, I),
accept = if (identical(mh.proposal, "slice")) NULL else n.accept / i,
kept = sprintf("%d/%d", kept_now, n.mcmc),
.verbose = verbose
)
}
if (i %% callback.every == 0L) {
safe_progress_callback(list(
event = "progress",
iter = as.integer(i),
total_iter = as.integer(I),
phase = if (i <= n.burn) "burn" else "keep",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = cursam.gamma,
kernel = mh.proposal,
accept = if (identical(mh.proposal, "slice")) NA_real_ else n.accept / i
))
}
# exAL parameters
tau = p.fn(p0,cursam.gamma)
a_tau = (1-2*tau)/(tau*(1-tau))
b_tau = (2)/(tau*(1-tau))
c_tau = (as.numeric(cursam.gamma>0)-tau)^(-1)
# sample theta
ex.f = cursam.sigma*c_tau*abs(cursam.gamma)*cursam.st + cursam.Ut*a_tau
ex.q = b_tau*cursam.Ut*cursam.sigma
theta_state_sched <- warmup_state(i, theta_ctrl, theta_postwarmup_update_count)
if (isTRUE(theta_state_sched$active)) {
if (any(!is.finite(cursam.theta))) {
stop(sprintf(
"exdqlm_mcmc_theta_warmup_state (iter=%d) invalid warmup theta state: nonfinite=%d",
i, sum(!is.finite(cursam.theta))
))
}
reg_theta_hold <- state_signal(FF, cursam.theta)
theta.out <- list(
standard.forecast.errors = rep(NA_real_, TT),
post.pred = vapply(seq_len(TT), function(t) {
rexal(
1, tau,
reg_theta_hold[t] + c_tau * cursam.sigma * abs(cursam.gamma) * cursam.st[t],
cursam.sigma, 0
)
}, numeric(1)),
sam.theta = cursam.theta,
fm = matrix(NA_real_, p, TT),
fC = array(NA_real_, c(p, p, TT))
)
} else {
theta.out <- ex_samp_theta(ex.f,ex.q,cursam.gamma,cursam.sigma,cursam.st,tau,c_tau)
cursam.theta = theta.out$sam.theta
}
theta_frozen_trace[i] <- isTRUE(theta_state_sched$active)
theta_update_reason_trace[i] <- theta_state_sched$reason
if (!isTRUE(theta_state_sched$active)) {
theta_update_performed_trace[i] <- TRUE
theta_update_count <- theta_update_count + 1L
theta_update_count_trace[i] <- theta_update_count
if (is.na(theta_first_active_iter)) theta_first_active_iter <- as.integer(i)
if (i <= n.burn) {
theta_updates_burn <- theta_updates_burn + 1L
} else {
theta_updates_keep <- theta_updates_keep + 1L
}
if (theta_ctrl$freeze_burnin_iters > 0L) {
theta_postwarmup_update_count <- theta_postwarmup_update_count + 1L
}
theta_forced_postwarmup_trace[i] <- isTRUE(theta_state_sched$force_now)
} else {
theta_update_count_trace[i] <- theta_update_count
theta_forced_postwarmup_trace[i] <- FALSE
}
# sample uts, sts
reg1 = state_signal(FF, cursam.theta)
latent_state_sched <- warmup_state(i, latent_ctrl, latent_postwarmup_update_count)
if (!is.finite(cursam.sigma) || cursam.sigma <= 0 || !is.finite(cursam.gamma) ||
cursam.gamma <= L || cursam.gamma >= U || any(!is.finite(reg1)) || any(!is.finite(cursam.theta))) {
stop(sprintf(
paste(
"exdqlm_mcmc_pre_latent (iter=%d mode=%s reason=%s) invalid state before Ut/st update:",
"sigma=%s gamma=%s reg1_finite=%s theta_finite=%s st_finite=%s ut_finite=%s",
"max_abs_reg1=%s max_abs_theta=%s"
),
i,
latent_ctrl$mode,
latent_state_sched$reason,
format(cursam.sigma, digits = 6),
format(cursam.gamma, digits = 6),
all(is.finite(reg1)),
all(is.finite(cursam.theta)),
all(is.finite(cursam.st)),
all(is.finite(cursam.Ut)),
max_abs_finite(reg1),
max_abs_finite(cursam.theta)
))
}
if (isTRUE(latent_state_sched$active)) {
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"exdqlm_mcmc_latent_warmup_state (iter=%d mode=%s) invalid warmup Ut state: nonfinite=%d nonpositive=%d",
i, latent_ctrl$mode, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut)
))
}
if (identical(latent_ctrl$mode, "u_only")) {
cursam.st <- ex_samp_sts(reg1, cursam.gamma, cursam.sigma, cursam.Ut, a_tau, b_tau, c_tau)
if (count_nonfinite(cursam.st) > 0L || count_nonpositive(cursam.st) > 0L) {
stop(sprintf(
"exdqlm_mcmc_post_sts (iter=%d mode=%s) invalid st draws after Ut warmup: nonfinite=%d nonpositive=%d",
i, latent_ctrl$mode, count_nonfinite(cursam.st), count_nonpositive(cursam.st)
))
}
} else if (count_nonfinite(cursam.st) > 0L || count_nonpositive(cursam.st) > 0L) {
stop(sprintf(
"exdqlm_mcmc_latent_warmup_state (iter=%d mode=%s) invalid warmup st state: nonfinite=%d nonpositive=%d",
i, latent_ctrl$mode, count_nonfinite(cursam.st), count_nonpositive(cursam.st)
))
}
} else {
cursam.Ut <- ex_samp_uts(reg1, cursam.gamma, cursam.sigma, cursam.st, a_tau, b_tau, c_tau)
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"exdqlm_mcmc_post_uts (iter=%d mode=%s) invalid Ut draws: nonfinite=%d nonpositive=%d sigma=%s gamma=%s max_abs_reg1=%s",
i, latent_ctrl$mode, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut),
format(cursam.sigma, digits = 6), format(cursam.gamma, digits = 6), max_abs_finite(reg1)
))
}
cursam.st <- ex_samp_sts(reg1, cursam.gamma, cursam.sigma, cursam.Ut, a_tau, b_tau, c_tau)
if (count_nonfinite(cursam.st) > 0L || count_nonpositive(cursam.st) > 0L) {
stop(sprintf(
"exdqlm_mcmc_post_sts (iter=%d mode=%s) invalid st draws: nonfinite=%d nonpositive=%d sigma=%s gamma=%s max_abs_reg1=%s",
i, latent_ctrl$mode, count_nonfinite(cursam.st), count_nonpositive(cursam.st),
format(cursam.sigma, digits = 6), format(cursam.gamma, digits = 6), max_abs_finite(reg1)
))
}
}
latent_frozen_trace[i] <- isTRUE(latent_state_sched$active)
latent_update_reason_trace[i] <- latent_state_sched$reason
if (!isTRUE(latent_state_sched$active)) {
latent_update_performed_trace[i] <- TRUE
latent_update_count <- latent_update_count + 1L
latent_update_count_trace[i] <- latent_update_count
if (is.na(latent_first_active_iter)) latent_first_active_iter <- as.integer(i)
if (i <= n.burn) {
latent_updates_burn <- latent_updates_burn + 1L
} else {
latent_updates_keep <- latent_updates_keep + 1L
}
if (latent_ctrl$freeze_burnin_iters > 0L) {
latent_postwarmup_update_count <- latent_postwarmup_update_count + 1L
}
latent_forced_postwarmup_trace[i] <- isTRUE(latent_state_sched$force_now)
} else {
latent_update_count_trace[i] <- latent_update_count
latent_forced_postwarmup_trace[i] <- FALSE
}
# sample sigma and gamma
sigmagam_sched <- warmup_state(i, sigmagam_ctrl, sigmagam_postwarmup_update_count)
sigmagam_warmup_active <- isTRUE(sigmagam_sched$active)
sigmagam_force_now <- isTRUE(sigmagam_sched$force_now)
sigmagam_update_reason <- sigmagam_sched$reason
if (isTRUE(sigmagam_warmup_active)) {
lsiglgam.out <- list(
log.sigma = cursam.log.sigma,
logit.gamma = cursam.logit.gamma,
accept = NA,
slice_evals = NA_integer_
)
} else if (identical(mh.proposal, "slice")) {
if (!fix.sigma) {
cursam.sigma <- samp_sigma_exact(reg1, cursam.sigma, cursam.gamma, cursam.st, cursam.Ut)
cursam.log.sigma <- log(cursam.sigma)
}
gamma_log_density <- make_logpost_gamma(reg1, cursam.sigma, cursam.st, cursam.Ut)
slice_evals <- NA_integer_
if (!fix.gamma) {
current_lp <- gamma_log_density(cursam.gamma)
if (!is.finite(current_lp)) {
cursam.gamma <- min(max(cursam.gamma, L + 1e-8), U - 1e-8)
current_lp <- gamma_log_density(cursam.gamma)
}
if (!is.finite(current_lp)) {
cursam.gamma <- min(max(0, L + 1e-8), U - 1e-8)
}
slice_out <- .exdqlm_uni_slice_bounded(
x0 = cursam.gamma,
log_density = gamma_log_density,
w = slice.width,
m = slice.max.steps,
lower = L + 1e-10,
upper = U - 1e-10
)
cursam.gamma <- as.numeric(slice_out$value)[1]
slice_evals <- as.integer(slice_out$evals)
}
cursam.logit.gamma <- logit(cursam.gamma)
lsiglgam.out <- list(
log.sigma = cursam.log.sigma,
logit.gamma = cursam.logit.gamma,
accept = NA,
slice_evals = slice_evals
)
} else {
lsiglgam.out <- ex_samp_lsiglgam(reg1, cursam.log.sigma, cursam.logit.gamma, cursam.st, cursam.Ut, chol_Sig.mh)
cursam.gamma <- inv.logit(lsiglgam.out$logit.gamma)
cursam.logit.gamma <- lsiglgam.out$logit.gamma
cursam.sigma <- exp(lsiglgam.out$log.sigma)
cursam.log.sigma <- lsiglgam.out$log.sigma
n.accept <- n.accept + lsiglgam.out$accept
}
sigmagam_frozen_trace[i] <- isTRUE(sigmagam_warmup_active)
sigmagam_update_reason_trace[i] <- sigmagam_update_reason
if (!isTRUE(sigmagam_warmup_active)) {
sigmagam_update_performed_trace[i] <- TRUE
sigmagam_update_count <- sigmagam_update_count + 1L
sigmagam_update_count_trace[i] <- sigmagam_update_count
if (is.na(sigmagam_first_active_iter)) sigmagam_first_active_iter <- as.integer(i)
if (i <= n.burn) {
sigmagam_updates_burn <- sigmagam_updates_burn + 1L
} else {
sigmagam_updates_keep <- sigmagam_updates_keep + 1L
}
if (sigmagam_ctrl$freeze_burnin_iters > 0L) {
sigmagam_postwarmup_update_count <- sigmagam_postwarmup_update_count + 1L
}
sigmagam_forced_postwarmup_trace[i] <- isTRUE(sigmagam_force_now)
} else {
sigmagam_update_count_trace[i] <- sigmagam_update_count
sigmagam_forced_postwarmup_trace[i] <- FALSE
}
if (trace.diagnostics && (i %% trace.every == 0L)) {
s_stats <- .exdqlm_trace_summary(cursam.st)
u_stats <- .exdqlm_trace_summary(cursam.Ut)
trace_idx <- trace_idx + 1L
trace_rows[[trace_idx]] <- data.frame(
iter = i,
phase = if (i <= n.burn) "burn" else "keep",
sigma = cursam.sigma,
gamma = cursam.gamma,
accepted = if (identical(mh.proposal, "slice")) NA else isTRUE(lsiglgam.out$accept),
mh_scale = if (identical(mh.proposal, "slice")) NA_real_ else mh.scale,
latent_frozen = isTRUE(latent_state_sched$active),
latent_mode = latent_ctrl$mode,
latent_update_reason = latent_state_sched$reason,
latent_forced_postwarmup = isTRUE(latent_forced_postwarmup_trace[i]),
latent_update_performed = isTRUE(latent_update_performed_trace[i]),
latent_update_count = as.integer(latent_update_count_trace[i]),
theta_frozen = isTRUE(theta_state_sched$active),
theta_update_reason = theta_state_sched$reason,
theta_forced_postwarmup = isTRUE(theta_forced_postwarmup_trace[i]),
theta_update_performed = isTRUE(theta_update_performed_trace[i]),
theta_update_count = as.integer(theta_update_count_trace[i]),
sigmagam_frozen = isTRUE(sigmagam_warmup_active),
sigmagam_update_reason = sigmagam_update_reason,
sigmagam_forced_postwarmup = isTRUE(sigmagam_forced_postwarmup_trace[i]),
sigmagam_update_performed = isTRUE(sigmagam_update_performed_trace[i]),
sigmagam_update_count = as.integer(sigmagam_update_count_trace[i]),
slice_evals = if (!is.null(lsiglgam.out$slice_evals)) lsiglgam.out$slice_evals else NA_integer_,
s_mean = s_stats[["mean"]],
s_sd = s_stats[["sd"]],
s_q05 = s_stats[["q05"]],
s_q50 = s_stats[["median"]],
s_q95 = s_stats[["q95"]],
s_min = s_stats[["min"]],
s_max = s_stats[["max"]],
u_mean = u_stats[["mean"]],
u_sd = u_stats[["sd"]],
u_q05 = u_stats[["q05"]],
u_q50 = u_stats[["median"]],
u_q95 = u_stats[["q95"]],
u_min = u_stats[["min"]],
u_max = u_stats[["max"]],
stringsAsFactors = FALSE
)
}
# save samples after burn
if(i <= n.burn){
if (!identical(mh.proposal, "slice") && !isTRUE(sigmagam_warmup_active)) {
n.trial.burn <- n.trial.burn + 1L
n.accept.burn <- n.accept.burn + as.integer(isTRUE(lsiglgam.out$accept))
window.accept <- window.accept + as.integer(isTRUE(lsiglgam.out$accept))
window.total <- window.total + 1L
}
init.log.sigma[i] = cursam.log.sigma
init.logit.gamma[i] = cursam.logit.gamma
laplace_refreshed <- FALSE
if (identical(mh.proposal, "laplace_rw") && !fix.gamma && !fix.sigma &&
(!isTRUE(sigmagam_ctrl$delay_laplace_refresh_until_after_warmup) || !isTRUE(sigmagam_warmup_active)) &&
i >= mh.laplace.refresh.start && i < n.burn &&
(i %% mh.laplace.refresh.interval == 0)) {
laplace_refresh_attempts <- laplace_refresh_attempts + 1L
cov.lap.step <- laplace_cov_init(reg1, cursam.log.sigma, cursam.logit.gamma, cursam.st, cursam.Ut)
if (!is.null(cov.lap.step) && all(is.finite(cov.lap.step))) {
cov.lap.step <- prep_Sig_mh(cov.lap.step)
Sig.mh <- prep_Sig_mh((1 - mh.laplace.refresh.weight) * Sig.mh + mh.laplace.refresh.weight * cov.lap.step)
laplace_refreshed <- TRUE
laplace_refresh_success <- laplace_refresh_success + 1L
}
}
if (!identical(mh.proposal, "slice") &&
(!isTRUE(sigmagam_ctrl$delay_adapt_until_after_warmup) || !isTRUE(sigmagam_warmup_active)) &&
mh.adapt && i >= mh.min_burn_adapt && i < n.burn && (i %% mh.adapt.interval == 0)) {
acc.win <- window.accept / pmax(window.total, 1L)
if (acc.win < mh.target.accept[1]) {
mh.scale <- mh.scale * (1 - mh.max_scale.step)
} else if (acc.win > mh.target.accept[2]) {
mh.scale <- mh.scale * (1 + mh.max_scale.step)
}
mh.scale <- min(max(mh.scale, mh.scale.bounds[1]), mh.scale.bounds[2])
Sig.scaled <- prep_Sig_mh(Sig.mh * (mh.scale^2))
chol_Sig.mh <- build_chol(Sig.scaled)
adapt.history <- rbind(
adapt.history,
data.frame(
iter = i,
window_accept = acc.win,
mh_scale = mh.scale,
sig11 = Sig.scaled[1, 1],
sig22 = Sig.scaled[2, 2],
laplace_refreshed = isTRUE(laplace_refreshed),
stringsAsFactors = FALSE
)
)
window.accept <- 0L
window.total <- 0L
}
if(!identical(mh.proposal, "slice") && i==n.burn && joint.sample && !isTRUE(sigmagam_warmup_active)){
Sig.mh = stats::cov(cbind(init.log.sigma[1:n.burn],init.logit.gamma[1:n.burn]))
Sig.mh <- prep_Sig_mh(Sig.mh)
chol_Sig.mh <- build_chol(Sig.mh * (mh.scale^2))
}
}else{
if (!identical(mh.proposal, "slice") && !isTRUE(sigmagam_warmup_active)) {
n.trial.keep <- n.trial.keep + 1L
n.accept.keep <- n.accept.keep + as.integer(isTRUE(lsiglgam.out$accept))
}
save.sigma[(i-n.burn)] = cursam.sigma
save.gamma[(i-n.burn)] = cursam.gamma
save.theta[,,(i-n.burn)] = cursam.theta
save.Ut[,(i-n.burn)] = cursam.Ut
save.st[,(i-n.burn)] = cursam.st
save.post.pred[,(i-n.burn)] = theta.out$post.pred
}
}
run.time = tictoc::toc(quiet = TRUE)
.exdqlm_progress(
"MCMC done",
model = "exDQLM",
status = "complete",
iter = I,
runtime_sec = run.time$toc - run.time$tic,
accept = if (identical(mh.proposal, "slice")) NULL else n.accept / I,
.verbose = verbose
)
safe_progress_callback(list(
event = "complete",
iter = as.integer(I),
total_iter = as.integer(I),
phase = "done",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = cursam.gamma,
kernel = mh.proposal,
accept = if (identical(mh.proposal, "slice")) NA_real_ else n.accept / I,
runtime_sec = as.numeric(run.time$toc - run.time$tic)
))
# exdqlm MAP standard forecast errors
map.gam = mean(save.gamma)
map.sig = mean(save.sigma)
map.st = rowMeans(save.st)
map.Ut = rowMeans(save.Ut)
tau = p.fn(p0,map.gam)
a_tau = (1-2*tau)/(tau*(1-tau))
b_tau = (2)/(tau*(1-tau))
c_tau = (as.numeric(map.gam>0)-tau)^(-1)
theta.out <- smoothed_theta(map.sig*c_tau*abs(map.gam)*map.st+map.Ut*a_tau,b_tau*map.Ut*map.sig)
map.standard.forecast.errors = theta.out$standard.forecast.errors
Sig.mh.final <- prep_Sig_mh(Sig.mh * (mh.scale^2))
ess_sigma <- tryCatch(as.numeric(coda::effectiveSize(coda::as.mcmc(save.sigma))), error = function(e) NA_real_)
ess_gamma <- tryCatch(as.numeric(coda::effectiveSize(coda::as.mcmc(save.gamma))), error = function(e) NA_real_)
chain_health_sigma <- .exdqlm_chain_health_metrics(save.sigma, n_keep = n.mcmc)
chain_health_gamma <- .exdqlm_chain_health_metrics(save.gamma, n_keep = n.mcmc)
accept_total <- if (identical(mh.proposal, "slice")) NA_real_ else n.accept / I
accept_burn <- if (identical(mh.proposal, "slice")) NA_real_ else if (n.trial.burn > 0) n.accept.burn / n.trial.burn else NA_real_
accept_keep <- if (identical(mh.proposal, "slice")) NA_real_ else if (n.trial.keep > 0) n.accept.keep / n.trial.keep else NA_real_
kernel_exact <- mh.proposal %in% c("rw", "laplace_rw", "slice")
mh.diag <- list(
proposal = mh.proposal,
adapt = if (identical(mh.proposal, "slice")) FALSE else mh.adapt,
joint_sample = isTRUE(joint.sample),
adapt_interval = if (identical(mh.proposal, "slice")) NA_integer_ else mh.adapt.interval,
target_accept = if (identical(mh.proposal, "slice")) c(NA_real_, NA_real_) else mh.target.accept,
scale_bounds = if (identical(mh.proposal, "slice")) c(NA_real_, NA_real_) else mh.scale.bounds,
scale_final = if (identical(mh.proposal, "slice")) NA_real_ else mh.scale,
joint_sigma_gamma = mh.proposal %in% c("rw", "laplace_rw"),
transformed_state = if (mh.proposal %in% c("rw", "laplace_rw")) c("log_sigma", "logit_gamma") else c("gamma"),
# Backward-compatible aliases used by some diagnostics scripts.
final_scale = if (identical(mh.proposal, "slice")) NA_real_ else mh.scale,
slice_width = if (identical(mh.proposal, "slice")) slice.width else NA_real_,
slice_max_steps = if (identical(mh.proposal, "slice")) slice.max.steps else NA_real_,
laplace_refresh = list(
enabled = identical(mh.proposal, "laplace_rw"),
interval = if (identical(mh.proposal, "laplace_rw")) as.integer(mh.laplace.refresh.interval) else NA_integer_,
start = if (identical(mh.proposal, "laplace_rw")) as.integer(mh.laplace.refresh.start) else NA_integer_,
weight = if (identical(mh.proposal, "laplace_rw")) as.numeric(mh.laplace.refresh.weight) else NA_real_,
attempts = if (identical(mh.proposal, "laplace_rw")) as.integer(laplace_refresh_attempts) else NA_integer_,
success = if (identical(mh.proposal, "laplace_rw")) as.integer(laplace_refresh_success) else NA_integer_
),
kernel_exact = kernel_exact,
signoff_ready = kernel_exact,
approximation_note = NA_character_,
accept = list(
total = accept_total,
burn = accept_burn,
kept = accept_keep,
n_accept = if (identical(mh.proposal, "slice")) NA_integer_ else n.accept,
n_total = if (identical(mh.proposal, "slice")) NA_integer_ else I
),
Sig.mh.initial = if (identical(mh.proposal, "slice")) matrix(NA_real_, 2, 2) else Sig.mh.initial,
Sig.mh.final = if (identical(mh.proposal, "slice")) matrix(NA_real_, 2, 2) else Sig.mh.final,
adaptation = if (identical(mh.proposal, "slice")) data.frame() else adapt.history,
adapt_trace = if (identical(mh.proposal, "slice")) data.frame() else adapt.history,
trace_enabled = trace.diagnostics,
trace_every = if (trace.diagnostics) trace.every else NA_integer_,
verbose_every = as.integer(verbose.every),
callback_every = as.integer(callback.every),
sigmagam = list(
freeze_burnin_iters = as.integer(sigmagam_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(sigmagam_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(sigmagam_ctrl$force_after_warmup),
delay_adapt_until_after_warmup = isTRUE(sigmagam_ctrl$delay_adapt_until_after_warmup),
delay_laplace_refresh_until_after_warmup = isTRUE(sigmagam_ctrl$delay_laplace_refresh_until_after_warmup),
first_active_iter = if (is.na(sigmagam_first_active_iter)) NA_integer_ else as.integer(sigmagam_first_active_iter),
updates_burn = as.integer(sigmagam_updates_burn),
updates_keep = as.integer(sigmagam_updates_keep),
update_count = as.integer(sigmagam_update_count),
postwarmup_update_count = as.integer(sigmagam_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(sigmagam_frozen_trace[seq_len(n.burn)]) else NA_real_
),
latent_state = list(
mode = latent_ctrl$mode,
freeze_burnin_iters = as.integer(latent_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(latent_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(latent_ctrl$force_after_warmup),
first_active_iter = if (is.na(latent_first_active_iter)) NA_integer_ else as.integer(latent_first_active_iter),
updates_burn = as.integer(latent_updates_burn),
updates_keep = as.integer(latent_updates_keep),
update_count = as.integer(latent_update_count),
postwarmup_update_count = as.integer(latent_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(latent_frozen_trace[seq_len(n.burn)]) else NA_real_
),
theta_state = list(
freeze_burnin_iters = as.integer(theta_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(theta_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(theta_ctrl$force_after_warmup),
first_active_iter = if (is.na(theta_first_active_iter)) NA_integer_ else as.integer(theta_first_active_iter),
updates_burn = as.integer(theta_updates_burn),
updates_keep = as.integer(theta_updates_keep),
update_count = as.integer(theta_update_count),
postwarmup_update_count = as.integer(theta_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(theta_frozen_trace[seq_len(n.burn)]) else NA_real_
),
trace = if (trace.diagnostics && trace_idx > 0L) {
do.call(rbind, trace_rows[seq_len(trace_idx)])
} else {
data.frame()
}
)
# exdqlm results
retlist = list(y=y,run.time=(run.time$toc-run.time$tic),model=model,p0=p0,df=df,dim.df=dim.df,dqlm.ind=dqlm.ind,
samp.theta = coda::as.mcmc(save.theta), theta.out = theta.out,
samp.post.pred = save.post.pred, map.standard.forecast.errors = map.standard.forecast.errors,
samp.sigma = coda::as.mcmc(save.sigma), samp.gamma = coda::as.mcmc(save.gamma),
init.log.sigma = coda::as.mcmc(init.log.sigma), init.logit.gamma = coda::as.mcmc(init.logit.gamma),
samp.vts = coda::as.mcmc(save.Ut), samp.sts = coda::as.mcmc(save.st),
accept.rate = accept_total,
accept.rate.burn = accept_burn,
accept.rate.keep = accept_keep,
Sig.mh = if (identical(mh.proposal, "slice")) matrix(NA_real_, 2, 2) else Sig.mh.final,
init.from.vb = init.from.vb,
vb.init.method = if (init.from.vb) vb.ctrl$method else NA_character_,
mh.diagnostics = mh.diag,
diagnostics = list(
mh = mh.diag,
progress = list(
verbose_every = as.integer(verbose.every),
callback_every = as.integer(callback.every)
),
ess = list(sigma = ess_sigma, gamma = ess_gamma),
chain_health = list(
sigma = chain_health_sigma,
gamma = chain_health_gamma
),
latent_state = mh.diag$latent_state,
theta_state = mh.diag$theta_state,
sigmagam = mh.diag$sigmagam,
s_block = list(
trace = mh.diag$trace,
final = if (is.data.frame(mh.diag$trace) && nrow(mh.diag$trace)) {
as.list(mh.diag$trace[nrow(mh.diag$trace), , drop = FALSE])
} else {
list()
}
),
rhat_ready = list(
sigma = as.numeric(save.sigma),
gamma = as.numeric(save.gamma)
),
sigmagam_trace = list(
frozen = sigmagam_frozen_trace,
update_reason = sigmagam_update_reason_trace,
forced_postwarmup = sigmagam_forced_postwarmup_trace,
update_performed = sigmagam_update_performed_trace,
update_count = sigmagam_update_count_trace
),
latent_state_trace = list(
frozen = latent_frozen_trace,
update_reason = latent_update_reason_trace,
forced_postwarmup = latent_forced_postwarmup_trace,
update_performed = latent_update_performed_trace,
update_count = latent_update_count_trace
),
theta_state_trace = list(
frozen = theta_frozen_trace,
update_reason = theta_update_reason_trace,
forced_postwarmup = theta_forced_postwarmup_trace,
update_performed = theta_update_performed_trace,
update_count = theta_update_count_trace
)
),
n.burn=n.burn,n.mcmc=n.mcmc)
}else{
######## DQLM
# fixed AL parameters
a_tau = (1-2*p0)/(p0*(1-p0))
b_tau = (2)/(p0*(1-p0))
# dqlm function sample theta ffbs
samp_theta<-function(ex.f,ex.q,sigma){
# initialize ffbs
m <- sam.theta <- matrix(NA,p,TT)
C <- array(NA,c(p,p,TT))
standard.forecast.errors <- post.pred <- rep(NA,TT)
## forward filter
# first iteration
a = as.vector(GG[,,1]%*%m0)
P = .exdqlm_regularize_cov(GG[,,1]%*%C0%*%t(GG[,,1]), context = "mcmc_al_sample_P_t1")
R = .exdqlm_regularize_cov(P + df.mat*P, context = "mcmc_al_sample_R_t1")
f = t(FF[,1])%*%a + ex.f[1]
q = .exdqlm_regularize_var(t(FF[,1])%*%R%*%FF[,1] + ex.q[1], context = "mcmc_al_sample_q_t1")
m[,1] = a + t(R)%*%FF[,1]%*%(y[1]-f)/q[1]
C[,,1] = .exdqlm_regularize_cov(
R - t(R)%*%FF[,1]%*%t(FF[,1])%*%R/q[1],
context = "mcmc_al_sample_C_t1"
)
standard.forecast.errors[1] = (y[1]-f)/sqrt(q)
# t = 2:TT
for(t in 2:TT){
a = as.vector(GG[,,t]%*%m[,(t-1)])
P = .exdqlm_regularize_cov(GG[,,t]%*%C[,,(t-1)]%*%t(GG[,,t]), context = sprintf("mcmc_al_sample_P_t%d", t))
R = .exdqlm_regularize_cov(P + df.mat*P, context = sprintf("mcmc_al_sample_R_t%d", t))
f = t(FF[,t])%*%a + ex.f[t]
fB = t(FF[,t])%*%R
q = .exdqlm_regularize_var(fB%*%FF[,t] + ex.q[t], context = sprintf("mcmc_al_sample_q_t%d", t))
m[,t] = a + t(fB)%*%(y[t]-f)/q[1]
C[,,t] = .exdqlm_regularize_cov(
R - t(fB)%*%fB/q[1],
context = sprintf("mcmc_al_sample_C_t%d", t)
)
standard.forecast.errors[t] = (y[t]-f)/sqrt(q)
}
## backwards sample
sC_TT = .exdqlm_regularize_cov(C[,,TT], context = "mcmc_al_sample_sC_TT")
svd.sC = svd(sC_TT)
sam.theta[,TT] = m[,TT] + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta <- numeric(TT)
reg_theta[TT] <- drop(crossprod(FF[,TT], sam.theta[,TT]))
post.pred[TT] = rexal(1,p0,reg_theta[TT],sigma,0)
for(t in (TT-1):1){
P = .exdqlm_regularize_cov(GG[,,(t+1)]%*%C[,,(t)]%*%t(GG[,,(t+1)]), context = sprintf("mcmc_al_back_P_t%d", t + 1L))
R.info = .exdqlm_cov_inverse(P + df.mat*P, context = sprintf("mcmc_al_back_R_t%d", t + 1L))
sB = C[,,t]%*%t(GG[,,(t+1)])%*%R.info$inverse
sm = m[,t] + sB%*%(sam.theta[,(t+1)]-as.vector(GG[,,(t+1)]%*%m[,(t)]))
sC = .exdqlm_regularize_cov(
C[,,t] - sB%*%GG[,,(t+1)]%*%C[,,t],
context = sprintf("mcmc_al_back_sC_t%d", t)
)
svd.sC = svd(sC)
sam.theta[,t] = sm + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta[t] <- drop(crossprod(FF[,t], sam.theta[,t]))
post.pred[t] = rexal(1,p0,reg_theta[t],sigma,0)
}
return(list(standard.forecast.errors=standard.forecast.errors,post.pred=post.pred,sam.theta=sam.theta,fm=m,fC=C))
}
if (use_cpp_mcmc) {
samp_theta <- function(ex.f, ex.q, sigma) {
out <- cpp_ffbs_sample(ex.f, ex.q)
sam.theta <- out$sam.theta
reg_theta <- state_signal(FF, sam.theta)
post.pred <- vapply(seq_len(TT), function(t) {
rexal(1, p0, reg_theta[t], sigma, 0)
}, numeric(1))
list(
standard.forecast.errors = out$standard.forecast.errors,
post.pred = post.pred,
sam.theta = sam.theta,
fm = out$fm,
fC = out$fC
)
}
}
# dqlm function sample uts
samp_uts<-function(reg1,sigma){
chi <- as.numeric(((y-reg1)^2)/(b_tau*sigma))
psi <- (a_tau^2)/(b_tau*sigma) + (2/sigma)
sample_gig_cpp_required(chi = chi, psi = psi, lambda = 0.5, context = "dqlm_mcmc_uts")
}
# dqlm function sample sigma
samp_sigma<-function(reg1,uts){
1/stats::rgamma(1, shape = PriorSigma$a_sig + 1.5*TT,
rate = PriorSigma$b_sig + 0.5*sum( ((as.vector(y) - reg1 - a_tau*uts)^2)/(b_tau*uts) ) + sum(uts) )
}
callback.every <- if (trace.diagnostics) {
trace.every
} else {
100L
}
callback.every <- max(1L, as.integer(callback.every)[1])
trace_rows <- if (trace.diagnostics) vector("list", ceiling(I / trace.every)) else NULL
trace_idx <- 0L
latent_frozen_trace <- rep(FALSE, I)
latent_forced_postwarmup_trace <- rep(FALSE, I)
latent_update_performed_trace <- rep(FALSE, I)
latent_update_reason_trace <- rep(NA_character_, I)
latent_update_count_trace <- integer(I)
latent_first_active_iter <- NA_integer_
latent_update_count <- 0L
latent_postwarmup_update_count <- 0L
latent_updates_burn <- 0L
latent_updates_keep <- 0L
theta_frozen_trace <- rep(FALSE, I)
theta_forced_postwarmup_trace <- rep(FALSE, I)
theta_update_performed_trace <- rep(FALSE, I)
theta_update_reason_trace <- rep(NA_character_, I)
theta_update_count_trace <- integer(I)
theta_first_active_iter <- NA_integer_
theta_update_count <- 0L
theta_postwarmup_update_count <- 0L
theta_updates_burn <- 0L
theta_updates_keep <- 0L
dqlm_sigma_frozen_trace <- rep(FALSE, I)
dqlm_sigma_forced_postwarmup_trace <- rep(FALSE, I)
dqlm_sigma_update_performed_trace <- rep(FALSE, I)
dqlm_sigma_update_reason_trace <- rep(NA_character_, I)
dqlm_sigma_update_count_trace <- integer(I)
dqlm_sigma_first_active_iter <- NA_integer_
dqlm_sigma_update_count <- 0L
dqlm_sigma_postwarmup_update_count <- 0L
dqlm_sigma_updates_burn <- 0L
dqlm_sigma_updates_keep <- 0L
# Sample from dqlm posterior
tictoc::tic()
.exdqlm_progress(
"MCMC start",
burn = n.burn,
keep = n.mcmc,
kernel = "conjugate",
warm_start = if (isTRUE(init.from.vb)) vb.ctrl$method else "none",
.verbose = verbose
)
safe_progress_callback(list(
event = "start",
iter = 0L,
total_iter = as.integer(I),
phase = "burn",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = NA_real_,
kernel = "conjugate",
accept = NA_real_
))
for (i in 1:I){
current_iter <- as.integer(i)
# counter
if (i %% verbose.every == 0L) {
kept_now <- if (i <= n.burn) 0L else (i - n.burn)
.exdqlm_progress(
"MCMC progress",
model = "DQLM",
phase = if (i <= n.burn) "burn" else "keep",
iter = sprintf("%d/%d", i, I),
kept = sprintf("%d/%d", kept_now, n.mcmc),
.verbose = verbose
)
}
if (i %% callback.every == 0L) {
safe_progress_callback(list(
event = "progress",
iter = as.integer(i),
total_iter = as.integer(I),
phase = if (i <= n.burn) "burn" else "keep",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = NA_real_,
kernel = "conjugate",
accept = NA_real_
))
}
# sample theta
ex.f = cursam.Ut*a_tau
ex.q = b_tau*cursam.Ut*cursam.sigma
theta_state_sched <- warmup_state(i, theta_ctrl, theta_postwarmup_update_count)
if (isTRUE(theta_state_sched$active)) {
if (any(!is.finite(cursam.theta))) {
stop(sprintf(
"dqlm_mcmc_theta_warmup_state (iter=%d) invalid warmup theta state: nonfinite=%d",
i, sum(!is.finite(cursam.theta))
))
}
reg_theta_hold <- state_signal(FF, cursam.theta)
theta.out <- list(
standard.forecast.errors = rep(NA_real_, TT),
post.pred = vapply(seq_len(TT), function(t) {
rexal(1, p0, reg_theta_hold[t], cursam.sigma, 0)
}, numeric(1)),
sam.theta = cursam.theta,
fm = matrix(NA_real_, p, TT),
fC = array(NA_real_, c(p, p, TT))
)
} else {
theta.out <- samp_theta(ex.f,ex.q,cursam.sigma)
cursam.theta = theta.out$sam.theta
}
theta_frozen_trace[i] <- isTRUE(theta_state_sched$active)
theta_update_reason_trace[i] <- theta_state_sched$reason
if (!isTRUE(theta_state_sched$active)) {
theta_update_performed_trace[i] <- TRUE
theta_update_count <- theta_update_count + 1L
theta_update_count_trace[i] <- theta_update_count
if (is.na(theta_first_active_iter)) theta_first_active_iter <- as.integer(i)
if (i <= n.burn) {
theta_updates_burn <- theta_updates_burn + 1L
} else {
theta_updates_keep <- theta_updates_keep + 1L
}
if (theta_ctrl$freeze_burnin_iters > 0L) {
theta_postwarmup_update_count <- theta_postwarmup_update_count + 1L
}
theta_forced_postwarmup_trace[i] <- isTRUE(theta_state_sched$force_now)
} else {
theta_update_count_trace[i] <- theta_update_count
theta_forced_postwarmup_trace[i] <- FALSE
}
# sample uts
reg1 = state_signal(FF, cursam.theta)
latent_state_sched <- warmup_state(i, latent_ctrl, latent_postwarmup_update_count)
if (!is.finite(cursam.sigma) || cursam.sigma <= 0 || any(!is.finite(reg1))) {
stop(sprintf(
"dqlm_mcmc_pre_uts (iter=%d reason=%s) invalid state before chi update: sigma=%s reg1_finite=%s max_abs_reg1=%s max_abs_theta=%s",
i,
latent_state_sched$reason,
format(cursam.sigma, digits = 6),
all(is.finite(reg1)),
if (all(is.finite(reg1))) format(max(abs(reg1), na.rm = TRUE), digits = 6) else "NA",
if (all(is.finite(cursam.theta))) format(max(abs(cursam.theta), na.rm = TRUE), digits = 6) else "NA"
))
}
if (isTRUE(latent_state_sched$active)) {
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"dqlm_mcmc_latent_warmup_state (iter=%d) invalid warmup Ut state: nonfinite=%d nonpositive=%d",
i, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut)
))
}
} else {
cursam.Ut <- samp_uts(reg1,cursam.sigma)
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"dqlm_mcmc_post_uts (iter=%d) invalid Ut draws: nonfinite=%d nonpositive=%d sigma=%s max_abs_reg1=%s",
i, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut),
format(cursam.sigma, digits = 6), max_abs_finite(reg1)
))
}
}
latent_frozen_trace[i] <- isTRUE(latent_state_sched$active)
latent_update_reason_trace[i] <- latent_state_sched$reason
if (!isTRUE(latent_state_sched$active)) {
latent_update_performed_trace[i] <- TRUE
latent_update_count <- latent_update_count + 1L
latent_update_count_trace[i] <- latent_update_count
if (is.na(latent_first_active_iter)) latent_first_active_iter <- as.integer(i)
if (i <= n.burn) {
latent_updates_burn <- latent_updates_burn + 1L
} else {
latent_updates_keep <- latent_updates_keep + 1L
}
if (latent_ctrl$freeze_burnin_iters > 0L) {
latent_postwarmup_update_count <- latent_postwarmup_update_count + 1L
}
latent_forced_postwarmup_trace[i] <- isTRUE(latent_state_sched$force_now)
} else {
latent_update_count_trace[i] <- latent_update_count
latent_forced_postwarmup_trace[i] <- FALSE
}
# sample sigma
dqlm_sigma_sched <- warmup_state(i, dqlm_sigma_ctrl, dqlm_sigma_postwarmup_update_count)
if(!fix.sigma){
if (!isTRUE(dqlm_sigma_sched$active)) {
cursam.sigma <- samp_sigma(reg1,cursam.Ut)
}
if (!is.finite(cursam.sigma) || cursam.sigma <= 0) {
stop(sprintf(
"dqlm_mcmc_post_sigma (iter=%d reason=%s) invalid sigma draw: sigma=%s max_abs_reg1=%s max_abs_ut=%s",
i, dqlm_sigma_sched$reason, format(cursam.sigma, digits = 6), max_abs_finite(reg1), max_abs_finite(cursam.Ut)
))
}
}
dqlm_sigma_frozen_trace[i] <- isTRUE(dqlm_sigma_sched$active)
dqlm_sigma_update_reason_trace[i] <- dqlm_sigma_sched$reason
if (!isTRUE(dqlm_sigma_sched$active)) {
dqlm_sigma_update_performed_trace[i] <- TRUE
dqlm_sigma_update_count <- dqlm_sigma_update_count + 1L
dqlm_sigma_update_count_trace[i] <- dqlm_sigma_update_count
if (is.na(dqlm_sigma_first_active_iter)) dqlm_sigma_first_active_iter <- as.integer(i)
if (i <= n.burn) {
dqlm_sigma_updates_burn <- dqlm_sigma_updates_burn + 1L
} else {
dqlm_sigma_updates_keep <- dqlm_sigma_updates_keep + 1L
}
if (dqlm_sigma_ctrl$freeze_burnin_iters > 0L) {
dqlm_sigma_postwarmup_update_count <- dqlm_sigma_postwarmup_update_count + 1L
}
dqlm_sigma_forced_postwarmup_trace[i] <- isTRUE(dqlm_sigma_sched$force_now)
} else {
dqlm_sigma_update_count_trace[i] <- dqlm_sigma_update_count
dqlm_sigma_forced_postwarmup_trace[i] <- FALSE
}
if (trace.diagnostics && (i %% trace.every == 0L)) {
u_stats <- .exdqlm_trace_summary(cursam.Ut)
trace_idx <- trace_idx + 1L
trace_rows[[trace_idx]] <- data.frame(
iter = i,
phase = if (i <= n.burn) "burn" else "keep",
sigma = cursam.sigma,
latent_frozen = isTRUE(latent_state_sched$active),
latent_mode = latent_ctrl$mode,
latent_update_reason = latent_state_sched$reason,
latent_forced_postwarmup = isTRUE(latent_forced_postwarmup_trace[i]),
latent_update_performed = isTRUE(latent_update_performed_trace[i]),
latent_update_count = as.integer(latent_update_count_trace[i]),
theta_frozen = isTRUE(theta_state_sched$active),
theta_update_reason = theta_state_sched$reason,
theta_forced_postwarmup = isTRUE(theta_forced_postwarmup_trace[i]),
theta_update_performed = isTRUE(theta_update_performed_trace[i]),
theta_update_count = as.integer(theta_update_count_trace[i]),
dqlm_sigma_frozen = isTRUE(dqlm_sigma_sched$active),
dqlm_sigma_update_reason = dqlm_sigma_sched$reason,
dqlm_sigma_forced_postwarmup = isTRUE(dqlm_sigma_forced_postwarmup_trace[i]),
dqlm_sigma_update_performed = isTRUE(dqlm_sigma_update_performed_trace[i]),
dqlm_sigma_update_count = as.integer(dqlm_sigma_update_count_trace[i]),
u_mean = u_stats[["mean"]],
u_sd = u_stats[["sd"]],
u_q05 = u_stats[["q05"]],
u_q50 = u_stats[["median"]],
u_q95 = u_stats[["q95"]],
u_min = u_stats[["min"]],
u_max = u_stats[["max"]],
stringsAsFactors = FALSE
)
}
# save samples after burn
if(i > n.burn){
save.sigma[(i-n.burn)] = cursam.sigma
save.theta[,,(i-n.burn)] = cursam.theta
save.Ut[,(i-n.burn)] = cursam.Ut
save.post.pred[,(i-n.burn)] = theta.out$post.pred
}
}
run.time = tictoc::toc(quiet = TRUE)
.exdqlm_progress(
"MCMC done",
model = "DQLM",
status = "complete",
iter = I,
runtime_sec = run.time$toc - run.time$tic,
.verbose = verbose
)
safe_progress_callback(list(
event = "complete",
iter = as.integer(I),
total_iter = as.integer(I),
phase = "done",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = NA_real_,
kernel = "conjugate",
accept = NA_real_,
runtime_sec = as.numeric(run.time$toc - run.time$tic)
))
# dqlm MAP standard forecast errors
map.sig = mean(save.sigma)
map.Ut = rowMeans(save.Ut)
theta.out <- smoothed_theta(map.Ut*a_tau,b_tau*map.Ut*map.sig)
map.standard.forecast.errors = theta.out$standard.forecast.errors
ess_sigma <- tryCatch(as.numeric(coda::effectiveSize(coda::as.mcmc(save.sigma))), error = function(e) NA_real_)
# dqlm results
retlist = list(y=y,run.time=(run.time$toc-run.time$tic),model=model,p0=p0,df=df,dim.df=dim.df,dqlm.ind=dqlm.ind,
samp.theta = coda::as.mcmc(save.theta), theta.out = theta.out,
samp.post.pred = save.post.pred, map.standard.forecast.errors = map.standard.forecast.errors,
samp.sigma = coda::as.mcmc(save.sigma),
samp.vts = coda::as.mcmc(save.Ut),
init.from.vb = init.from.vb,
vb.init.method = if (init.from.vb) vb.ctrl$method else NA_character_,
diagnostics = list(
progress = list(
verbose_every = as.integer(verbose.every),
callback_every = as.integer(callback.every)
),
latent_state = list(
mode = latent_ctrl$mode,
freeze_burnin_iters = as.integer(latent_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(latent_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(latent_ctrl$force_after_warmup),
first_active_iter = if (is.na(latent_first_active_iter)) NA_integer_ else as.integer(latent_first_active_iter),
updates_burn = as.integer(latent_updates_burn),
updates_keep = as.integer(latent_updates_keep),
update_count = as.integer(latent_update_count),
postwarmup_update_count = as.integer(latent_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(latent_frozen_trace[seq_len(n.burn)]) else NA_real_
),
theta_state = list(
freeze_burnin_iters = as.integer(theta_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(theta_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(theta_ctrl$force_after_warmup),
first_active_iter = if (is.na(theta_first_active_iter)) NA_integer_ else as.integer(theta_first_active_iter),
updates_burn = as.integer(theta_updates_burn),
updates_keep = as.integer(theta_updates_keep),
update_count = as.integer(theta_update_count),
postwarmup_update_count = as.integer(theta_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(theta_frozen_trace[seq_len(n.burn)]) else NA_real_
),
dqlm_sigma = list(
freeze_burnin_iters = as.integer(dqlm_sigma_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(dqlm_sigma_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(dqlm_sigma_ctrl$force_after_warmup),
first_active_iter = if (is.na(dqlm_sigma_first_active_iter)) NA_integer_ else as.integer(dqlm_sigma_first_active_iter),
updates_burn = as.integer(dqlm_sigma_updates_burn),
updates_keep = as.integer(dqlm_sigma_updates_keep),
update_count = as.integer(dqlm_sigma_update_count),
postwarmup_update_count = as.integer(dqlm_sigma_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(dqlm_sigma_frozen_trace[seq_len(n.burn)]) else NA_real_
),
latent_state_trace = list(
frozen = latent_frozen_trace,
update_reason = latent_update_reason_trace,
forced_postwarmup = latent_forced_postwarmup_trace,
update_performed = latent_update_performed_trace,
update_count = latent_update_count_trace
),
theta_state_trace = list(
frozen = theta_frozen_trace,
update_reason = theta_update_reason_trace,
forced_postwarmup = theta_forced_postwarmup_trace,
update_performed = theta_update_performed_trace,
update_count = theta_update_count_trace
),
dqlm_sigma_trace = list(
frozen = dqlm_sigma_frozen_trace,
update_reason = dqlm_sigma_update_reason_trace,
forced_postwarmup = dqlm_sigma_forced_postwarmup_trace,
update_performed = dqlm_sigma_update_performed_trace,
update_count = dqlm_sigma_update_count_trace
),
trace = if (trace.diagnostics && trace_idx > 0L) {
do.call(rbind, trace_rows[seq_len(trace_idx)])
} else {
data.frame()
},
ess = list(sigma = ess_sigma, gamma = NA_real_),
rhat_ready = list(
sigma = as.numeric(save.sigma),
gamma = numeric(0)
)
),
n.burn=n.burn,n.mcmc=n.mcmc)
}
retlist$backend <- list(mcmc = mcmc_backend, mode = cpp_mcmc_mode, gig = gig_backend)
# return results
class(retlist) <- .exdqlm_fit_class("exdqlmMCMC")
return(retlist)
}
# Internal helper for diagnostics-only multichain validation.
.exdqlm_mcmc_multichain_diag <- function(n.chains = 4L, seeds = NULL, mcmc_args = list()) {
n.chains <- suppressWarnings(as.integer(n.chains)[1])
if (!is.finite(n.chains) || n.chains < 2L) {
stop("n.chains must be >= 2 for multichain diagnostics.")
}
if (is.null(seeds)) {
seeds <- seq_len(n.chains) + 20260300L
}
seeds <- as.integer(seeds)
if (length(seeds) != n.chains) {
stop("Length of seeds must match n.chains.")
}
fits <- vector("list", n.chains)
for (i in seq_len(n.chains)) {
set.seed(seeds[i])
args_i <- utils::modifyList(mcmc_args, list(verbose = FALSE))
fits[[i]] <- do.call(exdqlmMCMC, args_i)
}
sigma_list <- coda::mcmc.list(lapply(fits, function(f) coda::as.mcmc(as.numeric(f$samp.sigma))))
sigma_rhat <- tryCatch(
as.numeric(coda::gelman.diag(sigma_list, autoburnin = FALSE)$psrf[1, "Point est."]),
error = function(e) NA_real_
)
sigma_ess <- tryCatch(as.numeric(coda::effectiveSize(sigma_list))[1], error = function(e) NA_real_)
has_gamma <- !isTRUE(fits[[1]]$dqlm.ind) && !is.null(fits[[1]]$samp.gamma)
if (has_gamma) {
gamma_list <- coda::mcmc.list(lapply(fits, function(f) coda::as.mcmc(as.numeric(f$samp.gamma))))
gamma_rhat <- tryCatch(
as.numeric(coda::gelman.diag(gamma_list, autoburnin = FALSE)$psrf[1, "Point est."]),
error = function(e) NA_real_
)
gamma_ess <- tryCatch(as.numeric(coda::effectiveSize(gamma_list))[1], error = function(e) NA_real_)
} else {
gamma_list <- NULL
gamma_rhat <- NA_real_
gamma_ess <- NA_real_
}
list(
fits = fits,
seeds = seeds,
diagnostics = list(
sigma = list(rhat = sigma_rhat, ess = sigma_ess, chains = sigma_list),
gamma = list(rhat = gamma_rhat, ess = gamma_ess, chains = gamma_list)
)
)
}
Try the exdqlm package in your browser
Any scripts or data that you put into this service are public.
exdqlm documentation built on July 10, 2026, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .exdqlm_mcmc_multichain_diag exdqlmMCMC .exdqlm_mcmc_warmup_state .exdqlm_theta_state_mcmc_controls .exdqlm_dqlm_sigma_mcmc_controls .exdqlm_latent_state_mcmc_controls
Documented in exdqlmMCMC
#' exDQLM - MCMC algorithm
#'
#' The function applies a Markov chain Monte Carlo (MCMC) algorithm to sample the posterior of an exDQLM.
#'
#' @param y A univariate time-series.
#' @param p0 The quantile of interest, a value between 0 and 1.
#' @param model List of the state-space model including \code{GG}, \code{FF},
#' prior parameters \code{m0} and \code{C0}.
#' @param df Discount factors for each block.
#' @param dim.df Dimension of each block of discount factors.
#' @param fix.gamma Logical value indicating whether to fix gamma at
#' \code{gam.init}. Default is \code{FALSE}.
#' @param gam.init Initial value for gamma (skewness parameter), or value at
#' which gamma will be fixed if \code{fix.gamma = TRUE}.
#' @param fix.sigma Logical value indicating whether to fix sigma at `sig.init`.
#' Default is `FALSE`.
#' @param sig.init Initial value for sigma (scale parameter), or value at which
#' sigma will be fixed if \code{fix.sigma = TRUE}.
#' @param dqlm.ind Logical value indicating whether to fix gamma at \code{0},
#' reducing the exDQLM to the AL/DQLM special case. Default is \code{FALSE}.
#' @param Sig.mh Covariance matrix used in the random walk MH step to jointly sample sigma and gamma.
#' @param joint.sample Logical value indicating whether or not to recompute `Sig.mh` based off the initial burn-in samples of gamma and sigma. Default is `FALSE`.
#' @param n.burn Number of MCMC iterations to burn. Default is `n.burn = 2000`.
#' @param n.mcmc Number of MCMC iterations to sample. Default is `n.mcmc = 1500`.
#' @param init.from.isvb Logical value indicating whether to use the legacy ISVB
#' warm start when `init.from.vb = TRUE`. Default is `FALSE`, which favors
#' `LDVB` as the default VB warm start. This flag only chooses the
#' warm-start source; it does not change the subsequent MCMC proposal kernel.
#' @param init.from.vb Optional logical. If `TRUE`, run a VB pre-initialization
#' step (`LDVB` by default, or `ISVB` when `init.from.isvb = TRUE`) and
#' initialize MCMC from converged VB moments. Default is `TRUE`. If
#' explicitly set to `NULL`, it falls back to `init.from.isvb` behavior for
#' backward compatibility.
#' @param vb_init_controls Optional list controlling VB warm start. Supported keys:
#' `method` (`"isvb"` or `"ldvb"`), `tol`, `n.IS`, `n.samp`, `max_iter`, `verbose`.
#' @param vb_init_fit Optional precomputed VB fit object. If supplied, warm start
#' uses this object directly and does not rerun VB internally.
#' @param mcmc_control Optional normalized MCMC control list, usually from
#' [exal_make_mcmc_control()]. When supplied, the core MCMC arguments and
#' warmup blocks are read from `mcmc_control` first and then merged with the
#' explicit function arguments. When omitted, exAL-style dynamic MCMC uses the
#' package's conservative default `(sigma, gamma)` warmup profile
#' automatically; explicit controls remain the advanced override path.
#' @param sigmagam_controls Optional list controlling warmup/freeze for the
#' exDQLM sigma/gamma block during MCMC.
#' @param latent_state_controls Optional list controlling early latent-state
#' warmup/freeze in dynamic MCMC. Supported keys include
#' `freeze_burnin_iters`, `freeze_only_during_burn`, `force_after_warmup`,
#' and `mode` (`"u_only"` or `"u_st_pair"`).
#' @param theta_state_controls Optional list controlling early theta-state
#' warmup/freeze in dynamic MCMC. Supported keys include
#' `freeze_burnin_iters`, `freeze_only_during_burn`, and
#' `force_after_warmup`.
#' @param dqlm_sigma_controls Optional list controlling sigma-only
#' warmup/freeze in the DQLM branch. Supported keys mirror
#' `sigmagam_controls`.
#' @param mh.proposal Character; proposal kernel for the exDQLM scale/skew block.
#' `"slice"` (default) uses
#' an exact sigma GIG update plus a bounded univariate slice sampler directly
#' on `gamma`; `"laplace_rw"` uses a Laplace-informed covariance then RW;
#' and `"rw"` uses joint random-walk MH on `(log sigma, logit gamma)`.
#' This choice is separate from the VB warm-start method.
#' @param mh.adapt Logical; adapt MH proposal scale during burn-in.
#' @param mh.adapt.interval Integer; adaptation interval (iterations).
#' @param mh.target.accept Numeric length-2 vector with lower/upper target acceptance rates.
#' @param mh.scale.bounds Numeric length-2 vector with min/max global scaling for MH covariance.
#' @param mh.max_scale.step Numeric in (0,1); maximum fractional scale change per adaptation step.
#' @param mh.min_burn_adapt Minimum burn-in iterations required to enable adaptation.
#' @param slice.width Positive numeric width for the bounded slice sampler when
#' `mh.proposal = "slice"`. Default `0.1` for parity with `bqrgal`.
#' @param slice.max.steps Positive integer or `Inf`; maximum stepping-out
#' expansions for the slice sampler.
#' @param trace.diagnostics Logical; if `TRUE`, retain per-iteration
#' sigma/gamma/s/u diagnostics under `mh.diagnostics$trace`. Set `FALSE` for
#' lighter-weight runs.
#' @param trace.every Positive integer; when `trace.diagnostics = TRUE`, record
#' one diagnostics row every `trace.every` iterations.
#' @param verbose.every Positive integer controlling how often console progress
#' is printed when `verbose = TRUE`. Default `500`, independent of
#' `trace.every`.
#' @param progress_callback Optional callback invoked with a named list at MCMC
#' start, at each progress checkpoint, and on completion. Intended for
#' workflow-level progress logging.
#' @param PriorSigma List of parameters for inverse gamma prior on sigma; shape
#' \code{a_sig} and scale \code{b_sig}. Default is an inverse gamma with mean
#' 1, or \code{sig.init} when supplied, and variance 10.
#' @param PriorGamma List of parameters for truncated Student-t prior on gamma;
#' center \code{m_gam}, scale \code{s_gam}, and degrees of freedom
#' \code{df_gam}. Default is a standard Student-t with 1 degree of freedom,
#' truncated to the support of gamma.
#' @param verbose Logical value indicating whether progress should be displayed.
#'
#' @return An object with classes "\code{exdqlmMCMC}" and "\code{exdqlmFit}"
#' containing the following:
#' \itemize{
#' \item `y` - Time-series data used to fit the model.
#' \item `run.time` - Algorithm run time in seconds.
#' \item `dqlm.ind` - Logical value indicating whether gamma was fixed at `0`, reducing the exDQLM to the special case of the DQLM.
#' \item `model` - List of the state-space model including `GG`, `FF`, prior parameters `m0` and `C0`.
#' \item `p0` - The quantile which was estimated.
#' \item `df` - Discount factors used for each block.
#' \item `dim.df` - Dimension used for each block of discount factors.
#' \item `samp.theta` - Posterior sample of the state vector.
#' \item `samp.post.pred` - Sample of the posterior predictive distributions.
#' \item `map.standard.forecast.errors` - MAP standardized one-step-ahead forecast errors.
#' \item `samp.sigma` - Posterior sample of scale parameter sigma.
#' \item `samp.vts` - Posterior sample of latent parameters, v_t.
#' \item `theta.out` - List containing the distributions of the state vector including filtered distribution parameters (`fm` and `fC`) and smoothed distribution parameters (`sm` and `sC`).
#' \item `n.burn` Number of MCMC iterations that were burned.
#' \item `n.mcmc` Number of MCMC iterations that were sampled.
#' }
#' If `dqlm.ind=FALSE`, the object also contains the following:
#' \itemize{
#' \item `samp.gamma` - Posterior sample of skewness parameter gamma.
#' \item `samp.sts` - Posterior sample of latent parameters, s_t.
#' \item `init.log.sigma` - Burned samples of log sigma from the random walk MH joint sampling of sigma and gamma.
#' \item `init.logit.gamma` - Burned samples of logit gamma from the random walk MH joint sampling of sigma and gamma.
#' \item `accept.rate` - Acceptance rate of the MH step.
#' \item `accept.rate.burn` - MH acceptance rate during burn-in.
#' \item `accept.rate.keep` - MH acceptance rate in kept MCMC samples.
#' \item `Sig.mh` - Covariance matrix used in MH step to jointly sample sigma and gamma.
#' \item `mh.diagnostics` - MH tuning diagnostics (proposal mode, scaling path, adaptation summary).
#' \item `diagnostics` - ESS and chain-ready summaries for sigma/gamma.
#' }
#' @export
#'
#' @examples
#' \donttest{
#' data("scIVTmag", package = "exdqlm")
#' y = scIVTmag[1:80]
#' trend.comp = polytrendMod(order = 1, m0 = stats::quantile(y, 0.85), C0 = 10)
#' seas.comp = seasMod(p = 365, h = c(1,2), C0 = 10*diag(4))
#' model = trend.comp + seas.comp
#' M2 = exdqlmMCMC(y, p0=0.85, model, df = c(1,1), dim.df = c(1,4),
#' gam.init = -3.5, sig.init = 15,
#' n.burn = 40, n.mcmc = 40,
#' init.from.vb = FALSE, verbose = FALSE)
#'
#' M2_al = exdqlmMCMC(y, p0=0.85, model, df = c(1,1), dim.df = c(1,4),
#' dqlm.ind = TRUE, sig.init = 15,
#' n.burn = 30, n.mcmc = 30,
#' init.from.vb = FALSE, verbose = FALSE)
#' }
#'
#' @name exdqlmMCMC
NULL
.exdqlm_latent_state_mcmc_controls <- function(latent_cfg = NULL, default_mode = c("u_only", "u_st_pair")) {
default_mode <- match.arg(default_mode)
latent_cfg <- latent_cfg %||% list()
freeze_burnin_iters <- suppressWarnings(as.integer(
latent_cfg$freeze_burnin_iters %||%
latent_cfg$freeze_latent_burnin_iters %||%
latent_cfg$freeze_ust_burnin_iters %||%
0L
)[1L])
if (!is.finite(freeze_burnin_iters) || freeze_burnin_iters < 0L) freeze_burnin_iters <- 0L
mode <- tolower(as.character(
latent_cfg$mode %||%
latent_cfg$latent_mode %||%
latent_cfg$freeze_mode %||%
default_mode
)[1L])
if (!(mode %in% c("u_only", "u_st_pair"))) mode <- default_mode
min_postwarmup_updates <- suppressWarnings(as.integer(
latent_cfg$min_postwarmup_updates %||%
latent_cfg$latent_min_postwarmup_updates %||%
0L
)[1L])
if (!is.finite(min_postwarmup_updates) || min_postwarmup_updates < 0L) {
min_postwarmup_updates <- 0L
}
list(
mode = mode,
freeze_burnin_iters = freeze_burnin_iters,
freeze_only_during_burn = if (is.null(latent_cfg$freeze_only_during_burn)) TRUE else isTRUE(latent_cfg$freeze_only_during_burn),
force_after_warmup = if (is.null(latent_cfg$force_after_warmup)) TRUE else isTRUE(latent_cfg$force_after_warmup),
min_postwarmup_updates = min_postwarmup_updates,
trace = if (is.null(latent_cfg$trace)) TRUE else isTRUE(latent_cfg$trace)
)
}
.exdqlm_dqlm_sigma_mcmc_controls <- function(sigma_cfg = NULL) {
sigma_cfg <- sigma_cfg %||% list()
freeze_burnin_iters <- suppressWarnings(as.integer(
sigma_cfg$freeze_burnin_iters %||%
sigma_cfg$freeze_sigma_burnin_iters %||%
sigma_cfg$freeze_dqlm_sigma_burnin_iters %||%
0L
)[1L])
if (!is.finite(freeze_burnin_iters) || freeze_burnin_iters < 0L) freeze_burnin_iters <- 0L
list(
freeze_burnin_iters = freeze_burnin_iters,
freeze_only_during_burn = if (is.null(sigma_cfg$freeze_only_during_burn)) TRUE else isTRUE(sigma_cfg$freeze_only_during_burn),
force_after_warmup = if (is.null(sigma_cfg$force_after_warmup)) TRUE else isTRUE(sigma_cfg$force_after_warmup),
trace = if (is.null(sigma_cfg$trace)) TRUE else isTRUE(sigma_cfg$trace)
)
}
.exdqlm_theta_state_mcmc_controls <- function(theta_cfg = NULL) {
theta_cfg <- theta_cfg %||% list()
freeze_burnin_iters <- suppressWarnings(as.integer(
theta_cfg$freeze_burnin_iters %||%
theta_cfg$freeze_theta_burnin_iters %||%
theta_cfg$freeze_theta_state_burnin_iters %||%
0L
)[1L])
if (!is.finite(freeze_burnin_iters) || freeze_burnin_iters < 0L) freeze_burnin_iters <- 0L
list(
freeze_burnin_iters = freeze_burnin_iters,
freeze_only_during_burn = if (is.null(theta_cfg$freeze_only_during_burn)) TRUE else isTRUE(theta_cfg$freeze_only_during_burn),
force_after_warmup = if (is.null(theta_cfg$force_after_warmup)) TRUE else isTRUE(theta_cfg$force_after_warmup)
)
}
.exdqlm_mcmc_warmup_state <- function(iter, ctrl, postwarmup_update_count, n_burn) {
active <- isTRUE(
ctrl$freeze_burnin_iters > 0L &&
iter <= ctrl$freeze_burnin_iters &&
(!ctrl$freeze_only_during_burn || iter <= n_burn)
)
force_now <- isTRUE(
!active &&
ctrl$freeze_burnin_iters > 0L &&
ctrl$force_after_warmup &&
postwarmup_update_count <= 0L
)
reason <- if (isTRUE(active)) {
"warmup"
} else if (isTRUE(force_now)) {
"force_after_warmup"
} else {
"scheduled"
}
list(active = active, force_now = force_now, reason = reason)
}
exdqlmMCMC <- function(y,p0,model,df,dim.df,fix.gamma=FALSE,gam.init=NA,fix.sigma=FALSE,sig.init=NA,dqlm.ind=FALSE,
Sig.mh,joint.sample=FALSE,n.burn=2000,n.mcmc=1500,init.from.isvb=FALSE,PriorSigma=NULL,PriorGamma=NULL,verbose=TRUE,
init.from.vb=TRUE,vb_init_controls=NULL,vb_init_fit=NULL,mcmc_control=NULL,sigmagam_controls=NULL,latent_state_controls=NULL,theta_state_controls=NULL,dqlm_sigma_controls=NULL,
mh.proposal=c("slice","laplace_rw","rw"),mh.adapt=TRUE,mh.adapt.interval=50L,
mh.target.accept=c(0.20,0.45),mh.scale.bounds=c(0.1,10),
mh.max_scale.step=0.35,mh.min_burn_adapt=50L,
slice.width=0.1,slice.max.steps=Inf,
trace.diagnostics=TRUE,trace.every=1L,verbose.every=500L,
progress_callback=NULL){
if (!is.null(mcmc_control)) {
mcmc_control <- exal_make_mcmc_control(control = mcmc_control)
if (!is.null(mcmc_control$n_burn)) n.burn <- as.integer(mcmc_control$n_burn)[1L]
if (!is.null(mcmc_control$n_mcmc)) n.mcmc <- as.integer(mcmc_control$n_mcmc)[1L]
if (!is.null(mcmc_control$verbose)) verbose <- isTRUE(mcmc_control$verbose)
if (!is.null(mcmc_control$progress_every)) verbose.every <- as.integer(mcmc_control$progress_every)[1L]
if (!is.null(mcmc_control$init_from_vb)) init.from.vb <- isTRUE(mcmc_control$init_from_vb)
if (!is.null(mcmc_control$vb_warm_start_control)) {
vb_init_controls <- utils::modifyList(vb_init_controls %||% list(), mcmc_control$vb_warm_start_control)
}
if (!is.null(mcmc_control$sigmagam)) sigmagam_controls <- utils::modifyList(sigmagam_controls %||% list(), mcmc_control$sigmagam)
if (!is.null(mcmc_control$latent_state)) latent_state_controls <- utils::modifyList(latent_state_controls %||% list(), mcmc_control$latent_state)
if (!is.null(mcmc_control$theta)) theta_state_controls <- utils::modifyList(theta_state_controls %||% list(), mcmc_control$theta)
if (!is.null(mcmc_control$dqlm_sigma)) dqlm_sigma_controls <- utils::modifyList(dqlm_sigma_controls %||% list(), mcmc_control$dqlm_sigma)
}
# check inputs
y = check_ts(y)
model = check_mod(model)
rv = check_logics(gam.init,sig.init,fix.gamma,fix.sigma,dqlm.ind)
gam.init = rv$gam.init
dqlm.ind = rv$dqlm.ind
fix.gamma = rv$fix.gamma
### MCMC iterations
if(n.mcmc<=0){
stop("number of mcmc samples must be positive")
}
if(verbose & n.burn<=0){
warning("mcmc will be sampled without burn-in, a burn-in is recommended even if initializing using the isvb algorithm")
n.burn=0
}
I = n.mcmc + n.burn
mh.proposal <- match.arg(mh.proposal)
if (is.null(init.from.vb)) {
init.from.vb <- isTRUE(init.from.isvb)
}
if (!is.null(vb_init_fit)) {
init.from.vb <- TRUE
}
init.from.vb <- isTRUE(init.from.vb)
vb.ctrl.default <- list(
method = if (isTRUE(init.from.isvb)) "isvb" else "ldvb",
tol = 0.5,
n.IS = 200L,
n.samp = 200L,
max_iter = getOption("exdqlm.max_iter", 200L),
verbose = FALSE,
ld_controls = NULL
)
if (is.null(vb_init_controls)) vb_init_controls <- list()
vb.ctrl <- utils::modifyList(vb.ctrl.default, vb_init_controls)
vb.ctrl$method <- tolower(as.character(vb.ctrl$method)[1])
if (!(vb.ctrl$method %in% c("isvb", "ldvb"))) vb.ctrl$method <- "isvb"
vb.ctrl$tol <- as.numeric(vb.ctrl$tol)[1]
if (!is.finite(vb.ctrl$tol) || vb.ctrl$tol <= 0) vb.ctrl$tol <- 0.5
vb.ctrl$n.IS <- suppressWarnings(as.integer(vb.ctrl$n.IS)[1])
if (!is.finite(vb.ctrl$n.IS) || vb.ctrl$n.IS < 20L) vb.ctrl$n.IS <- 200L
vb.ctrl$n.samp <- suppressWarnings(as.integer(vb.ctrl$n.samp)[1])
if (!is.finite(vb.ctrl$n.samp) || vb.ctrl$n.samp < 20L) vb.ctrl$n.samp <- 200L
vb.ctrl$max_iter <- suppressWarnings(as.integer(vb.ctrl$max_iter)[1])
if (!is.finite(vb.ctrl$max_iter) || vb.ctrl$max_iter < 5L) vb.ctrl$max_iter <- 200L
vb.ctrl$verbose <- isTRUE(vb.ctrl$verbose)
if (!is.null(vb.ctrl$ld_controls) && !is.list(vb.ctrl$ld_controls)) {
stop("vb_init_controls$ld_controls must be a list or NULL")
}
mh.adapt <- isTRUE(mh.adapt)
mh.adapt.interval <- suppressWarnings(as.integer(mh.adapt.interval)[1])
if (!is.finite(mh.adapt.interval) || mh.adapt.interval < 5L) mh.adapt.interval <- 50L
mh.min_burn_adapt <- suppressWarnings(as.integer(mh.min_burn_adapt)[1])
if (!is.finite(mh.min_burn_adapt) || mh.min_burn_adapt < 20L) mh.min_burn_adapt <- 50L
if (length(mh.target.accept) != 2L) mh.target.accept <- c(0.20, 0.45)
mh.target.accept <- as.numeric(mh.target.accept)
mh.target.accept <- sort(pmin(pmax(mh.target.accept, 0.01), 0.99))
if (length(mh.scale.bounds) != 2L) mh.scale.bounds <- c(0.1, 10)
mh.scale.bounds <- sort(as.numeric(mh.scale.bounds))
if (!all(is.finite(mh.scale.bounds)) || mh.scale.bounds[1] <= 0 || mh.scale.bounds[2] <= mh.scale.bounds[1]) {
mh.scale.bounds <- c(0.1, 10)
}
mh.max_scale.step <- as.numeric(mh.max_scale.step)[1]
if (!is.finite(mh.max_scale.step) || mh.max_scale.step <= 0 || mh.max_scale.step >= 1) {
mh.max_scale.step <- 0.35
}
mh.laplace.refresh.interval <- suppressWarnings(as.integer(getOption("exdqlm.mcmc.laplace_refresh_interval", mh.adapt.interval))[1])
if (!is.finite(mh.laplace.refresh.interval) || mh.laplace.refresh.interval < 5L) {
mh.laplace.refresh.interval <- mh.adapt.interval
}
mh.laplace.refresh.start <- suppressWarnings(as.integer(getOption("exdqlm.mcmc.laplace_refresh_start", mh.min_burn_adapt))[1])
if (!is.finite(mh.laplace.refresh.start) || mh.laplace.refresh.start < 1L) {
mh.laplace.refresh.start <- mh.min_burn_adapt
}
mh.laplace.refresh.weight <- as.numeric(getOption("exdqlm.mcmc.laplace_refresh_weight", 0.60))[1]
if (!is.finite(mh.laplace.refresh.weight) || mh.laplace.refresh.weight <= 0 || mh.laplace.refresh.weight > 1) {
mh.laplace.refresh.weight <- 0.60
}
sigmagam_ctrl <- .exal_sigmagam_mcmc_controls(sigmagam_controls)
sigmagam_ctrl <- .exal_clamp_mcmc_sigmagam_control(sigmagam_ctrl, n_burn = n.burn)
latent_ctrl <- .exdqlm_latent_state_mcmc_controls(
latent_state_controls,
default_mode = if (isTRUE(dqlm.ind)) "u_only" else "u_st_pair"
)
theta_ctrl <- .exdqlm_theta_state_mcmc_controls(theta_state_controls)
dqlm_sigma_ctrl <- .exdqlm_dqlm_sigma_mcmc_controls(dqlm_sigma_controls)
slice.width <- as.numeric(slice.width)[1]
if (!is.finite(slice.width) || slice.width <= 0) slice.width <- 0.1
slice.max.steps <- as.numeric(slice.max.steps)[1]
if (!(is.infinite(slice.max.steps) || (is.finite(slice.max.steps) && slice.max.steps >= 1 && floor(slice.max.steps) == slice.max.steps))) {
slice.max.steps <- Inf
}
trace.diagnostics <- isTRUE(trace.diagnostics)
trace.every <- suppressWarnings(as.integer(trace.every)[1])
if (!is.finite(trace.every) || trace.every < 1L) trace.every <- 1L
verbose.every <- suppressWarnings(as.integer(verbose.every)[1])
if (!is.finite(verbose.every) || verbose.every < 1L) verbose.every <- 500L
verbose_every_env <- suppressWarnings(as.integer(Sys.getenv("EXDQLM_MCMC_PROGRESS_EVERY", NA_character_))[1])
if (is.finite(verbose_every_env) && !is.na(verbose_every_env) && verbose_every_env >= 1L) {
verbose.every <- verbose_every_env
}
safe_progress_callback <- function(info) {
if (!is.function(progress_callback)) return(invisible(NULL))
try(progress_callback(info), silent = TRUE)
invisible(NULL)
}
state_signal <- function(FF_local, theta_mat) {
drop(colSums(FF_local * theta_mat))
}
warmup_state <- function(iter, ctrl, postwarmup_update_count) {
.exdqlm_mcmc_warmup_state(
iter = iter,
ctrl = ctrl,
postwarmup_update_count = postwarmup_update_count,
n_burn = n.burn
)
}
max_abs_finite <- function(x) {
x <- as.numeric(x)
if (!length(x) || any(!is.finite(x))) return("NA")
format(max(abs(x), na.rm = TRUE), digits = 6)
}
count_nonfinite <- function(x) {
sum(!is.finite(as.numeric(x)))
}
count_nonpositive <- function(x) {
x <- as.numeric(x)
sum(is.finite(x) & x <= 0)
}
if (n.burn < mh.min_burn_adapt) mh.adapt <- FALSE
### Define L and U
bounds = .gamma_bounds(p0)
L = bounds["L"]; U = bounds["U"]
if(!is.na(gam.init)){
if(gam.init < L | gam.init > U){
stop(sprintf("gam.init must be between %s and %s for %s quantile",round(L,3),round(U,3),p0))
}
}
### sigma and gamma priors
# sigma ~ IG(a_sig,b_sig)
if(is.null(PriorSigma)){
m_sigma = 1
v_sigma = 10
PriorSigma$a_sig = (m_sigma^2)/(v_sigma) + 2
PriorSigma$b_sig = (m_sigma^3)/(v_sigma) + m_sigma
}else{
if(!is.list(PriorSigma) | any( is.na( match(c("a_sig", "b_sig"),names(PriorSigma)) ) )){
stop("`PriorSigma` must be a list containing `a_sig` and `b_sig`")
}
}
PriorSigmaDens<-function(sigma){ LaplacesDemon::dinvgamma(sigma,shape=PriorSigma$a_sig,scale=PriorSigma$b_sig) }
# gamma ~ truncated student t on L,U
PriorGamma <- .normalize_gamma_prior_trunc_t(PriorGamma)
PriorGammaDens <- function(gamma) {
.gamma_prior_density_trunc_t(gamma, bounds = c(L, U), PriorGamma = PriorGamma, log = FALSE)
}
### state-space model
## prior, theta ~ N(m0,C0)
m0 = model$m0
C0 = model$C0
#
TT = length(y)
p = length(m0)
if(!is.na(dim(model$GG)[3])){
if(dim(model$GG)[3] != TT){stop("time-varying dimension of GG does not match length of y")}
}
GG = array(model$GG,c(p,p,TT)); model$GG = GG
if(ncol(model$FF)>1){
if(ncol(model$FF) != TT){stop("time-varying dimension of FF does not match length of y")}
}
FF = matrix(model$FF,p,TT); model$FF = FF
## discount factor blocking
if(!methods::hasArg(dim.df)){
if(length(df)!=1){
stop("length of component discount factors does not match length of component dimensions")
}
dim.df = p
}
df.mat = make_df_mat(df,dim.df,p)
### backend controls (MCMC-specific)
use_cpp_mcmc_opt <- isTRUE(getOption("exdqlm.use_cpp_mcmc", FALSE))
cpp_mcmc_mode <- tolower(as.character(getOption("exdqlm.cpp_mcmc_mode", "strict")))
if (!(cpp_mcmc_mode %in% c("strict", "fast"))) {
warning("Invalid exdqlm.cpp_mcmc_mode; using 'strict'.")
cpp_mcmc_mode <- "strict"
}
has_cpp_mcmc <- exists("mcmc_ffbs_smooth_cpp", mode = "function") &&
exists("mcmc_ffbs_sample_cpp", mode = "function")
if (use_cpp_mcmc_opt && !has_cpp_mcmc) {
warning("exdqlm.use_cpp_mcmc=TRUE but C++ MCMC FFBS kernels not available; using R backend.")
}
# strict mode keeps R kernels to preserve exact legacy path; fast enables C++ FFBS.
use_cpp_mcmc <- isTRUE(use_cpp_mcmc_opt && has_cpp_mcmc && identical(cpp_mcmc_mode, "fast"))
mcmc_backend <- if (use_cpp_mcmc) "C++" else "R"
if (verbose) {
cat(sprintf("MCMC backend: %s (mode=%s)\n", mcmc_backend, cpp_mcmc_mode))
}
cpp_ffbs_smooth <- function(ex.f, ex.q) {
out <- mcmc_ffbs_smooth_cpp(
GG = GG,
m0 = as.numeric(m0),
C0 = C0,
FF = FF,
y = as.numeric(y),
ex_f = as.numeric(ex.f),
ex_q = as.numeric(ex.q),
df_mat = df.mat
)
out$standard.forecast.errors <- as.numeric(out$standard.forecast.errors)
out$sm <- as.matrix(out$sm)
out$fm <- as.matrix(out$fm)
out$sC <- array(out$sC, dim = c(p, p, TT))
out$fC <- array(out$fC, dim = c(p, p, TT))
out
}
cpp_ffbs_sample <- function(ex.f, ex.q) {
out <- mcmc_ffbs_sample_cpp(
GG = GG,
m0 = as.numeric(m0),
C0 = C0,
FF = FF,
y = as.numeric(y),
ex_f = as.numeric(ex.f),
ex_q = as.numeric(ex.q),
df_mat = df.mat
)
out$standard.forecast.errors <- as.numeric(out$standard.forecast.errors)
out$sam.theta <- as.matrix(out$sam.theta)
out$fm <- as.matrix(out$fm)
out$fC <- array(out$fC, dim = c(p, p, TT))
out
}
# function to produce smoothed estimates for return value
smoothed_theta<-function(ex.f,ex.q){
# initialize ffbs
m <- sm <- matrix(NA,p,TT)
C <- sC <- array(NA,c(p,p,TT))
standard.forecast.errors <- rep(NA,TT)
## forward filter
# first iteration
a = as.vector(GG[,,1]%*%m0)
P = .exdqlm_regularize_cov(GG[,,1]%*%C0%*%t(GG[,,1]), context = "mcmc_smooth_P_t1")
R = .exdqlm_regularize_cov(P + df.mat*P, context = "mcmc_smooth_R_t1")
f = t(FF[,1])%*%a + ex.f[1]
q = .exdqlm_regularize_var(t(FF[,1])%*%R%*%FF[,1] + ex.q[1], context = "mcmc_smooth_q_t1")
m[,1] = a + t(R)%*%FF[,1]%*%(y[1]-f)/q[1]
C[,,1] = .exdqlm_regularize_cov(
R - t(R)%*%FF[,1]%*%t(FF[,1])%*%R/q[1],
context = "mcmc_smooth_C_t1"
)
standard.forecast.errors[1] = (y[1]-f)/sqrt(q)
# t = 2:TT
for(t in 2:TT){
a = as.vector(GG[,,t]%*%m[,(t-1)])
P = .exdqlm_regularize_cov(GG[,,t]%*%C[,,(t-1)]%*%t(GG[,,t]), context = sprintf("mcmc_smooth_P_t%d", t))
R = .exdqlm_regularize_cov(P + df.mat*P, context = sprintf("mcmc_smooth_R_t%d", t))
f = t(FF[,t])%*%a + ex.f[t]
fB = t(FF[,t])%*%R
q = .exdqlm_regularize_var(fB%*%FF[,t] + ex.q[t], context = sprintf("mcmc_smooth_q_t%d", t))
m[,t] = a + t(fB)%*%(y[t]-f)/q[1]
C[,,t] = .exdqlm_regularize_cov(
R - t(fB)%*%fB/q[1],
context = sprintf("mcmc_smooth_C_t%d", t)
)
standard.forecast.errors[t] = (y[t]-f)/sqrt(q)
}
## backwards smoothing
sC[,,TT] = C[,,TT]
sm[,TT] = m[,TT]
for(t in (TT-1):1){
P = .exdqlm_regularize_cov(GG[,,(t+1)]%*%C[,,(t)]%*%t(GG[,,(t+1)]), context = sprintf("mcmc_smooth_back_P_t%d", t + 1L))
R.info = .exdqlm_cov_inverse(P + df.mat*P, context = sprintf("mcmc_smooth_back_R_t%d", t + 1L))
sB = C[,,t]%*%t(GG[,,(t+1)])%*%R.info$inverse
sm[,t] = m[,t] + sB%*%(sm[,(t+1)]-as.vector(GG[,,(t+1)]%*%m[,(t)]))
sC[,,t] = .exdqlm_regularize_cov(
C[,,t] + sB%*%(sC[,,(t+1)]-R.info$Sigma)%*%t(sB),
context = sprintf("mcmc_smooth_back_C_t%d", t)
)
}
return(list(standard.forecast.errors=standard.forecast.errors,sm=sm,sC=sC,fm=m,fC=C))
}
if (use_cpp_mcmc) {
smoothed_theta <- function(ex.f, ex.q) cpp_ffbs_smooth(ex.f, ex.q)
}
### Initialize MCMC
init.log.sigma <- init.logit.gamma <- rep(NA,n.burn)
save.sigma <- save.gamma <- rep(NA,n.mcmc)
save.Ut <- save.st <- matrix(NA,TT,n.mcmc)
save.theta <- array(NA,c(p,TT,n.mcmc))
save.post.pred <- matrix(NA,TT,n.mcmc)
vb.out <- NULL
gig_backend <- "cpp_devroye_required"
gig_eps <- 1e-12
gig_b_floor <- .gig_b_floor()
current_iter <- NA_integer_
sample_gig_cpp_required <- function(chi, psi, lambda = 0.5, context = "gig") {
if (!exists("sample_gig_devroye_vector", mode = "function")) {
stop(sprintf("%s requires sample_gig_devroye_vector(), but it is not available", context))
}
chi <- as.numeric(chi)
psi <- as.numeric(psi)[1]
lambda <- as.numeric(lambda)[1]
iter_suffix <- if (is.finite(current_iter)) sprintf(" (iter=%d)", current_iter) else ""
bad <- which(!is.finite(chi))
if (length(bad)) {
stop(sprintf("%s%s chi has %d non-finite values (first index=%d)", context, iter_suffix, length(bad), bad[1]))
}
badneg <- which(chi < 0)
if (length(badneg)) {
stop(sprintf("%s%s chi has %d negative values (first index=%d, value=%.6g)", context, iter_suffix, length(badneg), badneg[1], chi[badneg[1]]))
}
if (!is.finite(psi)) {
stop(sprintf("%s%s psi must be finite; got %.6g", context, iter_suffix, psi))
}
if (!is.finite(lambda)) {
stop(sprintf("%s%s lambda must be finite; got %.6g", context, iter_suffix, lambda))
}
chi <- pmax(chi, gig_b_floor)
psi <- max(psi, gig_eps)
draws <- as.numeric(sample_gig_devroye_vector(
1L, p = lambda, a = psi, b_vec = chi
)[1, ])
bad_draws <- which(!is.finite(draws) | draws <= 0)
if (length(bad_draws)) {
stop(sprintf("%s%s sample_gig_devroye_vector returned %d invalid draws (first index=%d, value=%.6g)",
context, iter_suffix, length(bad_draws), bad_draws[1], draws[bad_draws[1]]))
}
pmax(draws, gig_eps)
}
run_vb_init <- function() {
old_opt <- options(exdqlm.max_iter = vb.ctrl$max_iter)
on.exit(options(old_opt), add = TRUE)
if (vb.ctrl$method == "ldvb") {
if (is.list(vb.ctrl$ld_controls) && length(vb.ctrl$ld_controls)) {
ld_opt_names <- paste0("exdqlm.dynamic.ldvb.", names(vb.ctrl$ld_controls))
old_ld_opt <- options(stats::setNames(vb.ctrl$ld_controls, ld_opt_names))
on.exit(options(old_ld_opt), add = TRUE)
}
exdqlmLDVB(
y = y, p0 = p0, model = model, df = df, dim.df = dim.df,
fix.gamma = fix.gamma, gam.init = gam.init,
fix.sigma = fix.sigma, sig.init = sig.init,
dqlm.ind = dqlm.ind,
tol = vb.ctrl$tol, n.samp = vb.ctrl$n.samp,
PriorSigma = PriorSigma, PriorGamma = PriorGamma,
verbose = vb.ctrl$verbose
)
} else {
exdqlmISVB(
y = y, p0 = p0, model = model, df = df, dim.df = dim.df,
fix.gamma = fix.gamma, gam.init = gam.init,
fix.sigma = fix.sigma, sig.init = sig.init,
dqlm.ind = dqlm.ind,
tol = vb.ctrl$tol, n.IS = vb.ctrl$n.IS, n.samp = vb.ctrl$n.samp,
PriorSigma = PriorSigma, PriorGamma = PriorGamma,
verbose = vb.ctrl$verbose
)
}
}
# Set initial values
if(init.from.vb){
if(verbose){
cat(sprintf("running %s algorithm to initialize mcmc\n", toupper(vb.ctrl$method)))
}
if (!is.null(vb_init_fit)) {
vb.out <- vb_init_fit
if (verbose) {
cat("using provided vb_init_fit object for MCMC initialization\n")
}
} else {
vb.out <- run_vb_init()
}
cursam.sigma <- ifelse(fix.sigma,sig.init,ifelse(dqlm.ind,vb.out$sig.out$E.sigma,vb.out$gammasig.out$E.sigma))
cursam.Ut <- vb.out$vts.out$E.uts
cursam.theta <- vb.out$theta.out$sm
}else{
cursam.sigma <- m_sigma
cursam.Ut <- rep(1/m_sigma,TT)
cursam.theta <- matrix(m0,p,TT)
}
if (verbose) {
cat("GIG backend: C++ Devroye (required)\n")
}
######## exDQLM
if(!dqlm.ind){
### Define logit and inverse logit functions
logit = function(x){log((x-L)/(U-x))}
inv.logit = function(x){(U*exp(x)+L)/(exp(x)+1)}
log_prior_gamma <- function(gamma) {
.gamma_log_prior_trunc_t(gamma, bounds = c(L, U), PriorGamma = PriorGamma)
}
### Additional initial values
if(!is.null(vb.out)){
cursam.st <- vb.out$sts.out$E.sts
cursam.gamma <- ifelse(fix.gamma,gam.init,vb.out$gammasig.out$E.gam)
cursam.logit.gamma <- logit(cursam.gamma)
cursam.log.sigma <- log(cursam.sigma)
}else{
cursam.st <- truncnorm::rtruncnorm(TT,a=0,b=Inf,mean=0,sd=1)
cursam.gamma <- ifelse(!is.na(gam.init),gam.init,(L+U)/2)
cursam.logit.gamma <- logit(cursam.gamma)
cursam.log.sigma <- log(cursam.sigma)
}
### Initialize MH
n.accept = 0
n.accept.burn = 0L
n.accept.keep = 0L
n.trial.burn = 0L
n.trial.keep = 0L
adapt.history <- data.frame(
iter = integer(0),
window_accept = numeric(0),
mh_scale = numeric(0),
sig11 = numeric(0),
sig22 = numeric(0),
laplace_refreshed = logical(0),
stringsAsFactors = FALSE
)
trace_rows <- if (trace.diagnostics) vector("list", ceiling(I / trace.every)) else NULL
trace_idx <- 0L
mh.scale <- 1
window.accept <- 0L
window.total <- 0L
laplace_refresh_attempts <- 0L
laplace_refresh_success <- 0L
prep_Sig_mh <- function(S) {
S <- suppressWarnings(as.matrix(S))
if (!all(dim(S) == c(2L, 2L))) {
S <- diag(c(ifelse(fix.sigma, 0, 0.005), ifelse(fix.gamma, 0, 0.005)))
}
S[!is.finite(S)] <- 0
S <- (S + t(S)) / 2
if (fix.sigma) {
S[1, ] <- 0
S[, 1] <- 0
}
if (fix.gamma) {
S[2, ] <- 0
S[, 2] <- 0
}
for (j in 1:2) {
if (!is.finite(S[j, j]) || S[j, j] < 0) S[j, j] <- 0
}
if (!fix.sigma && S[1, 1] <= 0) S[1, 1] <- 0.005
if (!fix.gamma && S[2, 2] <= 0) S[2, 2] <- 0.005
S
}
build_chol <- function(S) {
S <- prep_Sig_mh(S)
if (fix.gamma || fix.sigma) {
sqrt(S)
} else {
out <- tryCatch(t(chol(S)), error = function(e) NULL)
if (is.null(out)) {
eig <- eigen(S, symmetric = TRUE)
vals <- pmax(eig$values, 1e-8)
out <- eig$vectors %*% diag(sqrt(vals), 2, 2) %*% t(eig$vectors)
}
out
}
}
if(!methods::hasArg(Sig.mh)){
if(!is.null(vb.out)){
sig.samples <- NULL
gam.samples <- NULL
if (!is.null(vb.out$gammasig.out$sigma.samples) && !is.null(vb.out$gammasig.out$gamma.samples)) {
sig.samples <- as.numeric(vb.out$gammasig.out$sigma.samples)
gam.samples <- as.numeric(vb.out$gammasig.out$gamma.samples)
} else if (!is.null(vb.out$samp.sigma) && !is.null(vb.out$samp.gamma)) {
sig.samples <- as.numeric(vb.out$samp.sigma)
gam.samples <- as.numeric(vb.out$samp.gamma)
}
if (!is.null(sig.samples) && !is.null(gam.samples) &&
all(is.finite(sig.samples)) && all(sig.samples > 0) &&
all(is.finite(gam.samples)) && all(gam.samples > L) && all(gam.samples < U)) {
Sig.mh <- stats::cov(cbind(log(sig.samples), logit(gam.samples)))
} else {
Sig.mh <- diag(c(ifelse(fix.sigma,0,0.005),ifelse(fix.gamma,0,0.005)))
}
}else{
Sig.mh = diag(c(ifelse(fix.sigma,0,0.005),ifelse(fix.gamma,0,0.005)))
}
}
Sig.mh <- prep_Sig_mh(Sig.mh)
Sig.mh.initial <- Sig.mh
chol_Sig.mh <- build_chol(Sig.mh)
# exdqlm function sample theta ffbs
ex_samp_theta<-function(ex.f,ex.q,gamma,sigma,sts,tau,c_tau){
# initialize ffbs
m <- sam.theta <- matrix(NA,p,TT)
C <- array(NA,c(p,p,TT))
standard.forecast.errors <- post.pred <- rep(NA,TT)
## forward filter
# first iteration
a = as.vector(GG[,,1]%*%m0)
P = .exdqlm_regularize_cov(GG[,,1]%*%C0%*%t(GG[,,1]), context = "mcmc_dqlm_sample_P_t1")
R = .exdqlm_regularize_cov(P + df.mat*P, context = "mcmc_dqlm_sample_R_t1")
f = t(FF[,1])%*%a + ex.f[1]
q = .exdqlm_regularize_var(t(FF[,1])%*%R%*%FF[,1] + ex.q[1], context = "mcmc_dqlm_sample_q_t1")
m[,1] = a + t(R)%*%FF[,1]%*%(y[1]-f)/q[1]
C[,,1] = .exdqlm_regularize_cov(
R - t(R)%*%FF[,1]%*%t(FF[,1])%*%R/q[1],
context = "mcmc_dqlm_sample_C_t1"
)
standard.forecast.errors[1] = (y[1]-f)/sqrt(q)
# t = 2:TT
for(t in 2:TT){
a = as.vector(GG[,,t]%*%m[,(t-1)])
P = .exdqlm_regularize_cov(GG[,,t]%*%C[,,(t-1)]%*%t(GG[,,t]), context = sprintf("mcmc_dqlm_sample_P_t%d", t))
R = .exdqlm_regularize_cov(P + df.mat*P, context = sprintf("mcmc_dqlm_sample_R_t%d", t))
f = t(FF[,t])%*%a + ex.f[t]
fB = t(FF[,t])%*%R
q = .exdqlm_regularize_var(fB%*%FF[,t] + ex.q[t], context = sprintf("mcmc_dqlm_sample_q_t%d", t))
m[,t] = a + t(fB)%*%(y[t]-f)/q[1]
C[,,t] = .exdqlm_regularize_cov(
R - t(fB)%*%fB/q[1],
context = sprintf("mcmc_dqlm_sample_C_t%d", t)
)
standard.forecast.errors[t] = (y[t]-f)/sqrt(q)
}
## backwards sample
sC_TT = .exdqlm_regularize_cov(C[,,TT], context = "mcmc_dqlm_sample_sC_TT")
svd.sC = svd(sC_TT)
sam.theta[,TT] = m[,TT] + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta <- numeric(TT)
reg_theta[TT] <- drop(crossprod(FF[,TT], sam.theta[,TT]))
post.pred[TT] = rexal(1,tau,reg_theta[TT]+c_tau*sigma*abs(gamma)*sts[TT],sigma,0)
for(t in (TT-1):1){
P = .exdqlm_regularize_cov(GG[,,(t+1)]%*%C[,,(t)]%*%t(GG[,,(t+1)]), context = sprintf("mcmc_dqlm_back_P_t%d", t + 1L))
R.info = .exdqlm_cov_inverse(P + df.mat*P, context = sprintf("mcmc_dqlm_back_R_t%d", t + 1L))
sB = C[,,t]%*%t(GG[,,(t+1)])%*%R.info$inverse
sm = m[,t] + sB%*%(sam.theta[,(t+1)]-as.vector(GG[,,(t+1)]%*%m[,(t)]))
sC = .exdqlm_regularize_cov(
C[,,t] - sB%*%GG[,,(t+1)]%*%C[,,t],
context = sprintf("mcmc_dqlm_back_sC_t%d", t)
)
svd.sC = svd(sC)
sam.theta[,t] = sm + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta[t] <- drop(crossprod(FF[,t], sam.theta[,t]))
post.pred[t] = rexal(1,tau,reg_theta[t]+c_tau*sigma*abs(gamma)*sts[t],sigma,0)
}
return(list(standard.forecast.errors=standard.forecast.errors,post.pred=post.pred,sam.theta=sam.theta,fm=m,fC=C))
}
if (use_cpp_mcmc) {
ex_samp_theta <- function(ex.f, ex.q, gamma, sigma, sts, tau, c_tau) {
out <- cpp_ffbs_sample(ex.f, ex.q)
sam.theta <- out$sam.theta
reg_theta <- state_signal(FF, sam.theta)
post.pred <- vapply(seq_len(TT), function(t) {
rexal(1, tau,
reg_theta[t] + c_tau * sigma * abs(gamma) * sts[t],
sigma, 0)
}, numeric(1))
list(
standard.forecast.errors = out$standard.forecast.errors,
post.pred = post.pred,
sam.theta = sam.theta,
fm = out$fm,
fC = out$fC
)
}
}
# exdqlm function sample uts
ex_samp_uts<-function(reg1,gamma,sigma,sts,a_tau,b_tau,c_tau){
chi <- as.numeric(((y-reg1-sigma*c_tau*abs(gamma)*sts)^2)/(b_tau*sigma))
psi <- (a_tau^2)/(b_tau*sigma) + (2/sigma)
sample_gig_cpp_required(chi = chi, psi = psi, lambda = 0.5, context = "exdqlm_mcmc_uts")
}
# exdqlm function sample sts
ex_samp_sts<-function(reg1,gamma,sigma,uts,a_tau,b_tau,c_tau){
s.sig2<-1/(1+c_tau^2*abs(gamma)^2*sigma/(b_tau*uts))
s.sig2<-pmax(s.sig2, 1e-12)
s.mu<-s.sig2*c_tau*abs(gamma)*(y-(reg1+a_tau*uts))/(b_tau*uts)
truncnorm::rtruncnorm(TT,rep(0,TT),rep(Inf,TT),s.mu,sqrt(s.sig2))
}
# exdqlm function sample sigma and gamma
logL<-function(reg1,log.sigma,logit.gamma,sts,uts){
sigma=exp(log.sigma); gamma=inv.logit(logit.gamma)
temp.p = p.fn(p0,gamma)
a = (1-2*temp.p)/(temp.p*(1-temp.p))
b = (2)/(temp.p*(1-temp.p))
c = (as.numeric(gamma>0)-temp.p)^(-1)
logJ<-logit.gamma-2*log(1+exp(logit.gamma))+log.sigma
PriorGamma<- PriorGammaDens(gamma)
PriorSigma<- PriorSigmaDens(sigma)
sum(stats::dnorm(y,reg1+sigma*c*abs(gamma)*sts+a*uts,sqrt(sigma*b*uts),log = TRUE)) +
sum(stats::dexp(uts,rate = 1/sigma,log=TRUE)) +
log(PriorSigma) + log(PriorGamma) + logJ
}
make_logpost_gamma <- function(reg1, sigma, sts, uts) {
sigma <- as.numeric(sigma)[1]
reg1 <- as.numeric(reg1)
sts <- as.numeric(sts)
uts <- as.numeric(uts)
valid_inputs <- is.finite(sigma) && sigma > 0 &&
all(is.finite(reg1)) && all(is.finite(sts)) && all(is.finite(uts)) &&
all(uts > 0)
if (!valid_inputs) {
return(function(gamma) -Inf)
}
y_center <- y - reg1
sigma_sts <- sigma * sts
sqrt_sigma_uts <- sqrt(sigma * uts)
function(gamma) {
gamma <- as.numeric(gamma)[1]
if (!is.finite(gamma) || gamma <= L || gamma >= U) return(-Inf)
temp.p <- p.fn(p0, gamma)
a <- (1 - 2 * temp.p) / (temp.p * (1 - temp.p))
b <- 2 / (temp.p * (1 - temp.p))
c <- (as.numeric(gamma > 0) - temp.p)^(-1)
if (!all(is.finite(c(a, b, c))) || b <= 0) return(-Inf)
mu_shift <- c * abs(gamma) * sigma_sts + a * uts
ll <- sum(stats::dnorm(
y_center,
mean = mu_shift,
sd = sqrt(b) * sqrt_sigma_uts,
log = TRUE
))
lp <- log_prior_gamma(gamma)
if (!is.finite(ll) || !is.finite(lp)) return(-Inf)
ll + lp
}
}
samp_sigma_exact <- function(reg1, sigma, gamma, sts, uts) {
temp.p <- p.fn(p0, gamma)
a <- (1 - 2 * temp.p) / (temp.p * (1 - temp.p))
b <- 2 / (temp.p * (1 - temp.p))
c <- (as.numeric(gamma > 0) - temp.p)^(-1)
r <- y - reg1 - a * uts
chi_sigma <- sum((r * r) / (b * uts)) + 2 * sum(uts) + 2 * PriorSigma$b_sig
psi_sigma <- ((c * abs(gamma))^2 / b) * sum((sts * sts) / uts)
k_sigma <- -(PriorSigma$a_sig + 1.5 * TT)
sigma_new <- sample_gig_cpp_required(
chi = chi_sigma,
psi = psi_sigma,
lambda = k_sigma,
context = "exdqlm_mcmc_sigma"
)[1]
if (is.finite(sigma_new) && sigma_new > 0) sigma_new else sigma
}
laplace_cov_init <- function(reg1, log.sigma, logit.gamma, sts, uts) {
fn <- function(z) {
val <- logL(reg1, z[1], z[2], sts, uts)
if (is.finite(val)) -val else 1e12
}
opt <- tryCatch(
stats::optim(c(log.sigma, logit.gamma), fn = fn, method = "BFGS",
control = list(maxit = 100), hessian = TRUE),
error = function(e) NULL
)
H <- if (!is.null(opt) && !is.null(opt$hessian)) opt$hessian else NULL
if (is.null(H) || any(!is.finite(H))) {
H <- tryCatch(numDeriv::hessian(fn, x = c(log.sigma, logit.gamma)), error = function(e) NULL)
}
if (is.null(H) || any(!is.finite(H))) return(NULL)
H <- (H + t(H)) / 2
eig <- eigen(H, symmetric = TRUE)
vals <- pmax(eig$values, 1e-6)
cov <- eig$vectors %*% diag(1 / vals, 2, 2) %*% t(eig$vectors)
cov
}
if (identical(mh.proposal, "laplace_rw") && !fix.gamma && !fix.sigma) {
reg1.init <- state_signal(FF, cursam.theta)
cov.lap <- laplace_cov_init(reg1.init, cursam.log.sigma, cursam.logit.gamma, cursam.st, cursam.Ut)
if (!is.null(cov.lap)) {
Sig.mh <- prep_Sig_mh(cov.lap)
}
}
chol_Sig.mh <- build_chol(Sig.mh * (mh.scale^2))
ex_samp_lsiglgam<-function(reg1,log.sigma,logit.gamma,sts,uts,chol_Sig){
prop<-c(log.sigma,logit.gamma)+chol_Sig%*%stats::rnorm(2)
if(inv.logit(prop[2]) < U && inv.logit(prop[2]) > L){
logr<-logL(reg1,prop[1],prop[2],sts,uts)-logL(reg1,log.sigma,logit.gamma,sts,uts)
accept=(log(stats::runif(1))<logr)
}else{
accept=FALSE
}
log.sigma.new<-accept*prop[1]+(1-accept)*log.sigma
logit.gamma.new<-accept*prop[2]+(1-accept)*logit.gamma
return(list(log.sigma=log.sigma.new,logit.gamma=logit.gamma.new,accept=accept))
}
callback.every <- if (trace.diagnostics) {
trace.every
} else {
100L
}
callback.every <- max(1L, as.integer(callback.every)[1])
sigmagam_frozen_trace <- rep(FALSE, I)
sigmagam_forced_postwarmup_trace <- rep(FALSE, I)
sigmagam_update_performed_trace <- rep(FALSE, I)
sigmagam_update_reason_trace <- rep(NA_character_, I)
sigmagam_update_count_trace <- integer(I)
sigmagam_first_active_iter <- NA_integer_
sigmagam_update_count <- 0L
sigmagam_postwarmup_update_count <- 0L
sigmagam_updates_burn <- 0L
sigmagam_updates_keep <- 0L
latent_frozen_trace <- rep(FALSE, I)
latent_forced_postwarmup_trace <- rep(FALSE, I)
latent_update_performed_trace <- rep(FALSE, I)
latent_update_reason_trace <- rep(NA_character_, I)
latent_update_count_trace <- integer(I)
latent_first_active_iter <- NA_integer_
latent_update_count <- 0L
latent_postwarmup_update_count <- 0L
latent_updates_burn <- 0L
latent_updates_keep <- 0L
theta_frozen_trace <- rep(FALSE, I)
theta_forced_postwarmup_trace <- rep(FALSE, I)
theta_update_performed_trace <- rep(FALSE, I)
theta_update_reason_trace <- rep(NA_character_, I)
theta_update_count_trace <- integer(I)
theta_first_active_iter <- NA_integer_
theta_update_count <- 0L
theta_postwarmup_update_count <- 0L
theta_updates_burn <- 0L
theta_updates_keep <- 0L
# Sample from exdqlm posterior
tictoc::tic()
.exdqlm_progress(
"MCMC start",
burn = n.burn,
keep = n.mcmc,
kernel = mh.proposal,
warm_start = if (isTRUE(init.from.vb)) vb.ctrl$method else "none",
.verbose = verbose
)
safe_progress_callback(list(
event = "start",
iter = 0L,
total_iter = as.integer(I),
phase = "burn",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = cursam.gamma,
kernel = mh.proposal,
accept = if (identical(mh.proposal, "slice")) NA_real_ else 0
))
for (i in 1:I){
current_iter <- as.integer(i)
# counter
if (i %% verbose.every == 0L) {
kept_now <- if (i <= n.burn) 0L else (i - n.burn)
.exdqlm_progress(
"MCMC progress",
model = "exDQLM",
phase = if (i <= n.burn) "burn" else "keep",
iter = sprintf("%d/%d", i, I),
accept = if (identical(mh.proposal, "slice")) NULL else n.accept / i,
kept = sprintf("%d/%d", kept_now, n.mcmc),
.verbose = verbose
)
}
if (i %% callback.every == 0L) {
safe_progress_callback(list(
event = "progress",
iter = as.integer(i),
total_iter = as.integer(I),
phase = if (i <= n.burn) "burn" else "keep",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = cursam.gamma,
kernel = mh.proposal,
accept = if (identical(mh.proposal, "slice")) NA_real_ else n.accept / i
))
}
# exAL parameters
tau = p.fn(p0,cursam.gamma)
a_tau = (1-2*tau)/(tau*(1-tau))
b_tau = (2)/(tau*(1-tau))
c_tau = (as.numeric(cursam.gamma>0)-tau)^(-1)
# sample theta
ex.f = cursam.sigma*c_tau*abs(cursam.gamma)*cursam.st + cursam.Ut*a_tau
ex.q = b_tau*cursam.Ut*cursam.sigma
theta_state_sched <- warmup_state(i, theta_ctrl, theta_postwarmup_update_count)
if (isTRUE(theta_state_sched$active)) {
if (any(!is.finite(cursam.theta))) {
stop(sprintf(
"exdqlm_mcmc_theta_warmup_state (iter=%d) invalid warmup theta state: nonfinite=%d",
i, sum(!is.finite(cursam.theta))
))
}
reg_theta_hold <- state_signal(FF, cursam.theta)
theta.out <- list(
standard.forecast.errors = rep(NA_real_, TT),
post.pred = vapply(seq_len(TT), function(t) {
rexal(
1, tau,
reg_theta_hold[t] + c_tau * cursam.sigma * abs(cursam.gamma) * cursam.st[t],
cursam.sigma, 0
)
}, numeric(1)),
sam.theta = cursam.theta,
fm = matrix(NA_real_, p, TT),
fC = array(NA_real_, c(p, p, TT))
)
} else {
theta.out <- ex_samp_theta(ex.f,ex.q,cursam.gamma,cursam.sigma,cursam.st,tau,c_tau)
cursam.theta = theta.out$sam.theta
}
theta_frozen_trace[i] <- isTRUE(theta_state_sched$active)
theta_update_reason_trace[i] <- theta_state_sched$reason
if (!isTRUE(theta_state_sched$active)) {
theta_update_performed_trace[i] <- TRUE
theta_update_count <- theta_update_count + 1L
theta_update_count_trace[i] <- theta_update_count
if (is.na(theta_first_active_iter)) theta_first_active_iter <- as.integer(i)
if (i <= n.burn) {
theta_updates_burn <- theta_updates_burn + 1L
} else {
theta_updates_keep <- theta_updates_keep + 1L
}
if (theta_ctrl$freeze_burnin_iters > 0L) {
theta_postwarmup_update_count <- theta_postwarmup_update_count + 1L
}
theta_forced_postwarmup_trace[i] <- isTRUE(theta_state_sched$force_now)
} else {
theta_update_count_trace[i] <- theta_update_count
theta_forced_postwarmup_trace[i] <- FALSE
}
# sample uts, sts
reg1 = state_signal(FF, cursam.theta)
latent_state_sched <- warmup_state(i, latent_ctrl, latent_postwarmup_update_count)
if (!is.finite(cursam.sigma) || cursam.sigma <= 0 || !is.finite(cursam.gamma) ||
cursam.gamma <= L || cursam.gamma >= U || any(!is.finite(reg1)) || any(!is.finite(cursam.theta))) {
stop(sprintf(
paste(
"exdqlm_mcmc_pre_latent (iter=%d mode=%s reason=%s) invalid state before Ut/st update:",
"sigma=%s gamma=%s reg1_finite=%s theta_finite=%s st_finite=%s ut_finite=%s",
"max_abs_reg1=%s max_abs_theta=%s"
),
i,
latent_ctrl$mode,
latent_state_sched$reason,
format(cursam.sigma, digits = 6),
format(cursam.gamma, digits = 6),
all(is.finite(reg1)),
all(is.finite(cursam.theta)),
all(is.finite(cursam.st)),
all(is.finite(cursam.Ut)),
max_abs_finite(reg1),
max_abs_finite(cursam.theta)
))
}
if (isTRUE(latent_state_sched$active)) {
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"exdqlm_mcmc_latent_warmup_state (iter=%d mode=%s) invalid warmup Ut state: nonfinite=%d nonpositive=%d",
i, latent_ctrl$mode, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut)
))
}
if (identical(latent_ctrl$mode, "u_only")) {
cursam.st <- ex_samp_sts(reg1, cursam.gamma, cursam.sigma, cursam.Ut, a_tau, b_tau, c_tau)
if (count_nonfinite(cursam.st) > 0L || count_nonpositive(cursam.st) > 0L) {
stop(sprintf(
"exdqlm_mcmc_post_sts (iter=%d mode=%s) invalid st draws after Ut warmup: nonfinite=%d nonpositive=%d",
i, latent_ctrl$mode, count_nonfinite(cursam.st), count_nonpositive(cursam.st)
))
}
} else if (count_nonfinite(cursam.st) > 0L || count_nonpositive(cursam.st) > 0L) {
stop(sprintf(
"exdqlm_mcmc_latent_warmup_state (iter=%d mode=%s) invalid warmup st state: nonfinite=%d nonpositive=%d",
i, latent_ctrl$mode, count_nonfinite(cursam.st), count_nonpositive(cursam.st)
))
}
} else {
cursam.Ut <- ex_samp_uts(reg1, cursam.gamma, cursam.sigma, cursam.st, a_tau, b_tau, c_tau)
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"exdqlm_mcmc_post_uts (iter=%d mode=%s) invalid Ut draws: nonfinite=%d nonpositive=%d sigma=%s gamma=%s max_abs_reg1=%s",
i, latent_ctrl$mode, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut),
format(cursam.sigma, digits = 6), format(cursam.gamma, digits = 6), max_abs_finite(reg1)
))
}
cursam.st <- ex_samp_sts(reg1, cursam.gamma, cursam.sigma, cursam.Ut, a_tau, b_tau, c_tau)
if (count_nonfinite(cursam.st) > 0L || count_nonpositive(cursam.st) > 0L) {
stop(sprintf(
"exdqlm_mcmc_post_sts (iter=%d mode=%s) invalid st draws: nonfinite=%d nonpositive=%d sigma=%s gamma=%s max_abs_reg1=%s",
i, latent_ctrl$mode, count_nonfinite(cursam.st), count_nonpositive(cursam.st),
format(cursam.sigma, digits = 6), format(cursam.gamma, digits = 6), max_abs_finite(reg1)
))
}
}
latent_frozen_trace[i] <- isTRUE(latent_state_sched$active)
latent_update_reason_trace[i] <- latent_state_sched$reason
if (!isTRUE(latent_state_sched$active)) {
latent_update_performed_trace[i] <- TRUE
latent_update_count <- latent_update_count + 1L
latent_update_count_trace[i] <- latent_update_count
if (is.na(latent_first_active_iter)) latent_first_active_iter <- as.integer(i)
if (i <= n.burn) {
latent_updates_burn <- latent_updates_burn + 1L
} else {
latent_updates_keep <- latent_updates_keep + 1L
}
if (latent_ctrl$freeze_burnin_iters > 0L) {
latent_postwarmup_update_count <- latent_postwarmup_update_count + 1L
}
latent_forced_postwarmup_trace[i] <- isTRUE(latent_state_sched$force_now)
} else {
latent_update_count_trace[i] <- latent_update_count
latent_forced_postwarmup_trace[i] <- FALSE
}
# sample sigma and gamma
sigmagam_sched <- warmup_state(i, sigmagam_ctrl, sigmagam_postwarmup_update_count)
sigmagam_warmup_active <- isTRUE(sigmagam_sched$active)
sigmagam_force_now <- isTRUE(sigmagam_sched$force_now)
sigmagam_update_reason <- sigmagam_sched$reason
if (isTRUE(sigmagam_warmup_active)) {
lsiglgam.out <- list(
log.sigma = cursam.log.sigma,
logit.gamma = cursam.logit.gamma,
accept = NA,
slice_evals = NA_integer_
)
} else if (identical(mh.proposal, "slice")) {
if (!fix.sigma) {
cursam.sigma <- samp_sigma_exact(reg1, cursam.sigma, cursam.gamma, cursam.st, cursam.Ut)
cursam.log.sigma <- log(cursam.sigma)
}
gamma_log_density <- make_logpost_gamma(reg1, cursam.sigma, cursam.st, cursam.Ut)
slice_evals <- NA_integer_
if (!fix.gamma) {
current_lp <- gamma_log_density(cursam.gamma)
if (!is.finite(current_lp)) {
cursam.gamma <- min(max(cursam.gamma, L + 1e-8), U - 1e-8)
current_lp <- gamma_log_density(cursam.gamma)
}
if (!is.finite(current_lp)) {
cursam.gamma <- min(max(0, L + 1e-8), U - 1e-8)
}
slice_out <- .exdqlm_uni_slice_bounded(
x0 = cursam.gamma,
log_density = gamma_log_density,
w = slice.width,
m = slice.max.steps,
lower = L + 1e-10,
upper = U - 1e-10
)
cursam.gamma <- as.numeric(slice_out$value)[1]
slice_evals <- as.integer(slice_out$evals)
}
cursam.logit.gamma <- logit(cursam.gamma)
lsiglgam.out <- list(
log.sigma = cursam.log.sigma,
logit.gamma = cursam.logit.gamma,
accept = NA,
slice_evals = slice_evals
)
} else {
lsiglgam.out <- ex_samp_lsiglgam(reg1, cursam.log.sigma, cursam.logit.gamma, cursam.st, cursam.Ut, chol_Sig.mh)
cursam.gamma <- inv.logit(lsiglgam.out$logit.gamma)
cursam.logit.gamma <- lsiglgam.out$logit.gamma
cursam.sigma <- exp(lsiglgam.out$log.sigma)
cursam.log.sigma <- lsiglgam.out$log.sigma
n.accept <- n.accept + lsiglgam.out$accept
}
sigmagam_frozen_trace[i] <- isTRUE(sigmagam_warmup_active)
sigmagam_update_reason_trace[i] <- sigmagam_update_reason
if (!isTRUE(sigmagam_warmup_active)) {
sigmagam_update_performed_trace[i] <- TRUE
sigmagam_update_count <- sigmagam_update_count + 1L
sigmagam_update_count_trace[i] <- sigmagam_update_count
if (is.na(sigmagam_first_active_iter)) sigmagam_first_active_iter <- as.integer(i)
if (i <= n.burn) {
sigmagam_updates_burn <- sigmagam_updates_burn + 1L
} else {
sigmagam_updates_keep <- sigmagam_updates_keep + 1L
}
if (sigmagam_ctrl$freeze_burnin_iters > 0L) {
sigmagam_postwarmup_update_count <- sigmagam_postwarmup_update_count + 1L
}
sigmagam_forced_postwarmup_trace[i] <- isTRUE(sigmagam_force_now)
} else {
sigmagam_update_count_trace[i] <- sigmagam_update_count
sigmagam_forced_postwarmup_trace[i] <- FALSE
}
if (trace.diagnostics && (i %% trace.every == 0L)) {
s_stats <- .exdqlm_trace_summary(cursam.st)
u_stats <- .exdqlm_trace_summary(cursam.Ut)
trace_idx <- trace_idx + 1L
trace_rows[[trace_idx]] <- data.frame(
iter = i,
phase = if (i <= n.burn) "burn" else "keep",
sigma = cursam.sigma,
gamma = cursam.gamma,
accepted = if (identical(mh.proposal, "slice")) NA else isTRUE(lsiglgam.out$accept),
mh_scale = if (identical(mh.proposal, "slice")) NA_real_ else mh.scale,
latent_frozen = isTRUE(latent_state_sched$active),
latent_mode = latent_ctrl$mode,
latent_update_reason = latent_state_sched$reason,
latent_forced_postwarmup = isTRUE(latent_forced_postwarmup_trace[i]),
latent_update_performed = isTRUE(latent_update_performed_trace[i]),
latent_update_count = as.integer(latent_update_count_trace[i]),
theta_frozen = isTRUE(theta_state_sched$active),
theta_update_reason = theta_state_sched$reason,
theta_forced_postwarmup = isTRUE(theta_forced_postwarmup_trace[i]),
theta_update_performed = isTRUE(theta_update_performed_trace[i]),
theta_update_count = as.integer(theta_update_count_trace[i]),
sigmagam_frozen = isTRUE(sigmagam_warmup_active),
sigmagam_update_reason = sigmagam_update_reason,
sigmagam_forced_postwarmup = isTRUE(sigmagam_forced_postwarmup_trace[i]),
sigmagam_update_performed = isTRUE(sigmagam_update_performed_trace[i]),
sigmagam_update_count = as.integer(sigmagam_update_count_trace[i]),
slice_evals = if (!is.null(lsiglgam.out$slice_evals)) lsiglgam.out$slice_evals else NA_integer_,
s_mean = s_stats[["mean"]],
s_sd = s_stats[["sd"]],
s_q05 = s_stats[["q05"]],
s_q50 = s_stats[["median"]],
s_q95 = s_stats[["q95"]],
s_min = s_stats[["min"]],
s_max = s_stats[["max"]],
u_mean = u_stats[["mean"]],
u_sd = u_stats[["sd"]],
u_q05 = u_stats[["q05"]],
u_q50 = u_stats[["median"]],
u_q95 = u_stats[["q95"]],
u_min = u_stats[["min"]],
u_max = u_stats[["max"]],
stringsAsFactors = FALSE
)
}
# save samples after burn
if(i <= n.burn){
if (!identical(mh.proposal, "slice") && !isTRUE(sigmagam_warmup_active)) {
n.trial.burn <- n.trial.burn + 1L
n.accept.burn <- n.accept.burn + as.integer(isTRUE(lsiglgam.out$accept))
window.accept <- window.accept + as.integer(isTRUE(lsiglgam.out$accept))
window.total <- window.total + 1L
}
init.log.sigma[i] = cursam.log.sigma
init.logit.gamma[i] = cursam.logit.gamma
laplace_refreshed <- FALSE
if (identical(mh.proposal, "laplace_rw") && !fix.gamma && !fix.sigma &&
(!isTRUE(sigmagam_ctrl$delay_laplace_refresh_until_after_warmup) || !isTRUE(sigmagam_warmup_active)) &&
i >= mh.laplace.refresh.start && i < n.burn &&
(i %% mh.laplace.refresh.interval == 0)) {
laplace_refresh_attempts <- laplace_refresh_attempts + 1L
cov.lap.step <- laplace_cov_init(reg1, cursam.log.sigma, cursam.logit.gamma, cursam.st, cursam.Ut)
if (!is.null(cov.lap.step) && all(is.finite(cov.lap.step))) {
cov.lap.step <- prep_Sig_mh(cov.lap.step)
Sig.mh <- prep_Sig_mh((1 - mh.laplace.refresh.weight) * Sig.mh + mh.laplace.refresh.weight * cov.lap.step)
laplace_refreshed <- TRUE
laplace_refresh_success <- laplace_refresh_success + 1L
}
}
if (!identical(mh.proposal, "slice") &&
(!isTRUE(sigmagam_ctrl$delay_adapt_until_after_warmup) || !isTRUE(sigmagam_warmup_active)) &&
mh.adapt && i >= mh.min_burn_adapt && i < n.burn && (i %% mh.adapt.interval == 0)) {
acc.win <- window.accept / pmax(window.total, 1L)
if (acc.win < mh.target.accept[1]) {
mh.scale <- mh.scale * (1 - mh.max_scale.step)
} else if (acc.win > mh.target.accept[2]) {
mh.scale <- mh.scale * (1 + mh.max_scale.step)
}
mh.scale <- min(max(mh.scale, mh.scale.bounds[1]), mh.scale.bounds[2])
Sig.scaled <- prep_Sig_mh(Sig.mh * (mh.scale^2))
chol_Sig.mh <- build_chol(Sig.scaled)
adapt.history <- rbind(
adapt.history,
data.frame(
iter = i,
window_accept = acc.win,
mh_scale = mh.scale,
sig11 = Sig.scaled[1, 1],
sig22 = Sig.scaled[2, 2],
laplace_refreshed = isTRUE(laplace_refreshed),
stringsAsFactors = FALSE
)
)
window.accept <- 0L
window.total <- 0L
}
if(!identical(mh.proposal, "slice") && i==n.burn && joint.sample && !isTRUE(sigmagam_warmup_active)){
Sig.mh = stats::cov(cbind(init.log.sigma[1:n.burn],init.logit.gamma[1:n.burn]))
Sig.mh <- prep_Sig_mh(Sig.mh)
chol_Sig.mh <- build_chol(Sig.mh * (mh.scale^2))
}
}else{
if (!identical(mh.proposal, "slice") && !isTRUE(sigmagam_warmup_active)) {
n.trial.keep <- n.trial.keep + 1L
n.accept.keep <- n.accept.keep + as.integer(isTRUE(lsiglgam.out$accept))
}
save.sigma[(i-n.burn)] = cursam.sigma
save.gamma[(i-n.burn)] = cursam.gamma
save.theta[,,(i-n.burn)] = cursam.theta
save.Ut[,(i-n.burn)] = cursam.Ut
save.st[,(i-n.burn)] = cursam.st
save.post.pred[,(i-n.burn)] = theta.out$post.pred
}
}
run.time = tictoc::toc(quiet = TRUE)
.exdqlm_progress(
"MCMC done",
model = "exDQLM",
status = "complete",
iter = I,
runtime_sec = run.time$toc - run.time$tic,
accept = if (identical(mh.proposal, "slice")) NULL else n.accept / I,
.verbose = verbose
)
safe_progress_callback(list(
event = "complete",
iter = as.integer(I),
total_iter = as.integer(I),
phase = "done",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = cursam.gamma,
kernel = mh.proposal,
accept = if (identical(mh.proposal, "slice")) NA_real_ else n.accept / I,
runtime_sec = as.numeric(run.time$toc - run.time$tic)
))
# exdqlm MAP standard forecast errors
map.gam = mean(save.gamma)
map.sig = mean(save.sigma)
map.st = rowMeans(save.st)
map.Ut = rowMeans(save.Ut)
tau = p.fn(p0,map.gam)
a_tau = (1-2*tau)/(tau*(1-tau))
b_tau = (2)/(tau*(1-tau))
c_tau = (as.numeric(map.gam>0)-tau)^(-1)
theta.out <- smoothed_theta(map.sig*c_tau*abs(map.gam)*map.st+map.Ut*a_tau,b_tau*map.Ut*map.sig)
map.standard.forecast.errors = theta.out$standard.forecast.errors
Sig.mh.final <- prep_Sig_mh(Sig.mh * (mh.scale^2))
ess_sigma <- tryCatch(as.numeric(coda::effectiveSize(coda::as.mcmc(save.sigma))), error = function(e) NA_real_)
ess_gamma <- tryCatch(as.numeric(coda::effectiveSize(coda::as.mcmc(save.gamma))), error = function(e) NA_real_)
chain_health_sigma <- .exdqlm_chain_health_metrics(save.sigma, n_keep = n.mcmc)
chain_health_gamma <- .exdqlm_chain_health_metrics(save.gamma, n_keep = n.mcmc)
accept_total <- if (identical(mh.proposal, "slice")) NA_real_ else n.accept / I
accept_burn <- if (identical(mh.proposal, "slice")) NA_real_ else if (n.trial.burn > 0) n.accept.burn / n.trial.burn else NA_real_
accept_keep <- if (identical(mh.proposal, "slice")) NA_real_ else if (n.trial.keep > 0) n.accept.keep / n.trial.keep else NA_real_
kernel_exact <- mh.proposal %in% c("rw", "laplace_rw", "slice")
mh.diag <- list(
proposal = mh.proposal,
adapt = if (identical(mh.proposal, "slice")) FALSE else mh.adapt,
joint_sample = isTRUE(joint.sample),
adapt_interval = if (identical(mh.proposal, "slice")) NA_integer_ else mh.adapt.interval,
target_accept = if (identical(mh.proposal, "slice")) c(NA_real_, NA_real_) else mh.target.accept,
scale_bounds = if (identical(mh.proposal, "slice")) c(NA_real_, NA_real_) else mh.scale.bounds,
scale_final = if (identical(mh.proposal, "slice")) NA_real_ else mh.scale,
joint_sigma_gamma = mh.proposal %in% c("rw", "laplace_rw"),
transformed_state = if (mh.proposal %in% c("rw", "laplace_rw")) c("log_sigma", "logit_gamma") else c("gamma"),
# Backward-compatible aliases used by some diagnostics scripts.
final_scale = if (identical(mh.proposal, "slice")) NA_real_ else mh.scale,
slice_width = if (identical(mh.proposal, "slice")) slice.width else NA_real_,
slice_max_steps = if (identical(mh.proposal, "slice")) slice.max.steps else NA_real_,
laplace_refresh = list(
enabled = identical(mh.proposal, "laplace_rw"),
interval = if (identical(mh.proposal, "laplace_rw")) as.integer(mh.laplace.refresh.interval) else NA_integer_,
start = if (identical(mh.proposal, "laplace_rw")) as.integer(mh.laplace.refresh.start) else NA_integer_,
weight = if (identical(mh.proposal, "laplace_rw")) as.numeric(mh.laplace.refresh.weight) else NA_real_,
attempts = if (identical(mh.proposal, "laplace_rw")) as.integer(laplace_refresh_attempts) else NA_integer_,
success = if (identical(mh.proposal, "laplace_rw")) as.integer(laplace_refresh_success) else NA_integer_
),
kernel_exact = kernel_exact,
signoff_ready = kernel_exact,
approximation_note = NA_character_,
accept = list(
total = accept_total,
burn = accept_burn,
kept = accept_keep,
n_accept = if (identical(mh.proposal, "slice")) NA_integer_ else n.accept,
n_total = if (identical(mh.proposal, "slice")) NA_integer_ else I
),
Sig.mh.initial = if (identical(mh.proposal, "slice")) matrix(NA_real_, 2, 2) else Sig.mh.initial,
Sig.mh.final = if (identical(mh.proposal, "slice")) matrix(NA_real_, 2, 2) else Sig.mh.final,
adaptation = if (identical(mh.proposal, "slice")) data.frame() else adapt.history,
adapt_trace = if (identical(mh.proposal, "slice")) data.frame() else adapt.history,
trace_enabled = trace.diagnostics,
trace_every = if (trace.diagnostics) trace.every else NA_integer_,
verbose_every = as.integer(verbose.every),
callback_every = as.integer(callback.every),
sigmagam = list(
freeze_burnin_iters = as.integer(sigmagam_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(sigmagam_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(sigmagam_ctrl$force_after_warmup),
delay_adapt_until_after_warmup = isTRUE(sigmagam_ctrl$delay_adapt_until_after_warmup),
delay_laplace_refresh_until_after_warmup = isTRUE(sigmagam_ctrl$delay_laplace_refresh_until_after_warmup),
first_active_iter = if (is.na(sigmagam_first_active_iter)) NA_integer_ else as.integer(sigmagam_first_active_iter),
updates_burn = as.integer(sigmagam_updates_burn),
updates_keep = as.integer(sigmagam_updates_keep),
update_count = as.integer(sigmagam_update_count),
postwarmup_update_count = as.integer(sigmagam_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(sigmagam_frozen_trace[seq_len(n.burn)]) else NA_real_
),
latent_state = list(
mode = latent_ctrl$mode,
freeze_burnin_iters = as.integer(latent_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(latent_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(latent_ctrl$force_after_warmup),
first_active_iter = if (is.na(latent_first_active_iter)) NA_integer_ else as.integer(latent_first_active_iter),
updates_burn = as.integer(latent_updates_burn),
updates_keep = as.integer(latent_updates_keep),
update_count = as.integer(latent_update_count),
postwarmup_update_count = as.integer(latent_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(latent_frozen_trace[seq_len(n.burn)]) else NA_real_
),
theta_state = list(
freeze_burnin_iters = as.integer(theta_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(theta_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(theta_ctrl$force_after_warmup),
first_active_iter = if (is.na(theta_first_active_iter)) NA_integer_ else as.integer(theta_first_active_iter),
updates_burn = as.integer(theta_updates_burn),
updates_keep = as.integer(theta_updates_keep),
update_count = as.integer(theta_update_count),
postwarmup_update_count = as.integer(theta_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(theta_frozen_trace[seq_len(n.burn)]) else NA_real_
),
trace = if (trace.diagnostics && trace_idx > 0L) {
do.call(rbind, trace_rows[seq_len(trace_idx)])
} else {
data.frame()
}
)
# exdqlm results
retlist = list(y=y,run.time=(run.time$toc-run.time$tic),model=model,p0=p0,df=df,dim.df=dim.df,dqlm.ind=dqlm.ind,
samp.theta = coda::as.mcmc(save.theta), theta.out = theta.out,
samp.post.pred = save.post.pred, map.standard.forecast.errors = map.standard.forecast.errors,
samp.sigma = coda::as.mcmc(save.sigma), samp.gamma = coda::as.mcmc(save.gamma),
init.log.sigma = coda::as.mcmc(init.log.sigma), init.logit.gamma = coda::as.mcmc(init.logit.gamma),
samp.vts = coda::as.mcmc(save.Ut), samp.sts = coda::as.mcmc(save.st),
accept.rate = accept_total,
accept.rate.burn = accept_burn,
accept.rate.keep = accept_keep,
Sig.mh = if (identical(mh.proposal, "slice")) matrix(NA_real_, 2, 2) else Sig.mh.final,
init.from.vb = init.from.vb,
vb.init.method = if (init.from.vb) vb.ctrl$method else NA_character_,
mh.diagnostics = mh.diag,
diagnostics = list(
mh = mh.diag,
progress = list(
verbose_every = as.integer(verbose.every),
callback_every = as.integer(callback.every)
),
ess = list(sigma = ess_sigma, gamma = ess_gamma),
chain_health = list(
sigma = chain_health_sigma,
gamma = chain_health_gamma
),
latent_state = mh.diag$latent_state,
theta_state = mh.diag$theta_state,
sigmagam = mh.diag$sigmagam,
s_block = list(
trace = mh.diag$trace,
final = if (is.data.frame(mh.diag$trace) && nrow(mh.diag$trace)) {
as.list(mh.diag$trace[nrow(mh.diag$trace), , drop = FALSE])
} else {
list()
}
),
rhat_ready = list(
sigma = as.numeric(save.sigma),
gamma = as.numeric(save.gamma)
),
sigmagam_trace = list(
frozen = sigmagam_frozen_trace,
update_reason = sigmagam_update_reason_trace,
forced_postwarmup = sigmagam_forced_postwarmup_trace,
update_performed = sigmagam_update_performed_trace,
update_count = sigmagam_update_count_trace
),
latent_state_trace = list(
frozen = latent_frozen_trace,
update_reason = latent_update_reason_trace,
forced_postwarmup = latent_forced_postwarmup_trace,
update_performed = latent_update_performed_trace,
update_count = latent_update_count_trace
),
theta_state_trace = list(
frozen = theta_frozen_trace,
update_reason = theta_update_reason_trace,
forced_postwarmup = theta_forced_postwarmup_trace,
update_performed = theta_update_performed_trace,
update_count = theta_update_count_trace
)
),
n.burn=n.burn,n.mcmc=n.mcmc)
}else{
######## DQLM
# fixed AL parameters
a_tau = (1-2*p0)/(p0*(1-p0))
b_tau = (2)/(p0*(1-p0))
# dqlm function sample theta ffbs
samp_theta<-function(ex.f,ex.q,sigma){
# initialize ffbs
m <- sam.theta <- matrix(NA,p,TT)
C <- array(NA,c(p,p,TT))
standard.forecast.errors <- post.pred <- rep(NA,TT)
## forward filter
# first iteration
a = as.vector(GG[,,1]%*%m0)
P = .exdqlm_regularize_cov(GG[,,1]%*%C0%*%t(GG[,,1]), context = "mcmc_al_sample_P_t1")
R = .exdqlm_regularize_cov(P + df.mat*P, context = "mcmc_al_sample_R_t1")
f = t(FF[,1])%*%a + ex.f[1]
q = .exdqlm_regularize_var(t(FF[,1])%*%R%*%FF[,1] + ex.q[1], context = "mcmc_al_sample_q_t1")
m[,1] = a + t(R)%*%FF[,1]%*%(y[1]-f)/q[1]
C[,,1] = .exdqlm_regularize_cov(
R - t(R)%*%FF[,1]%*%t(FF[,1])%*%R/q[1],
context = "mcmc_al_sample_C_t1"
)
standard.forecast.errors[1] = (y[1]-f)/sqrt(q)
# t = 2:TT
for(t in 2:TT){
a = as.vector(GG[,,t]%*%m[,(t-1)])
P = .exdqlm_regularize_cov(GG[,,t]%*%C[,,(t-1)]%*%t(GG[,,t]), context = sprintf("mcmc_al_sample_P_t%d", t))
R = .exdqlm_regularize_cov(P + df.mat*P, context = sprintf("mcmc_al_sample_R_t%d", t))
f = t(FF[,t])%*%a + ex.f[t]
fB = t(FF[,t])%*%R
q = .exdqlm_regularize_var(fB%*%FF[,t] + ex.q[t], context = sprintf("mcmc_al_sample_q_t%d", t))
m[,t] = a + t(fB)%*%(y[t]-f)/q[1]
C[,,t] = .exdqlm_regularize_cov(
R - t(fB)%*%fB/q[1],
context = sprintf("mcmc_al_sample_C_t%d", t)
)
standard.forecast.errors[t] = (y[t]-f)/sqrt(q)
}
## backwards sample
sC_TT = .exdqlm_regularize_cov(C[,,TT], context = "mcmc_al_sample_sC_TT")
svd.sC = svd(sC_TT)
sam.theta[,TT] = m[,TT] + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta <- numeric(TT)
reg_theta[TT] <- drop(crossprod(FF[,TT], sam.theta[,TT]))
post.pred[TT] = rexal(1,p0,reg_theta[TT],sigma,0)
for(t in (TT-1):1){
P = .exdqlm_regularize_cov(GG[,,(t+1)]%*%C[,,(t)]%*%t(GG[,,(t+1)]), context = sprintf("mcmc_al_back_P_t%d", t + 1L))
R.info = .exdqlm_cov_inverse(P + df.mat*P, context = sprintf("mcmc_al_back_R_t%d", t + 1L))
sB = C[,,t]%*%t(GG[,,(t+1)])%*%R.info$inverse
sm = m[,t] + sB%*%(sam.theta[,(t+1)]-as.vector(GG[,,(t+1)]%*%m[,(t)]))
sC = .exdqlm_regularize_cov(
C[,,t] - sB%*%GG[,,(t+1)]%*%C[,,t],
context = sprintf("mcmc_al_back_sC_t%d", t)
)
svd.sC = svd(sC)
sam.theta[,t] = sm + svd.sC$u%*%diag(sqrt(svd.sC$d),p)%*%stats::rnorm(p,0,1)
reg_theta[t] <- drop(crossprod(FF[,t], sam.theta[,t]))
post.pred[t] = rexal(1,p0,reg_theta[t],sigma,0)
}
return(list(standard.forecast.errors=standard.forecast.errors,post.pred=post.pred,sam.theta=sam.theta,fm=m,fC=C))
}
if (use_cpp_mcmc) {
samp_theta <- function(ex.f, ex.q, sigma) {
out <- cpp_ffbs_sample(ex.f, ex.q)
sam.theta <- out$sam.theta
reg_theta <- state_signal(FF, sam.theta)
post.pred <- vapply(seq_len(TT), function(t) {
rexal(1, p0, reg_theta[t], sigma, 0)
}, numeric(1))
list(
standard.forecast.errors = out$standard.forecast.errors,
post.pred = post.pred,
sam.theta = sam.theta,
fm = out$fm,
fC = out$fC
)
}
}
# dqlm function sample uts
samp_uts<-function(reg1,sigma){
chi <- as.numeric(((y-reg1)^2)/(b_tau*sigma))
psi <- (a_tau^2)/(b_tau*sigma) + (2/sigma)
sample_gig_cpp_required(chi = chi, psi = psi, lambda = 0.5, context = "dqlm_mcmc_uts")
}
# dqlm function sample sigma
samp_sigma<-function(reg1,uts){
1/stats::rgamma(1, shape = PriorSigma$a_sig + 1.5*TT,
rate = PriorSigma$b_sig + 0.5*sum( ((as.vector(y) - reg1 - a_tau*uts)^2)/(b_tau*uts) ) + sum(uts) )
}
callback.every <- if (trace.diagnostics) {
trace.every
} else {
100L
}
callback.every <- max(1L, as.integer(callback.every)[1])
trace_rows <- if (trace.diagnostics) vector("list", ceiling(I / trace.every)) else NULL
trace_idx <- 0L
latent_frozen_trace <- rep(FALSE, I)
latent_forced_postwarmup_trace <- rep(FALSE, I)
latent_update_performed_trace <- rep(FALSE, I)
latent_update_reason_trace <- rep(NA_character_, I)
latent_update_count_trace <- integer(I)
latent_first_active_iter <- NA_integer_
latent_update_count <- 0L
latent_postwarmup_update_count <- 0L
latent_updates_burn <- 0L
latent_updates_keep <- 0L
theta_frozen_trace <- rep(FALSE, I)
theta_forced_postwarmup_trace <- rep(FALSE, I)
theta_update_performed_trace <- rep(FALSE, I)
theta_update_reason_trace <- rep(NA_character_, I)
theta_update_count_trace <- integer(I)
theta_first_active_iter <- NA_integer_
theta_update_count <- 0L
theta_postwarmup_update_count <- 0L
theta_updates_burn <- 0L
theta_updates_keep <- 0L
dqlm_sigma_frozen_trace <- rep(FALSE, I)
dqlm_sigma_forced_postwarmup_trace <- rep(FALSE, I)
dqlm_sigma_update_performed_trace <- rep(FALSE, I)
dqlm_sigma_update_reason_trace <- rep(NA_character_, I)
dqlm_sigma_update_count_trace <- integer(I)
dqlm_sigma_first_active_iter <- NA_integer_
dqlm_sigma_update_count <- 0L
dqlm_sigma_postwarmup_update_count <- 0L
dqlm_sigma_updates_burn <- 0L
dqlm_sigma_updates_keep <- 0L
# Sample from dqlm posterior
tictoc::tic()
.exdqlm_progress(
"MCMC start",
burn = n.burn,
keep = n.mcmc,
kernel = "conjugate",
warm_start = if (isTRUE(init.from.vb)) vb.ctrl$method else "none",
.verbose = verbose
)
safe_progress_callback(list(
event = "start",
iter = 0L,
total_iter = as.integer(I),
phase = "burn",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = NA_real_,
kernel = "conjugate",
accept = NA_real_
))
for (i in 1:I){
current_iter <- as.integer(i)
# counter
if (i %% verbose.every == 0L) {
kept_now <- if (i <= n.burn) 0L else (i - n.burn)
.exdqlm_progress(
"MCMC progress",
model = "DQLM",
phase = if (i <= n.burn) "burn" else "keep",
iter = sprintf("%d/%d", i, I),
kept = sprintf("%d/%d", kept_now, n.mcmc),
.verbose = verbose
)
}
if (i %% callback.every == 0L) {
safe_progress_callback(list(
event = "progress",
iter = as.integer(i),
total_iter = as.integer(I),
phase = if (i <= n.burn) "burn" else "keep",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = NA_real_,
kernel = "conjugate",
accept = NA_real_
))
}
# sample theta
ex.f = cursam.Ut*a_tau
ex.q = b_tau*cursam.Ut*cursam.sigma
theta_state_sched <- warmup_state(i, theta_ctrl, theta_postwarmup_update_count)
if (isTRUE(theta_state_sched$active)) {
if (any(!is.finite(cursam.theta))) {
stop(sprintf(
"dqlm_mcmc_theta_warmup_state (iter=%d) invalid warmup theta state: nonfinite=%d",
i, sum(!is.finite(cursam.theta))
))
}
reg_theta_hold <- state_signal(FF, cursam.theta)
theta.out <- list(
standard.forecast.errors = rep(NA_real_, TT),
post.pred = vapply(seq_len(TT), function(t) {
rexal(1, p0, reg_theta_hold[t], cursam.sigma, 0)
}, numeric(1)),
sam.theta = cursam.theta,
fm = matrix(NA_real_, p, TT),
fC = array(NA_real_, c(p, p, TT))
)
} else {
theta.out <- samp_theta(ex.f,ex.q,cursam.sigma)
cursam.theta = theta.out$sam.theta
}
theta_frozen_trace[i] <- isTRUE(theta_state_sched$active)
theta_update_reason_trace[i] <- theta_state_sched$reason
if (!isTRUE(theta_state_sched$active)) {
theta_update_performed_trace[i] <- TRUE
theta_update_count <- theta_update_count + 1L
theta_update_count_trace[i] <- theta_update_count
if (is.na(theta_first_active_iter)) theta_first_active_iter <- as.integer(i)
if (i <= n.burn) {
theta_updates_burn <- theta_updates_burn + 1L
} else {
theta_updates_keep <- theta_updates_keep + 1L
}
if (theta_ctrl$freeze_burnin_iters > 0L) {
theta_postwarmup_update_count <- theta_postwarmup_update_count + 1L
}
theta_forced_postwarmup_trace[i] <- isTRUE(theta_state_sched$force_now)
} else {
theta_update_count_trace[i] <- theta_update_count
theta_forced_postwarmup_trace[i] <- FALSE
}
# sample uts
reg1 = state_signal(FF, cursam.theta)
latent_state_sched <- warmup_state(i, latent_ctrl, latent_postwarmup_update_count)
if (!is.finite(cursam.sigma) || cursam.sigma <= 0 || any(!is.finite(reg1))) {
stop(sprintf(
"dqlm_mcmc_pre_uts (iter=%d reason=%s) invalid state before chi update: sigma=%s reg1_finite=%s max_abs_reg1=%s max_abs_theta=%s",
i,
latent_state_sched$reason,
format(cursam.sigma, digits = 6),
all(is.finite(reg1)),
if (all(is.finite(reg1))) format(max(abs(reg1), na.rm = TRUE), digits = 6) else "NA",
if (all(is.finite(cursam.theta))) format(max(abs(cursam.theta), na.rm = TRUE), digits = 6) else "NA"
))
}
if (isTRUE(latent_state_sched$active)) {
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"dqlm_mcmc_latent_warmup_state (iter=%d) invalid warmup Ut state: nonfinite=%d nonpositive=%d",
i, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut)
))
}
} else {
cursam.Ut <- samp_uts(reg1,cursam.sigma)
if (count_nonfinite(cursam.Ut) > 0L || count_nonpositive(cursam.Ut) > 0L) {
stop(sprintf(
"dqlm_mcmc_post_uts (iter=%d) invalid Ut draws: nonfinite=%d nonpositive=%d sigma=%s max_abs_reg1=%s",
i, count_nonfinite(cursam.Ut), count_nonpositive(cursam.Ut),
format(cursam.sigma, digits = 6), max_abs_finite(reg1)
))
}
}
latent_frozen_trace[i] <- isTRUE(latent_state_sched$active)
latent_update_reason_trace[i] <- latent_state_sched$reason
if (!isTRUE(latent_state_sched$active)) {
latent_update_performed_trace[i] <- TRUE
latent_update_count <- latent_update_count + 1L
latent_update_count_trace[i] <- latent_update_count
if (is.na(latent_first_active_iter)) latent_first_active_iter <- as.integer(i)
if (i <= n.burn) {
latent_updates_burn <- latent_updates_burn + 1L
} else {
latent_updates_keep <- latent_updates_keep + 1L
}
if (latent_ctrl$freeze_burnin_iters > 0L) {
latent_postwarmup_update_count <- latent_postwarmup_update_count + 1L
}
latent_forced_postwarmup_trace[i] <- isTRUE(latent_state_sched$force_now)
} else {
latent_update_count_trace[i] <- latent_update_count
latent_forced_postwarmup_trace[i] <- FALSE
}
# sample sigma
dqlm_sigma_sched <- warmup_state(i, dqlm_sigma_ctrl, dqlm_sigma_postwarmup_update_count)
if(!fix.sigma){
if (!isTRUE(dqlm_sigma_sched$active)) {
cursam.sigma <- samp_sigma(reg1,cursam.Ut)
}
if (!is.finite(cursam.sigma) || cursam.sigma <= 0) {
stop(sprintf(
"dqlm_mcmc_post_sigma (iter=%d reason=%s) invalid sigma draw: sigma=%s max_abs_reg1=%s max_abs_ut=%s",
i, dqlm_sigma_sched$reason, format(cursam.sigma, digits = 6), max_abs_finite(reg1), max_abs_finite(cursam.Ut)
))
}
}
dqlm_sigma_frozen_trace[i] <- isTRUE(dqlm_sigma_sched$active)
dqlm_sigma_update_reason_trace[i] <- dqlm_sigma_sched$reason
if (!isTRUE(dqlm_sigma_sched$active)) {
dqlm_sigma_update_performed_trace[i] <- TRUE
dqlm_sigma_update_count <- dqlm_sigma_update_count + 1L
dqlm_sigma_update_count_trace[i] <- dqlm_sigma_update_count
if (is.na(dqlm_sigma_first_active_iter)) dqlm_sigma_first_active_iter <- as.integer(i)
if (i <= n.burn) {
dqlm_sigma_updates_burn <- dqlm_sigma_updates_burn + 1L
} else {
dqlm_sigma_updates_keep <- dqlm_sigma_updates_keep + 1L
}
if (dqlm_sigma_ctrl$freeze_burnin_iters > 0L) {
dqlm_sigma_postwarmup_update_count <- dqlm_sigma_postwarmup_update_count + 1L
}
dqlm_sigma_forced_postwarmup_trace[i] <- isTRUE(dqlm_sigma_sched$force_now)
} else {
dqlm_sigma_update_count_trace[i] <- dqlm_sigma_update_count
dqlm_sigma_forced_postwarmup_trace[i] <- FALSE
}
if (trace.diagnostics && (i %% trace.every == 0L)) {
u_stats <- .exdqlm_trace_summary(cursam.Ut)
trace_idx <- trace_idx + 1L
trace_rows[[trace_idx]] <- data.frame(
iter = i,
phase = if (i <= n.burn) "burn" else "keep",
sigma = cursam.sigma,
latent_frozen = isTRUE(latent_state_sched$active),
latent_mode = latent_ctrl$mode,
latent_update_reason = latent_state_sched$reason,
latent_forced_postwarmup = isTRUE(latent_forced_postwarmup_trace[i]),
latent_update_performed = isTRUE(latent_update_performed_trace[i]),
latent_update_count = as.integer(latent_update_count_trace[i]),
theta_frozen = isTRUE(theta_state_sched$active),
theta_update_reason = theta_state_sched$reason,
theta_forced_postwarmup = isTRUE(theta_forced_postwarmup_trace[i]),
theta_update_performed = isTRUE(theta_update_performed_trace[i]),
theta_update_count = as.integer(theta_update_count_trace[i]),
dqlm_sigma_frozen = isTRUE(dqlm_sigma_sched$active),
dqlm_sigma_update_reason = dqlm_sigma_sched$reason,
dqlm_sigma_forced_postwarmup = isTRUE(dqlm_sigma_forced_postwarmup_trace[i]),
dqlm_sigma_update_performed = isTRUE(dqlm_sigma_update_performed_trace[i]),
dqlm_sigma_update_count = as.integer(dqlm_sigma_update_count_trace[i]),
u_mean = u_stats[["mean"]],
u_sd = u_stats[["sd"]],
u_q05 = u_stats[["q05"]],
u_q50 = u_stats[["median"]],
u_q95 = u_stats[["q95"]],
u_min = u_stats[["min"]],
u_max = u_stats[["max"]],
stringsAsFactors = FALSE
)
}
# save samples after burn
if(i > n.burn){
save.sigma[(i-n.burn)] = cursam.sigma
save.theta[,,(i-n.burn)] = cursam.theta
save.Ut[,(i-n.burn)] = cursam.Ut
save.post.pred[,(i-n.burn)] = theta.out$post.pred
}
}
run.time = tictoc::toc(quiet = TRUE)
.exdqlm_progress(
"MCMC done",
model = "DQLM",
status = "complete",
iter = I,
runtime_sec = run.time$toc - run.time$tic,
.verbose = verbose
)
safe_progress_callback(list(
event = "complete",
iter = as.integer(I),
total_iter = as.integer(I),
phase = "done",
n_burn = as.integer(n.burn),
n_mcmc = as.integer(n.mcmc),
sigma = cursam.sigma,
gamma = NA_real_,
kernel = "conjugate",
accept = NA_real_,
runtime_sec = as.numeric(run.time$toc - run.time$tic)
))
# dqlm MAP standard forecast errors
map.sig = mean(save.sigma)
map.Ut = rowMeans(save.Ut)
theta.out <- smoothed_theta(map.Ut*a_tau,b_tau*map.Ut*map.sig)
map.standard.forecast.errors = theta.out$standard.forecast.errors
ess_sigma <- tryCatch(as.numeric(coda::effectiveSize(coda::as.mcmc(save.sigma))), error = function(e) NA_real_)
# dqlm results
retlist = list(y=y,run.time=(run.time$toc-run.time$tic),model=model,p0=p0,df=df,dim.df=dim.df,dqlm.ind=dqlm.ind,
samp.theta = coda::as.mcmc(save.theta), theta.out = theta.out,
samp.post.pred = save.post.pred, map.standard.forecast.errors = map.standard.forecast.errors,
samp.sigma = coda::as.mcmc(save.sigma),
samp.vts = coda::as.mcmc(save.Ut),
init.from.vb = init.from.vb,
vb.init.method = if (init.from.vb) vb.ctrl$method else NA_character_,
diagnostics = list(
progress = list(
verbose_every = as.integer(verbose.every),
callback_every = as.integer(callback.every)
),
latent_state = list(
mode = latent_ctrl$mode,
freeze_burnin_iters = as.integer(latent_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(latent_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(latent_ctrl$force_after_warmup),
first_active_iter = if (is.na(latent_first_active_iter)) NA_integer_ else as.integer(latent_first_active_iter),
updates_burn = as.integer(latent_updates_burn),
updates_keep = as.integer(latent_updates_keep),
update_count = as.integer(latent_update_count),
postwarmup_update_count = as.integer(latent_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(latent_frozen_trace[seq_len(n.burn)]) else NA_real_
),
theta_state = list(
freeze_burnin_iters = as.integer(theta_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(theta_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(theta_ctrl$force_after_warmup),
first_active_iter = if (is.na(theta_first_active_iter)) NA_integer_ else as.integer(theta_first_active_iter),
updates_burn = as.integer(theta_updates_burn),
updates_keep = as.integer(theta_updates_keep),
update_count = as.integer(theta_update_count),
postwarmup_update_count = as.integer(theta_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(theta_frozen_trace[seq_len(n.burn)]) else NA_real_
),
dqlm_sigma = list(
freeze_burnin_iters = as.integer(dqlm_sigma_ctrl$freeze_burnin_iters),
freeze_only_during_burn = isTRUE(dqlm_sigma_ctrl$freeze_only_during_burn),
force_after_warmup = isTRUE(dqlm_sigma_ctrl$force_after_warmup),
first_active_iter = if (is.na(dqlm_sigma_first_active_iter)) NA_integer_ else as.integer(dqlm_sigma_first_active_iter),
updates_burn = as.integer(dqlm_sigma_updates_burn),
updates_keep = as.integer(dqlm_sigma_updates_keep),
update_count = as.integer(dqlm_sigma_update_count),
postwarmup_update_count = as.integer(dqlm_sigma_postwarmup_update_count),
frozen_burn_rate = if (n.burn > 0L) mean(dqlm_sigma_frozen_trace[seq_len(n.burn)]) else NA_real_
),
latent_state_trace = list(
frozen = latent_frozen_trace,
update_reason = latent_update_reason_trace,
forced_postwarmup = latent_forced_postwarmup_trace,
update_performed = latent_update_performed_trace,
update_count = latent_update_count_trace
),
theta_state_trace = list(
frozen = theta_frozen_trace,
update_reason = theta_update_reason_trace,
forced_postwarmup = theta_forced_postwarmup_trace,
update_performed = theta_update_performed_trace,
update_count = theta_update_count_trace
),
dqlm_sigma_trace = list(
frozen = dqlm_sigma_frozen_trace,
update_reason = dqlm_sigma_update_reason_trace,
forced_postwarmup = dqlm_sigma_forced_postwarmup_trace,
update_performed = dqlm_sigma_update_performed_trace,
update_count = dqlm_sigma_update_count_trace
),
trace = if (trace.diagnostics && trace_idx > 0L) {
do.call(rbind, trace_rows[seq_len(trace_idx)])
} else {
data.frame()
},
ess = list(sigma = ess_sigma, gamma = NA_real_),
rhat_ready = list(
sigma = as.numeric(save.sigma),
gamma = numeric(0)
)
),
n.burn=n.burn,n.mcmc=n.mcmc)
}
retlist$backend <- list(mcmc = mcmc_backend, mode = cpp_mcmc_mode, gig = gig_backend)
# return results
class(retlist) <- .exdqlm_fit_class("exdqlmMCMC")
return(retlist)
}
# Internal helper for diagnostics-only multichain validation.
.exdqlm_mcmc_multichain_diag <- function(n.chains = 4L, seeds = NULL, mcmc_args = list()) {
n.chains <- suppressWarnings(as.integer(n.chains)[1])
if (!is.finite(n.chains) || n.chains < 2L) {
stop("n.chains must be >= 2 for multichain diagnostics.")
}
if (is.null(seeds)) {
seeds <- seq_len(n.chains) + 20260300L
}
seeds <- as.integer(seeds)
if (length(seeds) != n.chains) {
stop("Length of seeds must match n.chains.")
}
fits <- vector("list", n.chains)
for (i in seq_len(n.chains)) {
set.seed(seeds[i])
args_i <- utils::modifyList(mcmc_args, list(verbose = FALSE))
fits[[i]] <- do.call(exdqlmMCMC, args_i)
}
sigma_list <- coda::mcmc.list(lapply(fits, function(f) coda::as.mcmc(as.numeric(f$samp.sigma))))
sigma_rhat <- tryCatch(
as.numeric(coda::gelman.diag(sigma_list, autoburnin = FALSE)$psrf[1, "Point est."]),
error = function(e) NA_real_
)
sigma_ess <- tryCatch(as.numeric(coda::effectiveSize(sigma_list))[1], error = function(e) NA_real_)
has_gamma <- !isTRUE(fits[[1]]$dqlm.ind) && !is.null(fits[[1]]$samp.gamma)
if (has_gamma) {
gamma_list <- coda::mcmc.list(lapply(fits, function(f) coda::as.mcmc(as.numeric(f$samp.gamma))))
gamma_rhat <- tryCatch(
as.numeric(coda::gelman.diag(gamma_list, autoburnin = FALSE)$psrf[1, "Point est."]),
error = function(e) NA_real_
)
gamma_ess <- tryCatch(as.numeric(coda::effectiveSize(gamma_list))[1], error = function(e) NA_real_)
} else {
gamma_list <- NULL
gamma_rhat <- NA_real_
gamma_ess <- NA_real_
}
list(
fits = fits,
seeds = seeds,
diagnostics = list(
sigma = list(rhat = sigma_rhat, ess = sigma_ess, chains = sigma_list),
gamma = list(rhat = gamma_rhat, ess = gamma_ess, chains = gamma_list)
)
)
}
Try the exdqlm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.