EnvelopeDispersionBuild: Builds Dispersion-Aware Envelope for Simulation
In glmbayes: Bayesian Generalized Linear Models (IID Samples)

EnvelopeDispersionBuild

R Documentation

Builds Dispersion-Aware Envelope for Simulation

Description

Constructs a dispersion-aware envelope for simulation in Gaussian models with uncertain variance. This function extrapolates the coefficient envelope across a high-probability interval for the dispersion parameter sigma^2, and builds a global upper bound for the log-posterior remainder. It also computes mixture weights for envelope faces and adjusts the Gamma proposal for precision.

The envelope is constructed using the slopes of the face constants with respect to dispersion, evaluated at an anchor point. The resulting structure supports exact i.i.d. sampling via accept-reject correction.

The procedure follows these steps:

Posterior precision (Gamma) parameters. Using the prior and posterior-predictive RSS, set

\mathrm{shape2} = \mathrm{Shape} + n_{\mathrm{obs}}/2, \quad \mathrm{rate3} = \mathrm{Rate} + \mathrm{RSS}_{\mathrm{post}}/2

These parameterize the posterior precision v = 1/\sigma^2 \sim \mathrm{Gamma}(\mathrm{shape2}, \mathrm{rate3}).
Central credible interval for dispersion (low, upp). Choose a central mass level max_disp_perc (e.g., 0.99) for precision, then invert the corresponding Gamma quantiles to dispersion:

\mathrm{low} = 1 / Q_{\Gamma}(\mathrm{max\_disp\_perc}; \mathrm{shape2}, \mathrm{rate3}), \quad \mathrm{upp} = 1 / Q_{\Gamma}(1 - \mathrm{max\_disp\_perc}; \mathrm{shape2}, \mathrm{rate3})

The interval [\mathrm{low}, \mathrm{upp}] is the domain over which all envelopes must dominate.
Face slopes at an anchor (dispstar).

\mathrm{dispstar} = \mathrm{rate3} / (\mathrm{shape2} - 1)

(posterior mean of \sigma^2). Compute \mathrm{New\_LL\_Slope}_j for each face j.
Linear extrapolation of face constants.

\theta^{\mathrm{low}}_{\bar{j}} = \theta^{\mathrm{base}}_{\bar{j}} + (\mathrm{low} - \mathrm{dispstar}) \cdot \mathrm{New\_LL\_Slope}_j

\theta^{\mathrm{upp}}_{\bar{j}} = \theta^{\mathrm{base}}_{\bar{j}} + (\mathrm{upp} - \mathrm{dispstar}) \cdot \mathrm{New\_LL\_Slope}_j
Global upper line and endpoint maxima.

\mathrm{max\_low} = \max_j \theta^{\mathrm{low}}_{\bar{j}}, \quad \mathrm{max\_upp} = \max_j \theta^{\mathrm{upp}}_{\bar{j}}

\mathrm{new\_slope} = (\mathrm{max\_upp} - \mathrm{max\_low}) / (\mathrm{upp} - \mathrm{low}), \quad \mathrm{new\_int} = \mathrm{max\_low} - \mathrm{new\_slope} \cdot \mathrm{low}
Face slack and mixture weights.

\mathrm{lg\_prob\_factor}_j = \max\big(\theta^{\mathrm{upp}}_{\bar{j}} - \mathrm{max\_upp},\; \theta^{\mathrm{low}}_{\bar{j}} - \mathrm{max\_low}\big)

Combine with \mathrm{New\_logP2}_j = \mathrm{logP}_j + \tfrac{1}{2}\|\bar{c}_j\|^2 to form mixture weights \mathrm{PLSD}_j \propto \exp\!\bigl(\mathrm{New\_logP2}_j + \mathrm{lg\_prob\_factor}_j\bigr).
Gamma tilt and dispersion-axis envelope.

\mathrm{dispstar} = (\mathrm{upp} - \mathrm{low}) / \log(\mathrm{upp}/\mathrm{low})

\mathrm{lm\_log2} = \mathrm{new\_slope} \cdot \mathrm{dispstar}, \quad \mathrm{lm\_log1} = \mathrm{new\_int} + \mathrm{new\_slope} \cdot \mathrm{dispstar} - \mathrm{new\_slope} \cdot \log(\mathrm{dispstar})

