Package 'lmebayesCore'

lmebayesCore: Core C++ Sampling Engine for lmebayes

Description

Core C++ engine for envelope-based iid samplers, two-block Gibbs mixed-model engines, and optional OpenCL acceleration. Full-featured developer backend for lmebayes and related extensions. End users should install lmebayes for lmer/glmer-style mixed-model workflows.

Details

Low-level entry points include envelope construction (EnvelopeBuild, EnvelopeOrchestrator), registered simulation pipelines (rNormalGLM_std, rIndepNormalGammaReg_std), matrix-input samplers (rglmb, rlmb), the two-block Gibbs mixed-model engines (rlmerb, rglmerb), and OpenCL kernel loaders. Formula interfaces lmerb and glmerb live in lmebayes; lmebayesCore supplies the sampling engine.

OpenCL startup checks

In interactive sessions, attaching the package with library(lmebayesCore) may emit a short packageStartupMessage when glmbayesCore_has_opencl() is FALSE but a GPU or OpenCL stack appears available on the host. Set options(glmbayes.quiet_opencl_startup = TRUE) to suppress attach notes (recommended for CI and R CMD check).

Author(s)

Maintainer: Kjell Nygren [email protected]

Other contributors:

• The R Core Team (R Mathlib sources, R stats modeling code, and derived/adapted routines) [contributor, copyright holder]

• The R Foundation (Portions of R Mathlib and R source code) [copyright holder]

• Ross Ihaka (R Mathlib and original R modeling infrastructure) [contributor, copyright holder]

• Robert Gentleman (Portions of R Mathlib) [contributor, copyright holder]

• Simon Davies (Original R glm implementation) [contributor]

• Morten Welinder (Portions of R Mathlib) [contributor, copyright holder]

• Martin Maechler (Portions of R Mathlib) [contributor]

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

See Also

lmebayes for the end-user modelling package.

---

Amitriptyline overdose data

Description

Data with information on 17 overdoses of the drug amitriptyline

Usage

data(AMI)

Format

This data frame contains the following columns:

TOT

total TCAD plasma level

AMI

amount of amitriptyline present in the TCAD plasma level

GEN

gender (male = 0, female = 1)

AMT

amount of drug taken at time of overdose

PR

PR wave measurement

DIAP

diastolic blood pressure

QRS

QRS wave measurement

Details

Each row is one overdose episode. Variables include total tricyclic antidepressant level, amitriptyline component, gender, reported amount ingested, and ECG-related measures (PR interval, QRS duration, diastolic blood pressure). The dataset is used in package examples for binomial and related regression; see (Dobson 1990) for analogous generalized linear modelling of clinical outcomes.

References

Dobson A~J (1990). An Introduction to Generalized Linear Models. Chapman and Hall, London.

Examples

############################### Start of AMI dataset example ####################

data(AMI)
summary(AMI)

###############################################################################
## End of AMI dataset example
###############################################################################

---

Bike Sharing Dataset (Processed)

Description

A processed version of the UCI Bike Sharing Dataset (hourly data). The data include derived variables for part of day, quarter, and Fourier terms for cyclic effects of hour and month.

Usage

BikeSharing

Format

A data frame with 17,379 observations and 25 variables (original plus derived):

instant

Record index.

dteday

Date.

season

Season (1: spring, 2: summer, 3: fall, 4: winter).

yr

Year (0: 2011, 1: 2012).

mnth

Month (1–12).

hr

Hour (0–23).

holiday

Whether the day is a holiday (0/1).

weekday

Day of week (0–6).

workingday

Working day (0/1).

weathersit

Weather situation (1–4).

temp

Normalized temperature.

atemp

Normalized feeling temperature.

hum

Normalized humidity.

windspeed

Normalized wind speed.

casual

Count of casual users.

registered

Count of registered users.

cnt

Total count (casual + registered).

hr_num

Hour as numeric (same as hr).

month_num

Month as integer.

part_of_day

Factor: Night (0–5h), Morning (6–11h), Afternoon (12–17h), Evening (18–23h).

quarter

Factor: Q1–Q4.

hr_sin

Sine term for 24-hour cycle.

hr_cos

Cosine term for 24-hour cycle.

mon_sin

Sine term for 12-month cycle.

mon_cos

Cosine term for 12-month cycle.

Source

UCI Machine Learning Repository: Bike Sharing Dataset. https://archive.ics.uci.edu/dataset/275/bike+sharing+dataset

Examples

############################### Start of BikeSharing dataset example ####################

data("BikeSharing")
head(BikeSharing)
dim(BikeSharing)

###############################################################################
## End of BikeSharing dataset example
###############################################################################

---

Conditionally independent block simulation (Gibbs / product likelihood)

Description

Draw from blockwise full conditionals when the posterior factorizes across observation blocks. block_rNormalReg uses .block_rNormalReg_cpp() (Gaussian; each block calls rNormalReg). block_rNormalGLM uses .block_rNormalGLM_cpp() (GLM envelope; C++ partition and prior payload; each block calls rNormalGLM). Typical use is block Gibbs (n = 1 per outer step); n > 1 gives iid draws from the product conditional.

Usage

block_rNormalGLM(
  n,
  y,
  x,
  block,
  prior_list = NULL,
  prior_lists = NULL,
  offset = NULL,
  weights = 1,
  family = gaussian(),
  Gridtype = 2L,
  n_envopt = NULL,
  use_parallel = TRUE,
  use_opencl = FALSE,
  verbose = FALSE,
  progbar = FALSE
)

block_rNormalReg(
  n,
  y,
  x,
  block,
  prior_list = NULL,
  prior_lists = NULL,
  offset = NULL,
  weights = 1,
  Gridtype = 2L
)

block_rNormalReg_update(
  mu_all,
  P = NULL,
  dispersion = NULL,
  y,
  x,
  block,
  prior_list = NULL,
  prior_lists = NULL,
  offset = NULL,
  weights = 1,
  Gridtype = 2L,
  coef_cols = NULL
)

block_rNormalGLM_update(
  mu_all,
  sigma_theta_sq = NULL,
  y,
  x,
  block,
  family = poisson(),
  prior_list = NULL,
  prior_lists = NULL,
  offset = NULL,
  weights = 1,
  Gridtype = 2L,
  n_envopt = 1L,
  use_parallel = TRUE,
  use_opencl = FALSE,
  verbose = FALSE,
  progbar = FALSE,
  theta_coef_col = 1L
)

normalize_block(block, l2)

Arguments

n	Number of iid draws per block (n = 1 typical for Gibbs).
y	Response vector of length nrow(x).
x	Design matrix nrow(x) by ncol(x); same ncol in every block.
block	Block partition: factor or integer vector of length l2, l2_blocks counts summing to l2, or a list of disjoint row-index vectors covering 1:l2.
prior_list	Single prior specification recycled to all blocks, or with mu as l1 by k matrix or blocks sublist.
prior_lists	Optional list of length k (or 1) of per-block prior_list objects.
offset	Optional numeric vector (length 1 or length(y)); partitioned across blocks like y.
weights	Optional weights; same recycling and blocking as offset.
family	GLM family (not gaussian()).
Gridtype	Passed to each block's sampler (Armadillo Gridtype).
n_envopt	Passed to each block; defaults to 1 when NULL.
use_parallel, use_opencl, verbose, progbar	Passed to each block's GLM sampler.
mu_all	Numeric matrix l1 x k or vector of length l1 (recycled to all blocks): prior means per block (e.g. $X_\text{hyper} \gamma$). Required unless prior_lists or prior_list is supplied.
P	Prior precision matrix l1 x l1 (inverse of $\Sigma_{\text{ranef}}$). Shared across all blocks when building per-block priors from mu_all.
dispersion	Residual variance $\sigma^2$. Shared across all blocks when building per-block priors from mu_all.
coef_cols	Column indices of coefficients to return as b_draws (default NULL returns all columns).
sigma_theta_sq	Shared prior variance for scalar blocks when building prior_lists from mu_all.
theta_coef_col	Column index of coefficients to return as theta (default 1 for scalar intercept blocks).
l2	Number of observations (rows) after model.frame.

