EnvelopeOrchestrator: Envelope Construction Orchestrator for Bayesian Gaussian...
In glmbayes: Bayesian Generalized Linear Models (IID Samples)

EnvelopeOrchestrator

R Documentation

Envelope Construction Orchestrator for Bayesian Gaussian Regression

Description

EnvelopeOrchestrator() provides a unified interface for constructing the fixed‑dispersion and dispersion‑aware envelopes used in likelihood‑subgradient simulation for Bayesian Gaussian regression with Normal–Gamma priors.

This function coordinates:

fixed‑dispersion envelope construction via EnvelopeBuild,
dispersion‑refined envelope construction via EnvelopeDispersionBuild,
envelope sorting and reindexing via EnvelopeSort, and
UB‑list alignment (reordered lg_prob_factor and UB2min).

It is typically used inside *.cpp routines such as rIndepNormalGammaReg(), but may also be called directly for diagnostics, envelope visualization, or custom simulation workflows.

Usage

EnvelopeOrchestrator(
  bstar2,
  A,
  y,
  x2,
  mu2,
  P2,
  alpha,
  wt,
  n,
  Gridtype,
  n_envopt,
  shape,
  rate,
  RSS_Post2,
  RSS_ML,
  max_disp_perc,
  disp_lower,
  disp_upper,
  use_parallel = TRUE,
  use_opencl = FALSE,
  verbose = FALSE
)

Arguments

`bstar2`	Numeric vector. Posterior mode of the standardized regression coefficients (from the standardized model).
`A`	Numeric matrix. Posterior precision matrix (Hessian) at the mode.
`y`	Numeric response vector of length `m`.
`x2`	Numeric matrix of standardized predictors (`m × p`).
`mu2`	Numeric vector. Standardized prior mean (typically a zero vector).
`P2`	Numeric matrix. Standardized prior precision component moved into the log‑likelihood.
`alpha`	Numeric vector. Offset‑adjusted mean component.
`wt`	Numeric vector of prior weights.
`n`	Integer. Number of envelope grid points or simulation draws.
`Gridtype`	Integer specifying the envelope grid construction method for API compatibility. The C++ orchestrator overrides this to `3L` (full `3^{p}` grid): unknown dispersion does not use smaller grids.
`n_envopt`	Optional integer. Effective sample size passed to `EnvelopeOpt` during grid construction. Larger values encourage tighter envelopes.
`shape`	Numeric. Shape parameter of the Gamma prior for the dispersion.
`rate`	Numeric. Rate parameter of the Gamma prior for the dispersion.
`RSS_Post2`	Numeric. Expected posterior weighted RSS used to anchor the dispersion axis (typically `centering_out$RSS_post` from `EnvelopeCentering` inside `rindepNormalGamma_reg`; see vignette `Chapter-A11`).
`RSS_ML`	Numeric. Maximum‑likelihood residual sum of squares.
`max_disp_perc`	Numeric in `(0,1)`. Tail probability used to determine dispersion bounds when not explicitly supplied.
`disp_lower`	Optional numeric. Lower bound for the dispersion (`\sigma^2`). If supplied, overrides quantile‑based bounds.
`disp_upper`	Optional numeric. Upper bound for the dispersion (`\sigma^2`). Must be strictly greater than `disp_lower`.
`use_parallel`	Logical. Whether to allow parallel computation inside EnvelopeDispersionBuild.
`use_opencl`	Logical. Whether to allow OpenCL acceleration inside EnvelopeBuild.
`verbose`	Logical. Whether to print detailed progress and timing messages.

Details

EnvelopeOrchestrator() is the envelope-construction stage for Bayesian Gaussian regression with an independent Normal–Gamma prior on (\beta, \phi) (dispersion \phi; precision \tau = 1/\phi in much of the theory). It is implemented in ‘src/EnvelopeOrchestrator.cpp’ and composes direct C++ calls to EnvelopeBuild and EnvelopeDispersionBuild with an R call to EnvelopeSort.

What this function does not do. It does not run the iterative dispersion centering loop (EnvelopeCentering), not optimize the posterior mode or Hessian, not standardize the model (glmb_Standardize_Model), and not draw posterior samples. Those steps are performed by rindepNormalGamma_reg (see vignette Chapter-A11) before and after the orchestrator. Inputs such as bstar2, A, x2, mu2, and P2 must therefore already be in standard form for the coefficient subproblem, exactly as passed from that workflow.

What the return value is for. The returned Env, gamma_list, and UB_list are consumed by the internal standardized samplers rIndepNormalGammaReg_std and rIndepNormalGammaReg_std_parallel in ‘src/rIndepNormalGammaReg.cpp’, which implement the joint accept–reject procedure over (\beta, \phi). Theory for the dispersion envelope and bounding arguments is in vignette Chapter-A07; the end-to-end implementation map is in Chapter-A11. The coefficient-only likelihood-subgradient envelope (\insertCiteNygren2006glmbayes) is documented under EnvelopeBuild and vignette Chapter-A08.

