Learn R Programming

exdqlm (version 0.4.0)

exalStaticMCMC: exAL (static) - MCMC algorithm

Description

The function applies a Gibbs sampler for static Extended Asymmetric Laplace regression (exAL). We update \(\beta, v, s\) from their full conditionals, then update \((\sigma,\gamma)\) either jointly on transformed coordinates \((\ell,\eta)=(\log \sigma,\mathrm{logit}((\gamma-L)/(U-L)))\) (for "rw" and "laplace_rw") or with univariate gamma kernels ("slice", "slice_eta", "laplace_local"). Optional multi-refresh and global-jump controls are available to improve exAL mixing in hard cases.

Usage

exalStaticMCMC(
  y,
  X,
  p0,
  b0 = NULL,
  V0 = NULL,
  beta_prior = c("ridge", "rhs", "rhs_ns"),
  beta_prior_controls = NULL,
  a_sigma = 1,
  b_sigma = 1,
  gamma_bounds = c(L.fn(p0), U.fn(p0)),
  log_prior_gamma = NULL,
  init = NULL,
  dqlm.ind = FALSE,
  al.ind = NULL,
  n.burn = 2000,
  n.mcmc = 1500,
  thin = 1,
  init.from.vb = FALSE,
  vb_init_controls = NULL,
  vb_init_fit = NULL,
  mcmc_control = NULL,
  sigmagam_controls = NULL,
  mh.proposal = c("slice", "laplace_rw", "rw", "slice_eta", "laplace_local"),
  mh.adapt = TRUE,
  mh.adapt.interval = 50L,
  mh.target.accept = c(0.2, 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,
  gamma.substeps = 1L,
  p.global.eta.jump = 0,
  global.eta.jump.scale = 1,
  trace.diagnostics = TRUE,
  trace.every = 1L,
  verbose = TRUE,
  progress_callback = NULL
)

Value

An object of class "exalStaticMCMC" containing:

  • run.time - total wall time in seconds.

  • X, p0, bounds - design, quantile, and (L, U).

  • samp.beta - posterior sample of beta as coda::mcmc (n.mcmc x p).

  • samp.sigma - posterior sample of sigma as coda::mcmc.

  • samp.gamma - posterior sample of gamma as coda::mcmc.

  • samp.v - latent v draws as coda::mcmc (n.mcmc x n).

  • samp.s - latent s draws as coda::mcmc (n.mcmc x n).

  • samp.lambda, samp.tau, samp.c2 - RHS latent draws when an RHS-family prior is used; otherwise NULL.

  • beta_prior - prior metadata and, for RHS-family priors, posterior summaries of the shrinkage block. For "rhs_ns", the state tracks lambda2, nu, tau2, xi, and zeta2.

  • mh.diagnostics - proposal kernel diagnostics for the exAL gamma update, including whether the saved kernel is exact/signoff-ready.

  • rhs.diagnostics - RHS latent summaries and optional trace metadata when an RHS-family prior is used, including the resolved preflight configuration used at fit start.

  • last - last state of the chain (useful for restarts).

Arguments

y

Numeric vector of length \(n\).

X

Numeric matrix \(n \times p\) (design).

p0

Quantile level in \((0,1)\).

b0, V0

Prior mean and covariance for \(\beta\) (Normal). Defaults: \(b_0=\mathbf{0}_p\), \(V_0=10^6 I_p\).

beta_prior

Coefficient prior type: "ridge" (default), "rhs" (regularized horseshoe), or "rhs_ns" (Nishimura-Suchard regularized horseshoe with a closed-form inverse-gamma hierarchy for static inference).

beta_prior_controls

Optional list of prior-specific controls. For RHS-family priors (that is, when beta_prior is "rhs" or "rhs_ns"), supported keys include: tau0, nu, s or s2, shrink_intercept, intercept_prec, n_inner, eta_bounds, var_floor, h_curv, verbose, init_lambda, init_log_lambda, init_tau, init_log_tau, init_c2, init_log_c2, collapse_tau_ratio_tol, collapse_beta_max_abs_tol, collapse_invV_med_tol, collapse_beta_l2_tol, collapse_small_beta_frac_tol, small_beta_abs_tol, slice_width, and slice_max_steps. For beta_prior = "rhs_ns", optional slab controls a_zeta, b_zeta, and zeta2_fixed are also supported. In this mode, the local-global-slab block is represented by \((\lambda_j^2,\nu_j,\tau^2,\xi,\zeta^2)\) and updated with closed-form Gibbs steps. When beta_prior is "rhs" or "rhs_ns", b0 and V0 are ignored for the shrunk coefficients and retained only for backward-compatible ridge behavior. If both init_log_tau and init_tau are omitted (or NULL), the RHS global scale initializes at tau = 1 (init_log_tau = 0) instead of tau0. By default (shrink_intercept = FALSE), the intercept is excluded from horseshoe shrinkage and uses intercept_prec.

a_sigma, b_sigma

Hyperparameters for an inverse-gamma prior on sigma, with density proportional to sigma^{-(a_sigma+1)} exp(-b_sigma/sigma).

gamma_bounds

Numeric length-2 vector (L, U) for gamma. Defaults to c(L.fn(p0), U.fn(p0)).

log_prior_gamma

Function function(g) log pi(g) for gamma on (L, U). Default is a truncated Student-t prior centered at 0 on the admissible support.

init

Optional list with starting values: beta, sigma, gamma, v (length \(n\)), s (length \(n\)), and for RHS-family priors optionally lambda, tau, and c2. For beta_prior = "rhs_ns", the squared-scale parameterization lambda2, tau2, zeta2, and optional auxiliaries nu, xi are also accepted. Missing pieces are filled sensibly.

dqlm.ind

Logical; if TRUE, fit the reduced AL model (gamma = 0), corresponding to Bayesian linear quantile regression under the AL working likelihood. This removes the gamma- and s-blocks and leaves conjugate Gibbs updates for beta, sigma, and v. This argument is retained for consistency with the dynamic exDQLM API; in static examples, al.ind = TRUE is the clearer spelling for this AL special case.

al.ind

Optional static-model alias for dqlm.ind. Prefer this argument when requesting the static AL special case. When supplied, this flag maps directly to dqlm.ind. If both are given, they must agree.

n.burn

Number of burn-in iterations. Default 2000.

n.mcmc

Number of kept MCMC iterations (after burn). Default 1500.

thin

Integer; save every thin-th iteration after burn. We internally run n.burn + n.mcmc * thin iterations to return exactly n.mcmc saved draws.

init.from.vb

Logical; if TRUE, run static VB first and use its posterior moments as MCMC initialization.

vb_init_controls

Optional list controlling VB warm start. Supported keys: max_iter, tol, n_samp_xi, verbose, and ld_controls (passed through to exalStaticLDVB()).

vb_init_fit

Optional precomputed static VB fit object used as the warm-start reference when init.from.vb = TRUE.

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, the package applies its conservative default exAL (sigma, gamma) warmup profile and keeps the RHS-family tau warmup defaults active through the shared prior layer.

sigmagam_controls

Optional list controlling warmup/freeze behavior for the nonconjugate (sigma, gamma) block. Prefer exal_make_mcmc_sigmagam_control() or the sigmagam block of exal_make_mcmc_control() for new code; this argument remains available as a direct advanced override.

mh.proposal

Character string controlling the exAL nonconjugate update kernel. "slice" (default) uses an exact bounded univariate slice sampler on gamma (with sigma updated from its conditional), and "slice_eta" does the same on transformed eta. "laplace_rw" uses a Laplace-informed adaptive random-walk MH update on the transformed joint block \((\eta,\ell)=(\mathrm{logit}((\gamma-L)/(U-L)), \log\sigma)\). "rw" uses the same exact joint update with identity base covariance. "laplace_local" reproduces the prior approximate local-Gaussian gamma draw retained for compatibility and not recommended for routine use. Only "laplace_local" is approximate.

mh.adapt

Logical; adapt the random-walk proposal scale during burn-in for "rw" and "laplace_rw". Ignored for "laplace_local", "slice", and "slice_eta".

mh.adapt.interval

Integer adaptation window for RW-based kernels.

mh.target.accept

Numeric length-2 target acceptance band.

mh.scale.bounds

Numeric length-2 lower/upper bounds for RW proposal scale multiplier.

mh.max_scale.step

Numeric multiplicative adaptation cap in (0,1).

mh.min_burn_adapt

Integer minimum burn-in before adaptation starts.

slice.width

Positive numeric width for bounded slice updates when mh.proposal = "slice" or "slice_eta".

slice.max.steps

Positive integer or Inf; maximum stepping-out expansions for the slice sampler.

gamma.substeps

Positive integer; number of consecutive gamma-kernel refreshes per outer MCMC iteration. Default 1.

p.global.eta.jump

Numeric in [0,1]; per-substep probability of proposing a global independence-MH move on eta (the logit transform of gamma) using a Laplace proposal with MH correction. Default 0.

global.eta.jump.scale

Positive numeric scale multiplier applied to the Laplace proposal SD used in global eta jumps.

trace.diagnostics

Logical; if TRUE, retain per-iteration gamma/s diagnostics under mh.diagnostics$trace. Set FALSE for lighter-weight runs.

trace.every

Positive integer; when trace.diagnostics=TRUE, record one diagnostics row every trace.every iterations.

verbose

Print progress every progress_every iterations.

progress_callback

Optional callback invoked with a named list at MCMC start, at each progress checkpoint, and on completion. Intended for workflow-level per-case progress logging.

Examples

Run this code
# \donttest{
set.seed(123)
n <- 60; p <- 3
X <- cbind(1, rnorm(n), rnorm(n))
beta0 <- c(0.5, -1, 0.8); sigma0 <- 1.2
y <- as.numeric(X %*% beta0 + rnorm(n, 0, sigma0))
fit <- exalStaticMCMC(
  y, X, p0 = 0.5, n.burn = 60, n.mcmc = 60, thin = 1, verbose = FALSE
)
summary(fit$samp.beta)

fit_rhs <- exalStaticMCMC(
  y, X, p0 = 0.5,
  beta_prior = "rhs",
  beta_prior_controls = list(tau0 = 0.5, nu = 3, s2 = 1, shrink_intercept = FALSE),
  n.burn = 50, n.mcmc = 50, thin = 1, mh.proposal = "slice", verbose = FALSE
)
fit_rhs$beta_prior$type

fit_rhs_ns <- exalStaticMCMC(
  y, X, p0 = 0.5,
  beta_prior = "rhs_ns",
  beta_prior_controls = list(tau0 = 0.5, a_zeta = 1.5, b_zeta = 1, zeta2_fixed = 1),
  n.burn = 40, n.mcmc = 40, thin = 1, mh.proposal = "slice", verbose = FALSE
)
fit_rhs_ns$beta_prior$type

fit_al <- exalStaticMCMC(
  y, X, p0 = 0.5,
  al.ind = TRUE,
  n.burn = 50, n.mcmc = 50, thin = 1, verbose = FALSE
)
fit_al$dqlm.ind
# }

Run the code above in your browser using DataLab