Details

Output layout: coefficients and coef.mode are matrices with rows = blocks and columns = predictors.

See inst/DESIGN_RGLM_BLOCKS.md.

Per-block prior mean: pass prior_list$mu as an l1 x k matrix (one column per block). Standard for Block~1 of the lmebayes Gibbs sampler, where $\mu_j = X_{\text{hyper}} \gamma$ depends on the current fixed-effects draw.

Dispersion: must be supplied via prior_list$dispersion (residual variance $\sigma^2$).

Value

A list with class "block_rNormalGLM" including:

coefficients

Matrix k * p; row b is the draw for block b.

coef.mode

Matrix k * p; posterior mode per block.

block_info

Block partition metadata.

block_results

List of length k with each block's sampler output.

A list with class "block_rNormalReg" including coefficients, coef.mode, dispersion, and block_info.

A list with b_draws, coefficients, coef.mode, and block_rNormalReg (full block sampler output).

A list with theta, coefficients, coef.mode, and block_rNormalGLM (full block sampler output).

List with k, ids, l2_blocks, starts, and rows (per-block row indices).

Functions

• block_rNormalReg(): Gaussian blockwise full conditionals via .block_rNormalReg_cpp() (C++ partition, prior payload, and rNormalReg() per block). This is the Gaussian counterpart of block_rNormalGLM.

• block_rNormalReg_update(): One Gibbs block update via block_rNormalReg: draw a single blockwise Gaussian posterior sample (n = 1) and return coefficient draws for two-block Gibbs samplers where each group is its own block (e.g. school-level random effects $b_j$). When prior_lists and prior_list are both omitted, per-block dNormal priors are built from mu_all, P, and dispersion.

• block_rNormalGLM_update(): One Gibbs block update via block_rNormalGLM: draw a single blockwise GLM posterior sample (n = 1) and return latent-level draws (e.g. updated theta) for two-block Gibbs samplers where each observation or group is its own block. When prior_lists and prior_list are both omitted, per-block dNormal priors are built from mu_all and sigma_theta_sq.

• normalize_block(): Normalize a row-block partition for BY-style fits and block Gibbs samplers.

See Also

rNormal_reg, simfunction, normalize_block, inst/DESIGN_RGLM_BLOCKS.md

Examples

## Conditionally independent block GLM draw (Poisson, 3 outcome groups)

set.seed(42)

## Dobson (1990) p. 93 RCT
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)
d.AD <- data.frame(outcome, treatment, counts)

ps <- Prior_Setup(counts ~ treatment, family = poisson(), data = d.AD)
y <- ps$y
x <- ps$x
block <- outcome

out <- block_rNormalGLM(
  n = 1L,
  y = y,
  x = x,
  block = block,
  prior_list = list(mu = ps$mu, Sigma = ps$Sigma),
  family = poisson(),
  use_parallel = FALSE
)

out$coefficients
out$coef.mode
## Conditionally independent block Gaussian draw (3 school groups)

set.seed(42)

## Simulate a simple two-level data set: 3 schools, ~10 students each
n_schools <- 3L
n_per     <- 10L
school    <- rep(seq_len(n_schools), each = n_per)
x         <- cbind(1, rnorm(n_schools * n_per))  # intercept + covariate
b_true    <- matrix(c(5, 0.5, 3, -0.2, 7, 0.3), nrow = n_schools, byrow = TRUE)
sigma2    <- 1.5    # residual variance

y <- rowSums(x * b_true[school, ]) + rnorm(nrow(x), sd = sqrt(sigma2))

## Flat prior (large Sigma) shared across all schools; dispersion = sigma2
l1         <- ncol(x)
prior_list <- list(
  mu         = rep(0, l1),
  Sigma      = diag(100, l1),
  dispersion = sigma2,
  ddef       = FALSE
)

out <- block_rNormalReg(
  n          = 1L,
  y          = y,
  x          = x,
  block      = school,
  prior_list = prior_list
)

out$coefficients   ## k x l1 matrix: one row of b_j draws per school
out$coef.mode
## One block_rNormalReg_update step: draw b_j given current fixed effects

set.seed(7)

## Simulate data: 3 schools, 10 students each
n_schools <- 3L
n_per     <- 10L
school    <- rep(seq_len(n_schools), each = n_per)
Z         <- cbind(1, rnorm(n_schools * n_per))  # within-school design
b_true    <- matrix(c(5, 0.5, 3, -0.2, 7, 0.3), nrow = n_schools, byrow = TRUE)
sigma2    <- 1.5

y <- rowSums(Z * b_true[school, ]) + rnorm(nrow(Z), sd = sqrt(sigma2))

## Precision matrix for the random-effect prior (shared across schools)
l1  <- ncol(Z)
P_b <- diag(0.01, l1)   # vague prior

## Current "fixed-effects" prior mean (one vector recycled to all blocks)
mu_current <- rep(0, l1)

out <- block_rNormalReg_update(
  mu_all     = mu_current,   ## recycled to all k blocks (prior_list path)
  P          = P_b,
  dispersion = sigma2,
  y          = y,
  x          = Z,
  block      = school
)

out$b_draws         ## k x l1 matrix of updated random effects
out$coefficients    ## same content (all columns)
## One Gibbs block update: scalar intercept per observation (Poisson)

set.seed(42)

n <- 5L
y <- rpois(n, 2)
x <- matrix(1, n, 1)
mu <- rep(0, n)

upd <- block_rNormalGLM_update(
  mu_all = mu,
  sigma_theta_sq = 1,
  y = y,
  x = x,
  block = seq_len(n),
  family = poisson(),
  use_parallel = FALSE
)

upd$theta
upd$coefficients
## Block partition helpers (factor, l2_blocks counts, or row-index list)

l2 <- 10L

## One level per row (length l2)
normalize_block(factor(rep(c("A", "B"), each = 5L)), l2)

## Contiguous counts summing to l2
normalize_block(c(4L, 6L), l2)

## Explicit row-index list
normalize_block(list(A = 1:4, B = 5:10), l2)

---

Boston housing data with mean-centered predictors

Description

A copy of Boston where all predictors (every column except medv) have been mean-centered (subtract column means, no scaling).

Usage

data("Boston_centered")

Format

A data frame with 506 observations and 14 variables (same names as Boston). See ?MASS::Boston for variable descriptions.

Source

Derived from MASS::Boston. Original data described in Harrison and Rubinfeld (1978); see ?Boston in MASS.

---

Build per-group random-effect prior means for Block 1 sampling

Description

Forms the mu_all matrix passed to block_rNormalReg_update: for each grouping level $j$ and each random-effect column $k$ of the level-1 design Z,

$\mu\_\text{all}[k, j] = X_{\text{hyper},k}[j,]^\top \gamma_k,$

where $\gamma_k$ is the current hyper-parameter vector for RE $k$ (Block 2 state).

Usage

build_mu_all(design, fixef, group_levels = NULL, use_cpp = TRUE)

Arguments