Tilt the Gamma proposal via \mathrm{shape3} = \mathrm{shape2} - \mathrm{lm\_log2}.

Usage

EnvelopeDispersionBuild(
  Env, Shape, Rate, P, y, x, alpha, n_obs, RSS_post, RSS_ML,
  mu, wt, max_disp_perc = 0.99,
  disp_lower = NULL, disp_upper = NULL,
  verbose = FALSE, use_parallel = TRUE
)

Arguments

`Env`	Envelope object from `EnvelopeBuild`, containing tangency points and gradients
`Shape`	Prior shape parameter for precision `v = 1 / sigma^2`
`Rate`	Prior rate parameter for precision
`P`	Prior precision matrix for coefficients
`y`	Numeric response vector of length `m`
`x`	a design matrix of dimension `m * p`
`alpha`	Numeric offset vector of length `m`
`n_obs`	Number of observations
`RSS_post`	Expected posterior weighted residual sum of squares (i.e., `\mathbb{E}[\mathrm{RSS}(\beta)\mid y,\phi]` under the Normal posterior for `\beta` at fixed dispersion `\phi = d`). This value is used for dispersion anchoring / Gamma updates.
`RSS_ML`	Residual sum of squares associated with MLE estimate
`mu`	Prior mean parameter
`wt`	weight vector
`max_disp_perc`	Truncation level for dispersion (default 0.99)
`disp_lower`	lower bound truncation for dispersion
`disp_upper`	upper bound truncation for dispersion
`verbose`	Option to have verbose output
`use_parallel`	Logical. Whether to use parallel processing.

Details

This function is designed to complement EnvelopeBuild for Gaussian models with Normal-Gamma priors. It enables exact sampling of both coefficients and dispersion by constructing a joint envelope that respects posterior curvature in both dimensions.

The dispersion anchor point is chosen as the log-scale center of the credible interval, and the Gamma proposal is tilted to match the envelope slope at this point. Theory and narrative: \insertCiteNygren2006glmbayes; vignettes Chapter-A07, Chapter-A11; \insertCiteglmbayesChapterA08,glmbayesIndNormGammaVignetteglmbayes.

Value

EnvelopeDispersionBuild()

A list containing:

Env_out

Envelope object with updated mixture weights (PLSD)

gamma_list

Posterior Gamma tilt parameters

shape3: Adjusted shape parameter after slope correction
rate2: Posterior rate parameter, defined as Rate + rss_min_global/2
disp_upper: Upper bound of the dispersion interval \sigma^2
disp_lower: Lower bound of the dispersion interval \sigma^2

UB_list

Upper-bound diagnostics

RSS_ML: Residual sum of squares at the maximum-likelihood estimate
RSS_Min: Minimum residual sum of squares across envelope faces
max_New_LL_UB: Maximum extrapolated face constant at the upper dispersion bound
max_LL_log_disp: Log-posterior upper bound evaluated at disp_upper
lm_log1: Intercept term of the global upper line approximation
lm_log2: Slope term of the global upper line approximation
lg_prob_factor: Per-face slack factors used in mixture weighting
lmc1: Linear extrapolation constant (intercept)
lmc2: Linear extrapolation constant (slope)
UB2min: Minimum UB2 value across faces, used for diagnostics

diagnostics

Internal diagnostic values

dispstar: Anchor dispersion value (posterior mean or geometric mean)
New_LL_Slope: Vector of slopes of face constants at dispstar
shape2: Posterior shape parameter before tilt correction
rate3: Posterior rate parameter before tilt correction
shape3: Adjusted shape parameter (same as in gamma_list)
max_low: Maximum extrapolated face constant at the lower dispersion bound
max_upp: Maximum extrapolated face constant at the upper dispersion bound
new_slope: Slope of the global upper line across dispersion bounds
new_int: Intercept of the global upper line across dispersion bounds
prob_factor: Normalized mixture weights across faces
UB2min: Minimum UB2 diagnostic value (duplicated for consistency)

EnvBuildLinBound()

Numeric vector of slopes of face constants with respect to dispersion, evaluated at the anchor dispstar

thetabar_const()

Numeric vector of base face constants computed from tangency points and gradient vectors under prior precision P

Inv_f3_with_disp()

Numeric matrix of inverse function evaluations at a given dispersion and face subset, returned by the C++ routine _glmbayes_Inv_f3_with_disp

UB2()