The function does not perform simulation. Simulation is carried out afterward via .rIndepNormalGammaReg_std_cpp() or .rIndepNormalGammaReg_std_parallel_cpp(), depending on use_parallel.

Value

A list with components:

Env: The fully constructed and sorted envelope, including the PLSD component inserted by the dispersion‑aware refinement step.
gamma_list: Updated Gamma‑prior parameters for the dispersion (shape, rate, and dispersion bounds).
UB_list: Updated UB‑list including reordered lg_prob_factor and UB2min.
diagnostics: Diagnostic quantities returned by EnvelopeDispersionBuild, useful for debugging or envelope visualization.
low: Lower dispersion bound used.
upp: Upper dispersion bound used.

Use of the envelope during sampling

After EnvelopeOrchestrator() returns, rindepNormalGamma_reg delegates iid simulation to rIndepNormalGammaReg_std (serial) or rIndepNormalGammaReg_std_parallel (parallel). These routines are not exported; they are the direct analogues of the fixed-dispersion path .rNormalGLM_std_cpp() for GLMs, but for the joint posterior \pi(\beta, \phi \mid y) under the independent Normal–Gamma prior.

Dominating proposal (conceptual). The envelope list Env still describes a mixture of restricted multivariate Normal proposal pieces for the standardized regression coefficients, with mixture weights \tilde{p}_j stored in PLSD. After EnvelopeDispersionBuild, those weights and the per-face constants are adjusted so that, together with a truncated inverse-Gamma (dispersion) proposal derived from gamma_list, the joint proposal dominates the target posterior on the truncated dispersion interval [low, upp]. vignette("Chapter-A07", package = "glmbayes") derives the dispersion-related bounds; vignette("Chapter-A11", package = "glmbayes") records how UB_list entries enter the code.

One accept–reject iteration (standardized coordinates) proceeds as follows:

Draw a mixture component (face) J. An index J is drawn from the discrete distribution with probabilities PLSD.
Propose coefficients \beta^\star. Conditional on J, each coordinate is drawn from the restricted Normal used in the fixed-dispersion construction: cumulative-normal tail probabilities loglt[J, ], logrt[J, ], and subgradient shift -cbars[J, ] (internal rnorm_ct truncated Normal sampling, same structural role as ctrnorm_cpp() in the GLM path).
Propose dispersion \phi. A draw is taken from the truncated inverse-Gamma / Gamma piece defined by shape3, rate2, disp_lower, and disp_upper in gamma_list (rinvgamma_ct_safe).
Re-weight the likelihood for \phi. Observation weights in the Gaussian log-likelihood are scaled by 1/\phi (wt2 = wt / dispersion in the C++ sources).
Dispersion-adjusted tangency. Because the tangency point for the linear upper bound depends on dispersion, the code recomputes a face-specific \bar{\theta}_J(\phi) via Inv_f3_with_disp (using a one-time cache from Inv_f3_precompute_disp built from cbars and the data). The negative log-likelihood at that point feeds the UB1 tangent term.
Log-likelihood at the proposal. Compute -\log f(y \mid \beta^\star, \phi) with the same \phi and scaled weights (output LL_Test in the serial implementation).

Acceptance inequality (structure). Write \ell(\beta,\phi) for the Gaussian log-likelihood (weighted, with offset). The serial sampler forms \mathrm{UB1} from the tangent to -\ell at \bar{\theta}_J(\phi) along subgradient c_J = cbars[J, ]:

\mathrm{UB1} = -\ell\!\big(\bar{\theta}_J(\phi), \phi\big) - c_J^\top \big(\beta^\star - \bar{\theta}_J(\phi)\big).

Additional terms bound RSS variation along \phi (UB2, using RSS_Min and UB2min from UB_list) and dispersion-axis majorization (UB3A, UB3B) built from lg_prob_factor, lmc1, lmc2, lm_log1, lm_log2, max_New_LL_UB, and max_LL_log_disp. With L_{\mathrm{test}} = -\ell(\beta^\star,\phi), define T_1 = L_{\mathrm{test}} - \mathrm{UB1} and T = T_1 - (\mathrm{UB2} + \mathrm{UB3A} + \mathrm{UB3B}). The code draws U_2 \sim \mathrm{Unif}(0,1) and accepts (\beta^\star, \phi) when

T - \log(U_2) \ge 0.