design	List with components X_hyper, re_coef_names, and groups (typically supplied by a downstream mixed-effects model setup step).
fixef	Named list of hyper-parameter vectors, one entry per RE column of Z. Names must match design$re_coef_names. Each fixef[[k]] is a numeric vector of length ncol(X_hyper[[k]]) with names matching colnames(X_hyper[[k]]).
group_levels	Character vector of grouping levels defining the column order of mu_all. Defaults to levels(design$groups), which is the canonical ordering used in two-block mixed models and consistent with lmer and block_rNormalReg (which preserves the input factor's level order).
use_cpp	Logical; if TRUE (default), compute mu_all in C++ via .two_block_build_mu_all_cpp; if FALSE, use the R reference implementation build_mu_all_r().

Value

A list with:

mu_all

Numeric matrix p_re x J (p_re = number of RE columns, J = number of groups). Row names are design$re_coef_names; column names are group_levels in the order supplied.

re_coef_names

Copy of design$re_coef_names.

group_levels

Grouping levels used for columns.

See Also

lmerb_posterior_mean, two_block_rNormal_reg, block_rNormalReg_update

---

Canadian Automobile Insurance Claims for 1957-1958

Description

The data give the Canadian automobile insurance experience for policy years 1956 and 1957 as of June 30, 1959. The data includes virtually every insurance company operating in Canada and was collated by the Statistical Agency (Canadian Underwriters' Association - Statistical Department) acting under instructions from the Superintendent of Insurance. The data given here is for private passenger automobile liability for non-farmers for all of Canada excluding Saskatchewan.

Usage

data(carinsca)

Format

A data frame with 20 observations on the following 6 variables:

Merit

Merit Rating:
3 - licensed and accident free 3 or more years
2 - licensed and accident free 2 years
1 - licensed and accident free 1 year
0 - all others

Class

1 - pleasure, no male operator under 25
2 - pleasure, non-principal male operator under 25
3 - business use
4 - unmarried owner or principal operator under 25
5 - married owner or principal operator under 25

Insured

Earned car years

Premium

Earned premium in 1000's
(adjusted to what the premium would have been had all cars been written at 01 rates)

Claims

Number of claims

Cost

Total cost of the claim in 1000's of dollars

Details

One could apply Poisson regression to the number of claims and gamma regression to the cost per claim.

Source

Bailey, R. A., and Simon, LeRoy J. (1960). Two studies in automobile insurance ratemaking. ASTIN Bulletin, 192-217.

References

Data downloaded from http://www.statsci.org/data/general/carinsca.html. That site also contains classical Poisson and Gamma regression examples.

Examples

############################### Start of carinsca dataset example ####################

data(carinsca)
str(carinsca)
head(carinsca)

###############################################################################
## End of carinsca dataset example
###############################################################################

---

Cleveland Heart Disease Dataset

Description

A cleaned version of the Cleveland heart disease dataset from the UCI Machine Learning Repository. This version contains only complete cases and includes a derived binary outcome variable hd indicating the presence ("Yes") or absence ("No") of heart disease.

Usage

data("Cleveland")

Format

A data frame with 297 observations and 15 variables:

age

Age in years (numeric).

sex

Sex (0 = female, 1 = male).

cp

Chest pain type (numeric code 1–4).

trestbps

Resting blood pressure (mm Hg).

chol

Serum cholesterol (mg/dl).

fbs

Fasting blood sugar > 120 mg/dl (1 = true, 0 = false).

restecg

Resting electrocardiographic results (numeric code).

thalach

Maximum heart rate achieved.

exang

Exercise-induced angina (1 = yes, 0 = no).

oldpeak

ST depression induced by exercise relative to rest.

slope

Slope of the peak exercise ST segment.

ca

Number of major vessels colored by fluoroscopy (0–3).

thal

Thalassemia status (numeric code).

num

Original UCI disease score (0–4).

hd

Binary heart disease indicator: "No" (num = 0) or "Yes" (num > 0).

Source

UCI Machine Learning Repository: Heart Disease Data Set. https://archive.ics.uci.edu/dataset/45/heart+disease

---

Compute Calibrated Gaussian Normal–Gamma Prior Components

Description

Internal Gaussian calibration routine used by Prior_Setup. Given weighted Gaussian regression inputs and a dispersion–independent coefficient–scale prior covariance $\Sigma_0$, this function computes all Normal–Gamma quantities required by the Gaussian prior families: dispersion, shape, shape_ING, rate, rate_gamma, and the calibrated coefficient–scale covariance Sigma. The input Sigma_0 is returned unchanged.

Usage

compute_gaussian_prior(
  X,
  Y,
  weights,
  offset,
  dispersion = NULL,
  n_effective,
  bhat,
  mu,
  Sigma_0,
  Sigma = NULL,
  n_prior,
  k = 1
)

Arguments

X	Numeric model matrix with nrow(X) == length(Y).
Y	Numeric response vector.
weights	Numeric vector of case weights.
offset	Numeric offset vector.
dispersion	Optional scalar dispersion. If supplied, overrides the calibrated value.
n_effective	Effective sample size (typically sum(weights)).
bhat	Numeric coefficient vector, usually the weighted least–squares estimate.
mu	Numeric prior mean vector.
Sigma_0	Dispersion–independent prior covariance matrix $[p \times p]$.
Sigma	Optional coefficient–scale covariance matrix. If supplied, overrides the calibrated Sigma.
n_prior	Effective prior sample size.
k	Non–negative scalar with $k + p \ge 2$. Controls tail behavior of the variance prior; does not affect posterior means.

Details

The computation follows a structured pipeline:

1. Validate dimensions and numeric inputs.

2. Compute the weighted residual sum of squares at bhat.

3. Form the weighted Gram matrix $X^\top W X$, invert it, and construct the marginal quadratic term $S_{marg}$.

4. Map n_prior and k to the Normal–Gamma shape and rate.

5. Calibrate the implied dispersion and coefficient covariance.

6. Return all calibrated prior components.

The function assumes the Chapter 11 convention that $\Sigma_0$ is dispersion–free: it encodes prior structure on the precision–weighted coefficient scale. The returned Sigma is the corresponding coefficient–scale covariance after calibration.

A common choice is the Zellner–type form

$\Sigma_0 = \frac{1 - \mathrm{pwt}}{\mathrm{pwt}} (X^\top W X)^{-1},$

where $\mathrm{pwt}$ is a scalar prior weight. More generally, Sigma_0 may be any positive–definite matrix.

The function computes:

• The marginal quadratic term

  $S_{marg} = RSS_w + (\hat\beta - \mu)^\top \left(\Sigma_0 + (X^\top W X)^{-1}\right)^{-1} (\hat\beta - \mu).$

• Prior Gamma shape: $a_0 = (n_{\mathrm{prior}} + k)/2$.

• Posterior Gamma shape: $a_n = (n_{\mathrm{prior}} + n_w + k)/2$, where $n_w = n_{\mathrm{effective}}$.

• Calibrated dispersion:

  $E[\sigma^2 \mid y] = \frac{S_{marg}}{n_w - p},$

  the usual weighted residual–df estimator.

• Prior Gamma rate:

  $b_0 = \frac{1}{2} S_{marg} \frac{n_{\mathrm{prior}} + k + p - 2}{n_w - p},$

  ensuring $E[\sigma^2 \mid y] = S_{marg}/(n_w - p)$.

• Calibrated coefficient covariance:

  $\Sigma = \frac{n_w}{n_{\mathrm{prior}}} E[\sigma^2 \mid y] (X^\top W X)^{-1}.$

Limiting behavior.

• As $n_{\mathrm{prior}} \to \infty$, the prior becomes increasingly concentrated and dominates the likelihood.

• As $n_{\mathrm{prior}} \to 0^+$ with $n_w > p$, $S_{marg} \to RSS_w$ and $E[\sigma^2 \mid y] \to RSS_w/(n_w - p)$, matching the classical weighted Gaussian estimator.

• Strict positivity of the prior rate requires $n_{\mathrm{prior}} + k + p > 2$.

• The prior mean $\mu$ is never altered; this function calibrates only scale parameters.

Value

A list with components:

• dispersion — calibrated Gaussian dispersion.

• shape — Gamma shape for residual precision.

• shape_ING — shape for the independent Normal–Gamma prior.

• rate — Gamma rate for residual precision.

• rate_gamma — Gamma rate for the fixed–$\beta$ path (dGamma).

• Sigma — calibrated coefficient–scale covariance.

• Sigma_0 — the input dispersion–free covariance.

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

---

Build a named list of dGamma measurement-dispersion priors

Description

Generic for constructing a named list of dGamma prior objects from a prior-specification object. For mixed models, one dGamma() per group level supplies Block~1 observation-level $\sigma^2_j$ priors for lmerb(..., dispersion_ranef = ...).

Usage

dGamma_list(object, ...)

Arguments

object	A prior-specification object.
...	Additional arguments passed to methods (e.g. max_disp_perc, disp_upper_anchor).

Value

A named list whose elements are objects of class "pfamily".

See Also

dGamma, Prior_Setup_lmebayes, pfamily_list

---

Build per-group dGamma priors from a Prior_Setup_lmebayes object

Description

Converts the per-group Block~1 measurement-dispersion calibration stored in a Prior_Setup_lmebayes object into a named list of dGamma pfamily objects, one per group level.

Usage

## S3 method for class 'lmebayes_prior_setup'
dGamma_list(
  object,
  max_disp_perc = NULL,
  disp_center = c("sigma2_hat", "dispersion2"),
  disp_upper_anchor = c("blup", "symmetric"),
  warn_asymmetric = NULL,
  print_asymmetric = NULL,
  asymmetric_R_lo = NULL,
  asymmetric_R_hi = NULL,
  n_rss_iter = 10L,
  ...
)

Arguments

object	An object of class "lmebayes_prior_setup" as returned by Prior_Setup_lmebayes (Gaussian models only).
max_disp_perc	Scalar in $(0.5, 1)$; defaults to object$max_disp_perc.
disp_center	Character: "sigma2_hat" (default) mean-matches the truncation window at each group's calibrated point estimate (disp_upper_anchor controls upper-tail widening); "dispersion2" mean-matches symmetrically at an envelope-centering dispersion estimate that integrates over random-effect uncertainty instead (Part III of inst/DGAMMA_LIST_MARGINAL_AND_BOUNDS.md; ignores disp_upper_anchor).
disp_upper_anchor	Character: "blup" (default) scales the upper-bound rate by $\mathrm{RSS}_{\mathrm{blup}}/\mathrm{RSS}_{\mathrm{ols}}$ per group; "symmetric" uses the same rate as the lower bound. Ignored when disp_center = "dispersion2".
warn_asymmetric	If TRUE (default), warn when any group's truncation window is flagged as asymmetric. Defaults to option glmbayesCore.dgamma_window_warn when NULL.
print_asymmetric	If TRUE, print a table of flagged groups after the warning. Defaults to option glmbayesCore.dgamma_window_print when NULL.
asymmetric_R_lo	Flag when R_lo = lo_pct_BLUP/lo_pct_OLS is below this threshold (default 0.25, or option glmbayesCore.dgamma_window_R_lo).
asymmetric_R_hi	Flag when R_hi = hi_pct_OLS/hi_pct_BLUP exceeds this threshold (default 4, or option glmbayesCore.dgamma_window_R_hi).
n_rss_iter	Positive integer; number of fixed-point iterations for the disp_center = "dispersion2" envelope-centering estimate (default 10L). Ignored when disp_center = "sigma2_hat".
...	Currently ignored.

Details

Prior density (shape_ING, rate from $S_{\mathrm{marg}}$ A12 3.3.4) comes from object$ing_prior_measurement_group (calibrated in Prior_Setup_lmebayes() via compute_gaussian_prior with shared population sd_tau). Truncation bounds use an approximate posterior at $n_{\mathrm{combined},j} = n_{\mathrm{prior},j} + n_j$, mean-matched at one of two per-group centers (disp_center):

"sigma2_hat" (default)

disp_lower is OLS-anchored at sigma2_hat_j; disp_upper widens by the BLUP/OLS residual RSS inflation ratio when disp_upper_anchor = "blup" (default).

"dispersion2"

Both bounds are mean-matched, symmetrically, at a dispersion estimate that integrates over the random effect's own posterior uncertainty (an EnvelopeCentering()-style trace correction; disp_upper_anchor is ignored). Tends to produce narrower, better-centered upper tails than "sigma2_hat" for groups with large BLUP/OLS inflation.

See inst/DGAMMA_LIST_MARGINAL_AND_BOUNDS.md (Parts II and III) for the full derivation of both constructions and their empirical validation.

The returned list carries attribute "window_diagnostics": a data frame with cross-percentiles of the asymmetric truncation window and flag asymmetric_window when $R_{\mathrm{lo}} = \mathrm{lo\_pct\_BLUP}/\mathrm{lo\_pct\_OLS} <$ asymmetric_R_lo (default 0.25) or $R_{\mathrm{hi}} = \mathrm{hi\_pct\_OLS}/\mathrm{hi\_pct\_BLUP} >$ asymmetric_R_hi (default 4). A warning is emitted for flagged groups when warn_asymmetric = TRUE.

Value

A named list of "pfamily" objects keyed by group levels, suitable for lmerb(..., dispersion_ranef = dGamma_list(ps)), with attributes "window_diagnostics" (data frame, one row per group), "dispersion_fit" (object$dispersion_fit, the glmmTMB reference fit used for calibration – lmerb()/glmerb() reuse it as their own dispersion_fit instead of re-fitting glmmTMB when dispersion_ranef carries this attribute), and "calibration_source" (object$calibration_source).

See Also

Prior_Setup_lmebayes, dGamma_list, dGamma

Examples

if (requireNamespace("bayesrules", quietly = TRUE)) {
  data(big_word_club, package = "bayesrules")
  dat <- big_word_club
  dat$school_id <- factor(dat$school_id)
  dat <- subset(dat, !is.na(score_ppvt))

  ps <- Prior_Setup_lmebayes(
    score_ppvt ~ private_school + (1 | school_id),
    data = dat,
    pwt_measurement = 0.01,
    dispformula = ~school_id
  )
  disp_pf <- dGamma_list(ps)
  print(disp_pf[[1L]])
}

---

GPU and OpenCL diagnostics for glmbayes

Description

Compile-time OpenCL probing for glmbayes, plus diagnose_glmbayes() — a readable report that combines opencltools host/runtime checks with this package's OpenCL build status (glmbayesCore_has_opencl).

Workstation probes (GPU vendor detection, drivers, ICD/PATH, and related helpers) live in opencltools; call opencltools::… or see ?opencltools.

Usage

diagnose_glmbayes()

glmbayesCore_has_opencl()

Details

GPU acceleration speeds up envelope construction and grid evaluation when you pass use_opencl = TRUE in rglmb, rlmb, and related functions. CPU-only builds remain fully usable for standard modelling.

Start with diagnose_glmbayes() for a single readable report; use glmbayesCore_has_opencl() for a quick boolean when scripting. Full install notes: vignette("Chapter-16", package = "glmbayes") ((Nygren 2025)).

Diagnostics exported from glmbayes

• diagnose_glmbayes() — full report including compile-time OpenCL status for this package.

• glmbayesCore_has_opencl() — TRUE if this build was compiled with OpenCL support.

• opencltools::has_opencl() — compile-time flag for opencltools (distinct).

Host / runtime checks (opencltools)

• detect_environment_and_gpus()

• detect_compute_runtimes()

• verify_opencl_runtime()

• check_runtime_env()

• get_opencl_core_count()

• load_kernel_source(), load_kernel_library() (pass package = "lmebayesCore" for kernels under inst/cl/)

• add_to_path_windows() and related PATH helpers

References

Nygren K (2025). “Chapter 16: Large models — GPU acceleration using OpenCL.” Vignette in the glmbayes R package. R vignette name: Chapter-16.

See Also

diagnose_glmbayes, glmbayesCore_has_opencl, opencltools, rglmb, rlmb.

---

GPU-Accelerated Envelope Construction for Posterior Simulation

Description

GPU-Accelerated Envelope Construction for Posterior Simulation

Usage

EnvelopeBuild(bStar,A,y,x,mu,P,alpha,wt,family = "binomial",link = "logit", 
Gridtype = 2L,n = 1L,n_envopt=NULL,sortgrid = FALSE,use_opencl = FALSE,verbose = FALSE)

EnvelopeSetGrid(GridIndex, cbars, Lint)

EnvelopeSetLogP(logP, NegLL, cbars, G3)

Arguments

bStar	Point at which envelope should be centered (typically posterior mode).
A	Diagonal precision matrix for the log-likelihood in standard form.
y	A vector of observations of length m.
x	A design matrix of dimension m * p.
mu	A vector giving the prior means of the variables.
P	Prior precision matrix of the variables (positive-definite).
alpha	Offset vector.
wt	A vector of weights.
family	Family for the envelope: binomial, quasibinomial, poisson, quasipoisson, or Gamma.
link	Link function ("logit", "probit", "cloglog" for binomial; "log" for Poisson/Gamma).
Gridtype	Method to determine the number of subgradient densities in the grid.
n	Number of draws from the posterior (used for grid sizing).
n_envopt	Effective sample size passed to EnvelopeOpt for grid construction. Defaults to match n. Larger values encourage tighter envelopes.
sortgrid	Logical; if TRUE, sort the envelope descending by component probability.
use_opencl	Logical; if TRUE, use OpenCL for gradient evaluations.
verbose	Logical; if TRUE, print progress messages.
GridIndex	A matrix indicating, for each grid component, whether the component lies in the left tail, center, or right tail of the density. Rows correspond to grid components; columns correspond to standardized variables.
cbars	A matrix containing the subgradient of the (adjusted) negative log-likelihood at each grid component.
Lint	A matrix storing the lower and upper bounds for each grid component, depending on whether sampling is from the left, center, or right.
logP	A matrix (typically two columns) with information for each grid component. The first column usually holds the output from EnvelopeSetGrid(), corresponding to the restricted normal density.
NegLL	A vector of negative log-likelihood evaluations at each grid component.
G3	A matrix of tangency points used in the grid.

Details

Constructs an enveloping function for posterior simulation using a grid of tangency points. The envelope is used in accept-reject sampling to guarantee iid draws from the posterior distribution. The implementation follows (Nygren and Nygren 2006), with extensions for GPU acceleration (via OpenCL), dynamic grid optimization, and parallelized evaluation.

The envelope is typically built around the posterior mode $\theta^\star$ for a model in standard form (which in this context means a model with a diagonal posterior precision matrix and prior identity precision matrix - glmb_Standardize_Model). It uses dimension-specific width parameters $\omega_i$ derived from the precision matrix. Tangency points are selected per dimension, and the full grid is formed via Cartesian expansion. Negative log-likelihood and gradient values are computed at each grid point, either on CPU or GPU depending on the use_opencl flag. These values are used to construct a piecewise envelope function that dominates the posterior density.

Value

EnvelopeBuild()

A list of envelope components used for accept-reject sampling:

GridIndex

Integer matrix encoding sampling type (tail, center, line) per dimension and region.

thetabars

Matrix of tangency points $\bar{\theta}_j$ for each grid region.

cbars

Matrix of subgradients $c(\bar{\theta}_j)$ of the negative log-likelihood at tangency.

loglt

Matrix of log left-tail probabilities per dimension and region.

logrt

Matrix of log right-tail probabilities per dimension and region.

logU

Matrix of selected per-dimension log-density contributions (tail/center) for each region.

logP

Matrix of total log-probabilities per region (first column); used to derive mixture weights.

PLSD

Vector of normalized mixture weights over grid regions used to draw region indices.

LLconst

Vector of acceptance-test constants per region used in the inequality for rejection sampling.

EnvelopeSetGrid()

A list of matrices computed for grid-based log-density evaluation:

Down

Lower bounds for truncated-normal evaluation per dimension and region.

Up

Upper bounds for truncated-normal evaluation per dimension and region.

lglt

Log left-tail probabilities (from $(-\infty, \mathrm{Up}]$) per dimension and region.

lgrt

Log right-tail probabilities (from $[\mathrm{Down}, \infty)$) per dimension and region.

lgct

Log central-interval probabilities (from $[\mathrm{Down}, \mathrm{Up}]$) per dimension and region.

logU

Selected log-probability per grid cell based on GridIndex (tail or center).

logP

Matrix with row-wise sums of logU (first column) used to form mixture weights.

EnvelopeSetLogP()

A list with updated mixture-weight and acceptance constants:

logP

Input logP with its second column populated by the log of unnormalized visit probabilities per region (mixture denominators).

LLconst

Vector of acceptance constants $-\log f(y \mid \bar{\theta}_j) - c(\bar{\theta}_j)^{T}\bar{\theta}_j$ used in the accept-reject test.

Models in standard form

The standard-form restriction and its closed-form truncated-normal integrals follow (Nygren and Nygren 2006). See (Nygren 2025) for the full theoretical details (standard form restriction, closed-form truncated-normal integrals, and the resulting log-scale tractability).

In the implementation, these standard-form quantities determine the grid-based tangency shifts and the precomputed region constants (e.g., the log-CDF pieces) that drive the mixture weights used by the envelope sampler.

Construction of restricted subgradient densities

For the full restricted density construction and the resulting envelope constants, see (Nygren 2025).

In the implementation, these theory objects become the precomputed region log-constants (via closed-form CDF pieces) that are used to normalize the envelope mixture weights for the accept-reject sampler.

Mixture construction and tractable probabilities

The mixture construction and its tractable region probabilities are derived in (Nygren 2025). In the implementation, these theory objects become the precomputed region log-constants and the mixture weights (PLSD) used by the envelope-based accept-reject sampler.

Log-scale properties of the envelope function

The log-scale form of the envelope factor and the subgradient inequality that imply envelope dominance are given in (Nygren 2025). In the implementation, these properties allow pointwise evaluation in the log-domain and provide the theoretical basis for the rejection test inside the sampler.

Use of the envelope during sampling

The standardized sampler called through .rNormalGLM_std_cpp() uses the envelope to generate posterior samples via rejection sampling. Although not exported, this routine is called internally by .rNormalGLM_cpp(), which in turn is invoked by the user-facing function rNormal_reg(). Together, these routines implement envelope-based sampling for generalized linear models with log-concave likelihood functions and multivariate normal priors.

The envelope provides a mixture of restricted likelihood-subgradient densities, each defined over a region $A_i$, with associated mixture weights $\tilde{p}_i$ stored in PLSD. The sampling proceeds as follows:

1. A region index $J(i)$ is drawn from the discrete distribution defined by PLSD.

2. A candidate $\theta_i$ is drawn from the restricted density $q^{\bar{\theta}_{J(i)}}_{A_{J(i)}}$, using the normal CDF bounds loglt and logrt, and subgradient vector cbars. Simulation for each dimension uses the internal C++ function ctrnorm_cpp(), which explicitly uses these inputs.

3. The log-likelihood $\log f(y \mid \theta_i)$ is computed and stored in testll[0] using the appropriate likelihood function f2.

The acceptance test is performed using the inequality

$\log(U_2) \le \mathrm{LLconst}[J(i)] + \mathrm{cbars}[J(i), ]^{T} \theta_i + \log f(y \mid \theta_i),$

which is equivalent to

$\log(U_2) \le \log f(y \mid \theta_i) - \left( \log f(y \mid \bar{\theta}_{J(i)}) - c(\bar{\theta}_{J(i)})^{T}(\theta_i - \bar{\theta}_{J(i)}) \right),$

where:

• LLconst[J(i)] stores the precomputed quantity $-\log f(y \mid \bar{\theta}_{J(i)}) - c(\bar{\theta}_{J(i)})^{T} \bar{\theta}_{J(i)}$, computed during envelope construction via EnvelopeSet_LogP_C2().

• cbars[J(i), ] is the precomputed subgradient vector $c(\bar{\theta}_{J(i)})$, extracted via cbars(J(i), _). It defines the exponential tilt direction used to evaluate the envelope.

• testll[0] is the log-likelihood at the candidate draw $\theta_i$, evaluated using the model specified by family and link.

• $-\log(U_2)$ is the threshold from a uniform draw $U_2 \sim \mathrm{Unif}(0,1)$.

The right-hand side of this inequality is always non-positive, and equals zero when $\theta_i = \bar{\theta}_{J(i)}$. This reflects the fact that the envelope is tangent to the log-likelihood at each $\bar{\theta}_j$, and lies above it elsewhere.

This procedure guarantees that accepted samples are drawn from the posterior $\pi(\theta \mid y)$. The envelope ensures bounded rejection probability, and the mixture structure allows efficient sampling across regions. The output out contains accepted draws, and draws records the number of attempts per sample.

The components returned by EnvelopeBuild() are used in specific steps of the sampling procedure as follows:

• PLSD is used to randomly select a region index $J(i)$ from the envelope mixture.

• loglt and logrt define the truncated normal bounds for each dimension, used together with cbars to generate candidate values $\theta_i$.

• cbars provides the subgradient vectors $c(\bar{\theta}_j)$ used both for candidate generation and for computing the acceptance test.

• LLconst stores precomputed constants used in the acceptance inequality, avoiding recomputation of posterior terms at tangency points.

• logU stores the per-dimension log-density contributions for each region, computed during envelope setup. These values are summed to produce logP, which determines the mixture weights PLSD.

• logP contains the total log-probabilities for each grid component, which are normalized to form the mixture weights PLSD.

• thetabars stores the tangency points $\bar{\theta}_j$ used to define subgradients and region-specific densities.

• GridIndex encodes the sampling type (tail, center, line) used for each dimension and region, guiding how each coordinate is simulated.

Algorithmic steps (linked to theory)

The implementation of EnvelopeBuild follows the envelope construction in (Nygren and Nygren 2006) for models in standard form (see Section 3–3.3 there). Each computational step corresponds to a theoretical guarantee:

1. Compute width parameters $\omega_i$ from the diagonal precision matrix. In particular, let $\theta^{\ast}$ denote the unique posterior mode. For each dimension $i$, define

   $\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}}}.$

   As seen from the above, the widths $\omega_i$ are derived from the local curvature of the log-likelihood at the posterior mode. This ensures that the three-interval construction per dimension below yields an envelope whose efficiency does not deteriorate with sample size.