Numeric scalar representing the UB2 upper-bound criterion for a given dispersion and face, defined as (1/dispersion) * (RSS - rss\_min\_global)

rss_face_at_disp()

Numeric scalar giving the residual sum of squares for a specified face at a given dispersion, computed from cached matrices and the inverse function evaluation

Use in accept/reject procedure

The accept/reject sampler relies on a decomposition of the log-posterior into a test statistic and several bounding terms. Each component is constructed so that its sign is controlled, ensuring the validity of the accept/reject step.

test1 (log-likelihood bound): Placeholder: explain how test1 is formed and why it is non-positive.
UB1 (if applicable): Placeholder: describe UB1's role and why it is non-negative.
UB2 (residual sum of squares bound): Placeholder: explain how UB2 is constructed from RSS differences and why it is non-negative.
UB3A (face-wise quadratic/linear envelope surplus): Placeholder: explain how lg_prob_factor, lmc1, and lmc2 are derived and why UB3A >= 0.
UB3B (dispersion-axis envelope surplus): Placeholder: explain how lm_log1, lm_log2, and max_New_LL_UB are used and why UB3B >= 0.

Together, these components define

test = test1 - UB2 - UB3A - UB3B,

with test1 \le 0 and each UB term \ge 0, ensuring the accept/reject procedure is valid and unbiased.

References

\insertAllCited

Examples

############################### Start of EnvelopeDispersionBuild example ####################

# This example mirrors the current C++ algorithm path for Gaussian regression
# with an independent Normal-Gamma prior:
#   rIndepNormalGammaReg:
#     - Step A: EnvelopeCentering (initial dispersion + dispersion anchoring loop)
#     - Step B: optimize posterior mode for coefficients (optim + f2/f3)
#     - Step C: standardize the model (glmb_Standardize_Model)
#     - Step D: build coefficient envelope (EnvelopeBuild)
#     - Step E: build dispersion-aware envelope (EnvelopeDispersionBuild)
#     - Step F: sort envelope components (EnvelopeSort)
# It stops after envelope construction (no standardized-envelope sampling).

ctl <- c(4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14)
trt <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69)
group <- gl(2, 10, 20, labels = c("Ctl", "Trt"))
weight <- c(ctl, trt)

ps <- Prior_Setup(weight ~ group, gaussian())

x <- as.matrix(ps$x)
y <- as.vector(ps$y)
mu <- ps$mu
Sigma <- ps$Sigma
shape <- ps$shape
rate <- ps$rate

n_obs <- length(y)
wt <- rep(1, n_obs)
offset2 <- rep(0, n_obs)

# Reconstruct coefficient precision P (matches rindepNormalGamma_reg)
Rchol <- chol(Sigma)
Pinv <- chol2inv(Rchol)
P <- 0.5 * (Pinv + t(Pinv))

famfunc <- glmbfamfunc(gaussian())
f2 <- famfunc$f2
f3 <- famfunc$f3

Gridtype_core <- as.integer(2)

###############################################################################
# Step A: EnvelopeCentering (initial dispersion + dispersion anchoring loop)
###############################################################################
centering <- EnvelopeCentering(
  y = y,
  x = x,
  mu = as.vector(mu),
  P = P,
  offset = offset2,
  wt = wt,
  shape = shape,
  rate = rate,
  Gridtype = Gridtype_core,
  verbose = FALSE
)


dispersion2 <- centering$dispersion
RSS_Post2 <- centering$RSS_post

n_w <- sum(wt)

###############################################################################
# Step B: Coefficient posterior mode optimization (optim + f2/f3)
###############################################################################
dispstar <- dispersion2

wt2_opt <- wt / dispstar
alpha <- as.vector(x %*% as.vector(mu) + offset2)

mu2 <- rep(0, length(as.vector(mu)))   # mu2 = 0 * mu (as in C++)
parin <- rep(0, length(as.vector(mu))) # parin = 0 vector (mu - mu)

opt_out <- optim(
  par = parin,
  fn = f2,
  gr = f3,
  y = as.vector(y),
  x = as.matrix(x),
  mu = as.vector(mu2),
  P = as.matrix(P),
  alpha = as.vector(alpha),
  wt = as.vector(wt2_opt),
  method = "BFGS",
  hessian = TRUE
)

bstar <- opt_out$par
A1 <- opt_out$hessian