Serial and parallel workers use the same logical decomposition up to implementation detail. Under the construction in EnvelopeDispersionBuild, the terms are arranged so that T_1 \le 0 and \mathrm{UB2}, \mathrm{UB3A}, \mathrm{UB3B} are nonnegative up to controlled numerical slack. Iteration counts are stored in iters_out.

Mapping orchestrator outputs to the sampler.

Env$PLSD: mixture probabilities over envelope faces for Step 1.
Env$loglt, Env$logrt, Env$cbars: restricted Normal proposal for \beta^\star in Step 2.
Env$GridIndex, Env$thetabars, Env$logU, Env$logP: same role as in EnvelopeBuild for the coefficient mixture; dispersion refinement may update PLSD before sorting.
gamma_list: truncated dispersion proposal parameters (shape3, rate2, bounds) for Step 3.
UB_list: global and per-face constants (RSS_Min, UB2min, lg_prob_factor, linear lmc/lm_log pieces) for \mathrm{UB2}, \mathrm{UB3A}, \mathrm{UB3B}.
low, upp: dispersion interval endpoints (duplicated from gamma_list for convenience).

Unlike the fixed-dispersion GLM sampler, this path does not apply the stored LLconst vector directly in the acceptance test; the tangent piece is recomputed as \mathrm{UB1} once \phi and \bar{\theta}_J(\phi) are known.

Algorithmic steps

The orchestrator implements the independent Normal–Gamma envelope pipeline: first a coefficient envelope at a dispersion anchor (\insertCiteNygren2006glmbayes; vignette Chapter-A08), then dispersion-aware refinement (Chapter-A07), then sorting. Steps 3–8 repeat the internal logic of EnvelopeBuild (same formulas on that help page); here the likelihood is Gaussian with identity link, weights are w_i / d_\star with d_\star from the anchor below, and the first pass uses sortgrid = FALSE so sorting runs after dispersion refinement.

Force full grid for unknown dispersion. The argument Gridtype is overridden to 3L so the coefficient grid always uses the full 3^{p} partition (implementation policy in ‘src/EnvelopeOrchestrator.cpp’).
Anchor dispersion and rescale weights for EnvelopeBuild. Let n_w = \sum_i w_i. With prior hyperparameters shape (a_0) and rate (b_0) and centered RSS RSS_Post2, define s = a_0 + n_w/2 and r = b_0 + \mathrm{RSS}_{\mathrm{post}}/2 where \mathrm{RSS}_{\mathrm{post}} denotes RSS_Post2 (the C++ code names the scalars shape2 and rate3). The dispersion anchor is d_\star = r/(s - 1), and observation weights w_i passed into the embedded EnvelopeBuild call are scaled by 1/d_\star. This ties the coefficient envelope to the Gamma posterior for the precision conditional on the centered RSS (Chapters A07, A11).
Compute width parameters \omega_i from the diagonal precision matrix. Let \theta^{\ast} be the standardized posterior mode. For each dimension i,

\omega_{i} := \frac{\sqrt{2} - \exp\!\big(-1.20491 - 0.7321\,\sqrt{0.5 - \partial^{2}\log f(\theta^{\ast}\mid y)/\partial\theta_{i}^{2}}\big)} {\sqrt{1 - \partial^{2}\log f(\theta^{\ast}\mid y)/\partial\theta_{i}^{2}}}.

Here f is the weighted Gaussian log-posterior for \beta at the anchored dispersion.
Construct intervals and the 3^{p} partition around \theta^\star. Set

\ell_{i,1} = \theta^{\ast}_{i} - 0.5\,\omega_{i}, \quad \ell_{i,2} = \theta^{\ast}_{i} + 0.5\,\omega_{i},

and

A_{i,1} = (-\infty,\ell_{i,1}), \quad A_{i,2} = [\ell_{i,1},\ell_{i,2}], \quad A_{i,3} = (\ell_{i,2},\infty).

With J = \prod_{i=1}^{p} \{1,2,3\} and j = (j_1,\ldots,j_p), A^{\ast}_{j} = \prod_{i=1}^{p} A_{i,j_i} partitions standardized coefficient space.
Select tangency points \bar{\theta}_j per cell (left / mode / right of each interval). For index sets C_{j1},C_{j2},C_{j3} by coordinate,