2. Use the width parameters to construct intervals around the posterior mode $\theta^\star$. Specifically, we set

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

   and construct three intervals per dimension:

   $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).$

For each dimension $i$, let $J_{i} = \{1,2,3\}$ and define $J = \prod_{i=1}^{p} J_{i}$, which has $3^{p}$ elements. Each $j \in J$ is a vector $(j_{1},\ldots,j_{p})$, and we define

$A^{\ast}_{j} = \prod_{i=1}^{p} A_{i,j_{i}}.$

The collection $A^{\ast} = \{A^{\ast}_{j} : j \in J\}$ forms a partition of $\Theta$.

1. For each member of the partition, select tangency points $\theta^\star \pm \omega_i$.

For each $j \in J$, define index sets

$C_{j1} = \{i : j_{i} = 1\}, \quad C_{j2} = \{i : j_{i} = 2\}, \quad C_{j3} = \{i : j_{i} = 3\}.$

The tangency points $\bar{\theta}_{j}$ are then defined componentwise by

$\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}$

The tangency points are hence chosen so that the envelope touches the log-likelihood at representative points in each interval, guaranteeing dominance and tightness.

1. Build the full grid of tangency points (Cartesian product across dimensions).

   The Cartesian product of per-dimension partitions yields the $3^p$ restricted densities described in the paper, ensuring coverage of the full parameter space.