###############################################################################
# Step C: Standardize model (glmb_Standardize_Model)
###############################################################################
Standard_Mod <- glmb_Standardize_Model(
  y = as.vector(y),
  x = as.matrix(x),
  P = as.matrix(P),
  bstar = as.matrix(bstar, ncol = 1),
  A1 = as.matrix(A1)
)

bstar2 <- Standard_Mod$bstar2
A <- Standard_Mod$A
x2_std <- Standard_Mod$x2
mu2_std <- Standard_Mod$mu2
P2_std <- Standard_Mod$P2

###############################################################################
# Step D: EnvelopeBuild (coefficient envelope at Gridtype = 3)
###############################################################################
max_disp_perc <- 0.99
n_env <- as.integer(200) # used by EnvelopeBuild for diagnostics/overhead
Gridtype_env <- as.integer(3) # EnvelopeOrchestrator overrides to 3

shape2_env <- shape + n_w / 2.0
rate3_env <- rate + RSS_Post2 / 2.0
d1_star <- rate3_env / (shape2_env - 1.0)

wt2_env <- wt / d1_star

Env2 <- EnvelopeBuild(
  bStar = as.vector(bstar2),
  A = as.matrix(A),
  y = as.vector(y),
  x = as.matrix(x2_std),
  mu = as.matrix(mu2_std, ncol = 1),
  P = as.matrix(P2_std),
  alpha = as.vector(alpha),
  wt = as.vector(wt2_env),
  family = "gaussian",
  link = "identity",
  Gridtype = Gridtype_env,
  n = n_env,
  n_envopt = as.integer(1),
  sortgrid = FALSE,
  use_opencl = FALSE,
  verbose = FALSE
)

###############################################################################
# Step E: EnvelopeDispersionBuild (dispersion-aware envelope)
###############################################################################
disp_env_out <- EnvelopeDispersionBuild(
  Env = Env2,
  Shape = shape,
  Rate = rate,
  P = as.matrix(P2_std),
  y = as.vector(y),
  x = as.matrix(x2_std),
  alpha = as.vector(alpha),
  n_obs = as.integer(n_obs),
  RSS_post = RSS_Post2,
  RSS_ML = NA_real_,
  mu = as.matrix(mu2_std, ncol = 1),
  wt = as.vector(wt),
  max_disp_perc = max_disp_perc,
  disp_lower = NULL,
  disp_upper = NULL,
  verbose = FALSE,
  use_parallel = TRUE
)

###############################################################################
# Step F: EnvelopeSort (mirror EnvelopeOrchestrator: disp_grid_type = 2)
###############################################################################
Env3_raw <- disp_env_out$Env_out
UB_list_new <- disp_env_out$UB_list
gamma_list_new <- disp_env_out$gamma_list

cbars <- Env3_raw$cbars
l1 <- ncol(cbars)
l2 <- nrow(cbars)

logP_vec <- Env3_raw$logP
logP_mat <- matrix(logP_vec, nrow = length(logP_vec), ncol = 1)

Env3 <- EnvelopeSort(
  l1 = l1,
  l2 = l2,
  GIndex = Env3_raw$GridIndex,
  G3 = Env3_raw$thetabars,
  cbars = cbars,
  logU = Env3_raw$logU,
  logrt = Env3_raw$logrt,
  loglt = Env3_raw$loglt,
  logP = logP_mat,
  LLconst = Env3_raw$LLconst,
  PLSD = Env3_raw$PLSD,
  a1 = Env3_raw$a1,
  E_draws = Env3_raw$E_draws,
  lg_prob_factor = UB_list_new$lg_prob_factor,
  UB2min = UB_list_new$UB2min
)

UB_list_final <- UB_list_new
UB_list_final$lg_prob_factor <- Env3$lg_prob_factor
UB_list_final$UB2min <- Env3$UB2min

env_final <- list(
  Env = Env3,
  gamma_list = gamma_list_new,
  UB_list = UB_list_final,
  diagnostics = disp_env_out$diagnostics,
  low = gamma_list_new$disp_lower,
  upp = gamma_list_new$disp_upper
)

print(env_final$low)
print(env_final$upp)
print(env_final$gamma_list[c("shape3", "rate2")])

env_final

###############################################################################
# End: envelope construction only
###############################################################################

glmbayes documentation built on Aug. 5, 2026, 1:07 a.m.