\bar{\theta}_{j,i} = \begin{cases} \theta^{\ast}_{i} - \omega_{i}, & i \in C_{j1}, \\ \theta^{\ast}_{i}, & i \in C_{j2}, \\ \theta^{\ast}_{i} + \omega_{i}, & i \in C_{j3}. \end{cases}
Evaluate negative log-likelihood and gradient at each grid point. Subgradients c(\bar{\theta}_j) and negative log-likelihoods define the likelihood-subgradient envelope pieces (\insertCiteNygren2006glmbayes; Chapter-A08). CPU: f2_f3_non_opencl; GPU (optional): f2_f3_opencl.
Call EnvelopeSet_Grid_C2_pointwise and EnvelopeSet_LogP_C2 (C++ pipeline) to obtain restricted Normal log-densities, mixture log-probabilities, and constants as in Remarks 5–6 of the JASA paper (same as EnvelopeBuild).
Normalize to PLSD without sorting. The embedded EnvelopeBuild call sets sortgrid = FALSE so an intermediate sort is not wasted before EnvelopeDispersionBuild revises mixture weights for the joint (\beta,\phi) target and EnvelopeSort runs once at the end.
Call EnvelopeDispersionBuild (C++). Pass the coefficient envelope list, prior shape, rate, standardized P2, data y, x2, alpha, mu2, wt, RSS_Post2, RSS_ML, dispersion controls (max_disp_perc, optional bounds), and use_parallel. This constructs the dispersion truncation interval, updates Gamma proposal parameters (gamma_list), computes UB_list, and returns Env_out with adjusted PLSD. See Chapter-A07 and Chapter-A11, Section 3.3.
Call EnvelopeSort (R). Reorder envelope components and align lg_prob_factor and UB2min with the sorted indexing. If sorting cannot allocate safely, the implementation falls back to unsorted Env_out with UB fields patched (‘src/EnvelopeOrchestrator.cpp’).
Return Env, gamma_list, UB_list, diagnostics, low, and upp for the standardized samplers.

References

\insertAllCited

Examples

############################### Start of EnvelopeOrchestrator example ####################

# This example demonstrates calling EnvelopeOrchestrator directly for Gaussian
# regression with an independent Normal-Gamma prior. It mirrors the algorithm
# path used inside rIndepNormalGammaReg:
#   - Step A: Initial dispersion via weighted lm.wfit residual variance
#   - Step B: EnvelopeCentering loop (closed-form expected RSS, update
#             dispersion2 via Gamma posterior) to anchor the envelope
#   - Step C: Coefficient posterior mode optimization (optim + f2/f3)
#   - Step D: Standardize the model (glmb_Standardize_Model)
#   - Step E: EnvelopeOrchestrator (EnvelopeBuild + EnvelopeDispersionBuild
#             + EnvelopeSort in one call)
# It stops after envelope construction (no 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/B: EnvelopeCentering (starting at weighted lm.wfit dispersion and
# iteratively refining it via closed-form expected RSS)
###############################################################################
centering <- EnvelopeCentering(
  y = as.vector(y),
  x = as.matrix(x),
  mu = as.vector(mu),
  P = as.matrix(P),
  offset = as.vector(offset2),
  wt = as.vector(wt),
  shape = shape,
  rate = rate,
  Gridtype = Gridtype_core,
  verbose = FALSE
)

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

###############################################################################
# Step C: Coefficient posterior mode optimization (optim + f2/f3)
###############################################################################
dispstar <- dispersion2
wt2_opt <- wt / dispstar
alpha <- as.vector(x %*% as.vector(mu) + offset2)

mu2_opt <- rep(0, length(as.vector(mu)))
parin <- rep(0, length(as.vector(mu)))

opt_out <- optim(
  par = parin,
  fn = f2,
  gr = f3,
  y = as.vector(y),
  x = as.matrix(x),
  mu = as.vector(mu2_opt),
  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 D: 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 E: EnvelopeOrchestrator (EnvelopeBuild + EnvelopeDispersionBuild
#         + EnvelopeSort in one call)
###############################################################################
max_disp_perc <- 0.99
n_env <- as.integer(200)
Gridtype_env <- as.integer(3)  # EnvelopeOrchestrator overrides to 3 for unknown dispersion

env_out <- EnvelopeOrchestrator(
  bstar2 = as.vector(bstar2),
  A = as.matrix(A),
  y = as.vector(y),
  x2 = as.matrix(x2_std),
  mu2 = as.matrix(mu2_std, ncol = 1),
  P2 = as.matrix(P2_std),
  alpha = as.vector(alpha),
  wt = as.vector(wt),
  n = n_env,
  Gridtype = Gridtype_env,
  n_envopt = as.integer(1),
  shape = shape,
  rate = rate,
  RSS_Post2 = RSS_Post2,
  RSS_ML = NA_real_,
  max_disp_perc = max_disp_perc,
  disp_lower = NULL,
  disp_upper = NULL,
  use_parallel = TRUE,
  use_opencl = FALSE,
  verbose = FALSE
)

# Output structure matches that from the step-by-step Ex_EnvelopeDispersionBuild
print(env_out$low)
print(env_out$upp)
print(env_out$gamma_list[c("shape3", "rate2")])

env_out

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

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