2. Evaluate negative log-likelihood and gradients at each grid point to construct the likelihood subgradient densities and to facilitate accept rejection sampling

   The subgradients $c(\bar{\theta})$ enter the likelihood-subgradient density construction ((Nygren and Nygren 2006); see also (Nygren 2025)), and both subgradients and negative log-likelihoods (through $h_{\bar{\theta}}(\cdot)$) are used in the accept-reject procedure. CPU and GPU routines compute these values efficiently.

   • On CPU: via f2_f3_non_opencl.

   • On GPU: via f2_f3_opencl, which computes these in parallel across faces

3. Call EnvelopeSet_Grid_C2_pointwise to evaluate restricted multivariate normal log-densities. Each restricted density corresponds to a subset of the partition, normalized as in Remark 5 of (Nygren and Nygren 2006).

4. Call EnvelopeSet_LogP_C2 to compute component log-probabilities and constants. The constants $\tilde{a}$ and mixture weights $\tilde{p}_i$ are computed explicitly as in Remark 6 of the paper, ensuring that the mixture envelope is properly normalized.

5. Normalize probabilities (PLSD) and optionally sort grid components. Normalization implements Claim 2 of the paper so the mixture forms a valid dominating density for the posterior. Sorting is an implementation detail to improve sampling efficiency.

Theory reference (JASA paper and vignette)

Definitions, claims, theorems, remarks, and examples through Remark 16 (including standard form, the $3^p$ partition, and sampling remarks) are in (Nygren and Nygren 2006). An expanded narrative is in vignette("Chapter-A08", package = "glmbayes").

Subgradient density formulation

Each grid component corresponds to a tilted multivariate normal density, normalized using the moment-generating function (MGF). In the single-point case, centered at the posterior mode $\theta^\star$, the density is:

$f(\theta) = \frac{1}{(2\pi)^{p/2} |A|^{-1/2} \cdot \text{MGF}_A(c)} \exp\left( -\frac{1}{2} (\theta - \mu)^T A (\theta - \mu) + c^T (\theta - \theta^\star) \right)$

where:

• $A$ is the precision matrix,

• $\mu$ is the prior mean vector,

• $c$ is the gradient of the log-likelihood at $\theta^\star$,

• $\text{MGF}_A(c)$ is the moment-generating function:

  $\text{MGF}_A(c) = \exp\left( \frac{1}{2} c^T A^{-1} c \right)$

This closed-form density dominates the posterior locally and is used when Gridtype = 1. For richer envelopes, multiple such components are constructed at tangency points $\theta_j$, each with its own gradient $c_j$, and combined into a mixture:

$f_{\text{env}}(\theta) = \sum_{j=1}^{K} p_j f_j(\theta)$

where the weights $p_j$ are computed using log-CDF differences and constants:

$\log p_j = \log \Phi(U_j) - \log \Phi(L_j) - \text{NegLL}_j + \text{LLconst}_j$

Gridtype logic

The Gridtype argument controls how many tangency points are used per dimension:

• 1: Threshold rule. If $1 + a_i \le 2/\sqrt{\pi}$, use a single-point envelope at the mode; otherwise use three points.

• 2: Dynamic optimization via EnvelopeOpt, which balances grid build cost and expected acceptance rate. Grid size is scaled by n and the number of OpenCL cores when GPU is enabled.

• 3: Always use three points per dimension.

• 4: Always use a single point (mode only).

Supported families and links

The following families and link functions are supported:

• Binomial: logit, probit, cloglog

• Quasibinomial: logit, probit

• Poisson: log

• Quasipoisson: log

• Gamma: log

• Gaussian: identity

GPU acceleration (use_opencl = TRUE) is available for all of the above except Gaussian, which is always evaluated on CPU.

GPU acceleration

When use_opencl = TRUE, likelihood and gradient evaluations are offloaded to the GPU using OpenCL. This can substantially reduce runtime for high-dimensional models or large grids. Results are mathematically equivalent to the CPU version, but small numerical differences may occur due to floating-point arithmetic. If reproducibility across hardware is critical, prefer the CPU path.

If OpenCL support was not detected at compile time, the flag is ignored and the CPU implementation is used. Diagnostic messages are printed when verbose = TRUE.

Verbose output

When verbose = TRUE, the function prints:

• Grid type, number of draws, OpenCL usage, and detected core count.

• Grid size after expansion.

• Time-stamped messages when entering the grid loop, starting likelihood evaluations, starting gradient evaluations, and invoking GPU kernels.

• Messages when setting grid values, computing log-probabilities, and sorting.

Any constants needed by the sampling are added to a list and returned.

References

Nygren K (2025). “Chapter A08: Overview of Envelope Related Functions.” Vignette in the glmbayes R package. R vignette name: Chapter-A08.

Nygren K~N, Nygren L~M (2006). “Likelihood Subgradient Densities.” Journal of the American Statistical Association, 101(475), 1144–1156. doi:10.1198/016214506000000357.

See Also

EnvelopeSize, EnvelopeEval, EnvelopeSort, glmb_Standardize_Model; rNormal_reg, rglmb, rlmb. Theory and vignettes: (Nygren and Nygren 2006); (Nygren 2025, 2025).

Examples

data(menarche,package="MASS")
Age2=menarche$Age-13

summary(menarche)
plot(Menarche/Total ~ Age, data=menarche)


x<-matrix(as.numeric(1.0),nrow=length(Age2),ncol=2)
x[,2]=Age2

y=menarche$Menarche/menarche$Total
wt=menarche$Total

mu<-matrix(as.numeric(0.0),nrow=2,ncol=1)
mu[2,1]=(log(0.9/0.1)-log(0.5/0.5))/3

V1<-1*diag(as.numeric(2.0))

# 2 standard deviations for prior estimate at age 13 between 0.1 and 0.9
## Specifies uncertainty around the point estimates

V1[1,1]<-((log(0.9/0.1)-log(0.5/0.5))/2)^2 
V1[2,2]=(3*mu[2,1]/2)^2  # Allows slope to be up to 1 times as large as point estimate 

famfunc<-glmbfamfunc(binomial(logit))

f1<-famfunc$f1
f2<-famfunc$f2
f3<-famfunc$f3
f5<-famfunc$f5
f6<-famfunc$f6

dispersion2<-as.numeric(1.0)
start <- mu
offset2=rep(as.numeric(0.0),length(y))
P=solve(V1)
n=1000


###### Adjust weight for dispersion

wt2=wt/dispersion2

######################### Shift mean vector to offset so that adjusted model has 0 mean

alpha=x%*%as.vector(mu)+offset2
mu2=0*as.vector(mu)
P2=P
x2=x


#####  Optimization step to find posterior mode and associated Precision

parin=start-mu

opt_out=optim(parin,f2,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),
              method="BFGS",hessian=TRUE
)

bstar=opt_out$par  ## Posterior mode for adjusted model
bstar
bstar+as.vector(mu)  # mode for actual model
A1=opt_out$hessian # Approximate Precision at mode

## 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=Standard_Mod$x2
mu2=Standard_Mod$mu2
P2=Standard_Mod$P2
L2Inv=Standard_Mod$L2Inv
L3Inv=Standard_Mod$L3Inv

Env2=EnvelopeBuild(as.vector(bstar2), as.matrix(A),y, as.matrix(x2),
as.matrix(mu2,ncol=1),as.matrix(P2),as.vector(alpha),as.vector(wt2),
family="binomial",link="logit",Gridtype=as.integer(3), n=as.integer(n),
sortgrid=TRUE)

## These now seem to match

Env2

---

Envelope Centering for Bayesian Gaussian Regression

Description

EnvelopeCentering() computes an initial dispersion and the expected posterior weighted RSS (closed form under the Normal posterior for coefficients) for use in envelope construction when the dispersion is unknown. The dispersion-anchoring loop updates dispersion from the Gamma posterior using that expected RSS each iteration. This step is typically called inside rIndepNormalGammaReg() before EnvelopeOrchestrator, but may be used directly for diagnostics or custom workflows.

Usage

EnvelopeCentering(
  y,
  x,
  mu,
  P,
  offset,
  wt,
  shape,
  rate,
  Gridtype = 2L,
  verbose = FALSE
)

Arguments

y	Numeric response vector of length m.
x	Numeric design matrix of dimension m * p.
mu	Numeric vector of prior means (length p).
P	Numeric matrix of prior precision (p * p).
offset	Numeric vector of length m. Use rep(0, m) for none.
wt	Numeric vector of prior weights.
shape	Numeric. Shape parameter of the Gamma prior for the dispersion.
rate	Numeric. Rate parameter of the Gamma prior for the dispersion.
Gridtype	Integer. Grid construction method (default 2).
verbose	Logical. Reserved for API compatibility; currently unused in C++.

Details

The function first obtains an initial dispersion via lm.wfit residual variance, then iteratively: (1) computes the expected weighted RSS under the Normal posterior (closed form), (2) updates the dispersion via the Gamma posterior using the expected RSS. The result is used as dispersion2 and RSS_Post2 in downstream envelope construction (e.g., EnvelopeOrchestrator).

This anchors the joint Normal–Gamma accept–reject construction in (Nygren and Nygren 2006); see vignettes Chapter-A07, Chapter-A11, and (Nygren 2025, 2025).

Value

A list with components:

dispersion

Numeric. Anchored dispersion value.

RSS_post

Numeric. Expected posterior weighted RSS (closed form; last iteration).

References

Nygren K (2025). “Chapter A08: Overview of Envelope Related Functions.” Vignette in the glmbayes R package. R vignette name: Chapter-A08.

Nygren K (2025). “Independent Normal–Gamma Regression Sampler.” Vignette in the glmbayes R package. R vignette name: independent-norm-gamma.

Nygren K~N, Nygren L~M (2006). “Likelihood Subgradient Densities.” Journal of the American Statistical Association, 101(475), 1144–1156. doi:10.1198/016214506000000357.

See Also

EnvelopeOrchestrator for envelope construction; EnvelopeBuild, EnvelopeDispersionBuild; rindepNormalGamma_reg for the full simulation routine; rlmb for the user-facing linear-model interface.

Examples

############################### Start of EnvelopeCentering example ####################

# This example demonstrates EnvelopeCentering in isolation. It computes an
# initial dispersion and posterior RSS for use in envelope construction when
# the dispersion is unknown (Gaussian regression with Normal-Gamma prior).
# This is Step A of the full pipeline in Ex_EnvelopeDispersionBuild and
# Ex_rIndepNormalGammaReg_std.

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))

Gridtype_core <- as.integer(2)

###############################################################################
# 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
)

centering$dispersion
centering$RSS_post

###############################################################################
# End of EnvelopeCentering example
###############################################################################

---

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:

1. 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})$.

2. 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.

3. 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$.

4. 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$

5. 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}$

6. 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)$.

7. 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: (Nygren and Nygren 2006); vignettes Chapter-A07, Chapter-A11; (Nygren 2025, 2025).

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

Nygren K (2025). “Chapter A08: Overview of Envelope Related Functions.” Vignette in the glmbayes R package. R vignette name: Chapter-A08.

Nygren K (2025). “Independent Normal–Gamma Regression Sampler.” Vignette in the glmbayes R package. R vignette name: independent-norm-gamma.

Nygren K~N, Nygren L~M (2006). “Likelihood Subgradient Densities.” Journal of the American Statistical Association, 101(475), 1144–1156. doi:10.1198/016214506000000357.

See Also

EnvelopeBuild, EnvelopeOrchestrator, EnvelopeCentering (for obtaining RSS_post and anchored dispersion), rindepNormalGamma_reg, rlmb; rglmb, glmbfamfunc.

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
###############################################################################

---

Evaluate Negative Log-Likelihood and Gradients

Description

EnvelopeEval() evaluates the negative log-likelihood and gradients at a grid of parameter values, optionally using OpenCL acceleration.

Usage

EnvelopeEval(G4, y, x, mu, P, alpha, wt,
                    family, link,
                    use_opencl = FALSE, verbose = FALSE)

Arguments

G4	Numeric matrix of parameter values (parameters * grid points).
y	Numeric response vector.
x	Numeric design matrix.
mu	Numeric matrix of offsets or prior means.
P	Numeric matrix representing the portion of the prior precision shifted into the likelihood.
alpha	Numeric offset vector of length m
wt	Numeric vector of weights.
family	Character string; model family (e.g. "gaussian").
link	Character string; link function (e.g. "identity").
use_opencl	Logical; if TRUE, attempt OpenCL acceleration.
verbose	Logical; if TRUE, print diagnostic output.

Details

The lower-level helpers f2_f3_non_opencl and f2_f3_opencl are internal C++ kernels used by the CPU and OpenCL backends. The internal routine run_opencl_pilot benchmarks OpenCL performance on a pilot subset of the grid to estimate runtime before full evaluation.

These functions implement the grid evaluation logic used in envelope construction for rejection sampling. They make use of the theory described in (Nygren and Nygren 2006) and the general implementation outlined in (Nygren 2025).

The evaluation workflow has several layers: 1. High-level dispatch (EnvelopeEval)

• EnvelopeEval() is the user-facing entry point. It accepts a grid of parameter values (G4) and the data (y, x, mu, P, alpha, wt).

• If the grid is large (>= 14 columns), it first calls run_opencl_pilot to benchmark OpenCL performance and optionally report estimated runtime.

• It then dispatches to either the CPU or GPU backend:

  • If use_opencl = TRUE and the family is not "gaussian", it calls f2_f3_opencl (an internal C++ kernel).

  • Otherwise, it calls f2_f3_non_opencl (the CPU kernel).

2. CPU backend (f2_f3_non_opencl)

• This function evaluates the negative log-likelihood and gradients using standard CPU routines.

• It inspects the family and link arguments and routes to the correct pair of kernels (⁠f2_*⁠ for the likelihood, ⁠f3_*⁠ for the gradient).

• For example:

  • "binomial" with "logit" calls f2_binomial_logit() and f3_binomial_logit().

  • "poisson" calls f2_poisson() and f3_poisson().

  • "gaussian" calls f2_gaussian() and f3_gaussian().

• These kernels ultimately rely on the same C math routines that R itself uses (from the nmath/rmath libraries), ensuring numerical consistency with base R functions like dnorm, dpois, etc.

3. GPU backend (f2_f3_opencl)

• This function mirrors the CPU backend but executes the likelihood and gradient calculations on an OpenCL device (GPU or CPU).

• It flattens the input matrices/vectors and allocates output buffers.

• It then constructs a full OpenCL program by concatenating:

  • a generic OpenCL support header (OPENCL.CL),

  • OpenCL ports of R's rmath, nmath, and dpq libraries,

  • and the family/link-specific kernel source (e.g. f2_f3_binomial_logit.cl).

• The resulting program is compiled and passed to a kernel runner (f2_f3_kernel_runner) which executes the likelihood and gradient calculations in parallel on the device.

• This ensures that the GPU backend produces results consistent with the CPU backend, but can scale to much larger grids efficiently.

4. Pilot timing (run_opencl_pilot)

• This helper runs a small subset of the grid through the OpenCL backend to estimate runtime.

• It is used by EnvelopeEval() to inform users (when verbose = TRUE) whether OpenCL acceleration is likely to be beneficial.

5. Returned values

• All backends return a list with:

  • NegLL: numeric vector of negative log-likelihood values.

  • cbars: numeric matrix of gradients (parameters * grid points).

6. Role of likelihood and gradients in sampling

• The outputs of EnvelopeEval() - the negative log-likelihood values (NegLL) and the gradient matrix (cbars) - are not endpoints in themselves. They form the envelope used in the rejection sampler implemented by internal functions such as .rNormalGLM_std_cpp().

• This routine is called by .rNormalGLM_cpp(), which underlies the user-facing function rNormal_reg(). Together they implement envelope-based posterior sampling for GLMs with log-concave likelihoods and multivariate normal priors.

7. Simulation execution (accept/reject procedure)

The acceptance test is performed using

$\log(U_2) \leq \log f(y \mid \theta_i) - \Big(\log f(y \mid \bar{\theta}_{J(i)}) - c(\bar{\theta}_{J(i)})^T(\theta_i - \bar{\theta}_{J(i)})\Big) \leq 0$

Connections between code and notation:

• The arguments G4 (in EnvelopeEval) and b (in ⁠f2_f3_*⁠) both represent the grid of tangency points $\bar{\theta}_j$.

• The output NegLL corresponds to $-\log f(y \mid \bar{\theta}_{J(i)})$, i.e. the negative log-likelihood evaluated at each tangency point.

• The output cbars corresponds to the subgradient vectors $c(\bar{\theta}_{J(i)})$, which define the tangent hyperplanes used in the envelope construction.

Precomputation for efficiency:

• Both NegLL and cbars are computed once during envelope construction, prior to the simulation stage.

• This means the sampler does not need to recompute likelihoods or gradients at every candidate draw - it simply reuses the stored values (NegLL, cbars, and LLconst) in the acceptance inequality.

This design ensures that the envelope is tangent to the log-likelihood at each $\bar{\theta}_j$, lies above it elsewhere, and that the accept-reject procedure can run efficiently while still producing samples from the true posterior $\pi(\theta \mid y)$.

Value

EnvelopeEval

List with components NegLL (numeric vector of negative log-likelihood values) and cbars (numeric matrix of gradients).

f2_f3_non_opencl

List with components qf (negative log-likelihood) and grad (gradients) from the CPU kernel.

f2_f3_opencl

List with components qf and grad from the OpenCL kernel.

run_opencl_pilot

Numeric scalar giving estimated runtime (seconds) for OpenCL evaluation on a pilot subset of the grid.

References

Nygren K (2025). “Chapter A05: Simulation Methods – Likelihood Subgradient Densities.” Vignette in the glmbayes R package. R vignette name: Chapter-A05.

Nygren K~N, Nygren L~M (2006). “Likelihood Subgradient Densities.” Journal of the American Statistical Association, 101(475), 1144–1156. doi:10.1198/016214506000000357.

See Also

EnvelopeBuild, EnvelopeSize, EnvelopeSort; rNormal_reg, rglmb. Vignettes: (Nygren 2025, 2025).

Examples

############################### Start of EnvelopeEval example ####################

# This example demonstrates EnvelopeEval in isolation. EnvelopeEval evaluates
# the negative log-likelihood and gradients at a grid of parameter values.
# It is called internally by EnvelopeBuild. Here we build the same inputs
# (grid G4, standardized model) using EnvelopeSize and expand.grid, then
# call EnvelopeEval directly. The setup mirrors Ex_EnvelopeBuild through
# the standardization step.

data(menarche, package = "MASS")
Age2 <- menarche$Age - 13

x <- matrix(as.numeric(1.0), nrow = length(Age2), ncol = 2)
x[, 2] <- Age2

y <- menarche$Menarche / menarche$Total
wt <- menarche$Total

mu <- matrix(as.numeric(0.0), nrow = 2, ncol = 1)
mu[2, 1] <- (log(0.9 / 0.1) - log(0.5 / 0.5)) / 3

V1 <- 1 * diag(as.numeric(2.0))
V1[1, 1] <- ((log(0.9 / 0.1) - log(0.5 / 0.5)) / 2)^2
V1[2, 2] <- (3 * mu[2, 1] / 2)^2

famfunc <- glmbfamfunc(binomial(logit))
f2 <- famfunc$f2
f3 <- famfunc$f3

dispersion2 <- as.numeric(1.0)
start <- mu
offset2 <- rep(as.numeric(0.0), length(y))
P <- solve(V1)
n <- 1000

wt2 <- wt / dispersion2
alpha <- x %*% as.vector(mu) + offset2
mu2 <- 0 * as.vector(mu)
P2 <- P
x2 <- x

parin <- start - mu
opt_out <- optim(parin, f2, 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),
  method = "BFGS", hessian = TRUE
)

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

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 <- Standard_Mod$x2
mu2 <- Standard_Mod$mu2
P2 <- Standard_Mod$P2

###############################################################################
# Build grid G4 via EnvelopeSize and expand.grid (as EnvelopeBuild does)
###############################################################################
a <- diag(A)
omega <- (sqrt(2) - exp(-1.20491 - 0.7321 * sqrt(0.5 + a))) / sqrt(1 + a)
b2 <- as.vector(bstar2)
G1 <- rbind(b2 - omega, b2, b2 + omega)

size_info <- EnvelopeSize(a, G1, Gridtype = 3L, n = n)
G2 <- size_info$G2

G3 <- as.matrix(do.call(expand.grid, G2))
G4 <- t(G3)

###############################################################################
# EnvelopeEval: negative log-likelihood and gradients at grid points
###############################################################################
eval_out <- EnvelopeEval(
  G4 = G4,
  y = y,
  x = as.matrix(x2),
  mu = as.matrix(mu2, ncol = 1),
  P = as.matrix(P2),
  alpha = as.vector(alpha),
  wt = as.vector(wt2),
  family = "binomial",
  link = "logit",
  use_opencl = FALSE,
  verbose = FALSE
)

eval_out$NegLL
eval_out$cbars

###############################################################################
# End of EnvelopeEval example
###############################################################################

---