Package 'glmbayes'

glmbayes: Bayesian Generalized Linear Models with iid Sampling

Description

glmbayes provides independent and identically distributed (iid) samples for Bayesian generalized linear models (GLMs), serving as a Bayesian analogue to the base glm() function. Supported likelihood families include Gaussian, Poisson, Binomial, and Gamma models with log-concave likelihoods.

Details

The main user-facing interface is glmb(), which mirrors the structure of glm() and supports prior specification through pfamily objects. Lower-level functions such as rglmb() and rGamma_reg() provide direct access to the underlying samplers and can be used in block Gibbs sampling or hierarchical model implementations.

For an introduction to the package, examples, and a complete set of vignettes, see:

• README: https://github.com/knygren/glmbayes#readme

• All vignettes: browseVignettes("glmbayes")

The package includes extensive documentation on model fitting, prior construction, diagnostics, and optional GPU acceleration using OpenCL.

Releases: This source tree is 0.9.6 (in development). The current CRAN release is 0.9.5 (install.packages("glmbayes")). Source is available from GitHub; R-Universe (https://knygren.r-universe.dev/glmbayes) also builds binaries from that source. Prebuilt CRAN and R-Universe binaries do not include OpenCL; GPU support requires a source install once the host OpenCL environment is ready (see vignette("Chapter-16", "glmbayes") for the three-step process).

IID posterior simulation for non-Gaussian GLMs and several non-conjugate linear-model setups uses the likelihood-subgradient envelope method of (Nygren and Nygren 2006). Introductory material and worked examples are in (Nygren 2025, 2025); estimation and simulation background in (Nygren 2025, 2025); prior derivations for Prior_Setup() in (Nygren 2025); GPU/OpenCL topics in (Nygren 2025, 2025).

OpenCL startup checks

In interactive sessions, attaching the package with library(glmbayes) may emit a short packageStartupMessage when has_opencl() is FALSE (typical for CRAN binaries) but a GPU or OpenCL stack appears available on the host. OpenCL modelling paths require a source install of glmbayes with OpenCL at compile time; has_opencl() then reports whether that build succeeded. The note confirms full CPU use and points to vignette("Chapter-16"). Machines without a detectable GPU stack stay silent. Set options(glmbayes.quiet_opencl_startup = TRUE) to suppress attach notes (recommended for CI and R CMD check).

Author(s)

Kjell Nygren

References

Nygren K (2025). “Chapter 00: Introduction.” Vignette in the glmbayes R package. R vignette name: Chapter-00.

Nygren K (2025). “Chapter A01: A detailed overview of the glmbayes package.” Vignette in the glmbayes R package. R vignette name: Chapter-A01.

Nygren K (2025). “Chapter A02: Overview of Estimation Procedures.” Vignette in the glmbayes R package. R vignette name: Chapter-A02.

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

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

Nygren K (2025). “Chapter A12: Technical Derivations for Priors Returned by Prior_ Setup.” Vignette in the glmbayes R package. R vignette name: Chapter-A12.

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

Nygren K (2025). “Chapter A10: Accelerated EnvelopeBuild Implementation using OpenCL.” Vignette in the glmbayes R package. R vignette name: Chapter-A10.

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

Main interfaces: glmb, lmb, rglmb, rlmb; low-level simulation API simfuncs; envelope construction EnvelopeBuild.

Useful links:

• CRAN: https://CRAN.R-project.org/package=glmbayes

• GitHub: https://github.com/knygren/glmbayes

• R-Universe: https://knygren.r-universe.dev/glmbayes

Examples

set.seed(333)
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))

## Call to glm
glm.D93 <- glm(counts ~ outcome + treatment, 
               family = poisson())

## Using glmb
## Step 1: Set up Prior
ps=Prior_Setup(counts ~ outcome + treatment,family = poisson())
mu=ps$mu
V=ps$Sigma

# Step 2: Call the glmb function
glmb.D93<-glmb(counts ~ outcome + treatment, family=poisson(), 
               pfamily=dNormal(mu=mu,Sigma=V))

summary(glmb.D93)

## ----Printed_Views------------------------------------------------------------
## Printed view of the output from the glm function 
print(glm.D93)
## Printed view of the output from the glmb function 
print(glmb.D93)

## ----Methods---------------------------------------------------------------
## Methods for class "lm"
methods(class="lm")

## Methods for class "glm"
methods(class="glm")

## Methods for class "glmb"
methods(class="glmb")

## ----summary--------------------------------------------------------------
## summary output for the "glm" class
summary(glm.D93)

## summary output for the "glm" class
summary(glmb.D93)

## ----fitted outputs-------------------------------------------------------
## fitted outputs for the glm function
fitted(glm.D93)

## ----glmb fitted outputs------------------------------------------------------
## mean of fitted outputs for the glm function
colMeans(fitted(glmb.D93))

## ----predictions----------------------------------------------------------
## predictions for the glm function
predict(glm.D93)

## predictions for the glmb function
colMeans(predict(glmb.D93)) 

## ----residuals------------------------------------------------------------
## residuals for the glm function
residuals(glm.D93)

## residuals for the glmb function
colMeans(residuals(glmb.D93))

## ----vcov-----------------------------------------------------------------
## vcov for the glm function
vcov(glm.D93)

## vcov for the glm function
vcov(glmb.D93)

## ----confint--------------------------------------------------------------
## confint for the glm function
confint(glm.D93)

## confint for the glm function
confint(glmb.D93)

## ----AIC/DIC------------------------------------------------------------------
## AIC for the glm function (equivalent degrees of freedom and the AIC)
extractAIC(glm.D93)

## DIC for the glmb function (estimated effective number of parameters and the DIC)
extractAIC(glmb.D93)

## ----Deviance-------------------------------------------------------------
## Deviance for the glm function
deviance(glm.D93)

## Deviance for the glmb function
mean(deviance(glmb.D93))

## ----logLik---------------------------------------------------------------
## Deviance for the glm function
logLik(glm.D93)

## Deviance for the glmb function
mean(logLik(glmb.D93))

## ----Model Frame----------------------------------------------------------
## Model Frame for the glm function
model.frame(glm.D93)

## Model Frame for the glmb function
model.frame(glmb.D93$glm)

## ----formula--------------------------------------------------------------
## formula for the glm function
formula(glm.D93)

## ----formula-------------------------------------------------------------
## formula for the glmb function
formula(glmb.D93)

## ----family--------------------------------------------------------------
## family for the glm function
family(glm.D93)

## family for the glmb function
family(glmb.D93$glm)

## ----nobs-----------------------------------------------------------------
## nobs for the glm function
nobs(glm.D93)

## nobs for the glmb function
nobs(glmb.D93)

## ----show-----------------------------------------------------------------
## show for the glm function
show(glm.D93)

## show for the glmb function
show(glmb.D93)

---

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

---

Analysis of Deviance for Bayesian Generalized Linear Model Fits

Description

Compute an analysis of deviance table for one (current implementation) or more (future) Bayesian generalized linear model fits. The structure follows the sequential analysis of deviance (McCullagh and Nelder 1989), with Bayesian extensions for DIC, pD, Mahalanobis shift, and directional tail probability (Spiegelhalter et al. 2002).

Usage

## S3 method for class 'glmb'
anova(object, ...)

Arguments

object	an object of class glmb, typically the result of a call to glmb
...	Other arguments passed to or from other methods.

Details

Specifying a single object (currently only implementation) gives a sequential analysis of deviance table for that fit. The reductions in residual deviance as each term of the formula is added in turn are given as the rows of a table, plus the residual deviances themselves. The Mahalanobis shift and pDirectional columns report prior-posterior disagreement diagnostics; see directional_tail and (Nygren 2025).

Value

An object of class "anova" inheriting from class "data.frame".

References

McCullagh P, Nelder J~A (1989). Generalized Linear Models. Chapman and Hall, London.

Nygren K (2025). “Chapter A04: Directional Tail Diagnostics for Prior-Posterior Disagreement.” Vignette in the glmbayes R package. R vignette name: Chapter-A04.

Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A (2002). “Bayesian Measures of Model Complexity and Fit.” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4), 583–639. doi:10.1111/1467-9868.00353.

See Also

directional_tail, summary.glmb, glmb, glmbayes-package; rglmb, rlmb, lmb; anova.glm

Examples

set.seed(333)
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)

ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson())
glmb.D93 <- glmb(
  counts ~ outcome + treatment,
  family = poisson(),
  pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma)
)
summary(glmb.D93)

# anova for Bayesian Model
# Sequential classical ANOVA (see glm anova below): adding outcome reduces residual deviance;
# adding treatment barely changes it.
# By DIC, the model with outcome but not treatment has the lowest DIC;
# treatment does not improve the criterion.
# Use of traditional anova is questionable in a Bayesian context.
# One may instead use Bayes factors or other approaches.
anova(glmb.D93)


glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
summary(glm.D93)

# corresponding
anova(glm.D93)

---

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

---

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.

Examples

############################### Boston_centered dataset example ####################

data("Boston_centered")
head(Boston_centered)
summary(Boston_centered)

## Predictors are mean-centered (column means ~0)
predictors <- setdiff(names(Boston_centered), "medv")
colMeans(Boston_centered[predictors])

form <- medv ~
  crim + zn +
  indus + chas + nox + age + dis + rad + tax + ptratio + black + lstat + rm

lm.boston <- lm(form, data = Boston_centered, x = TRUE, y = TRUE)
summary(lm.boston)

ps <- Prior_Setup(form, gaussian(), data = Boston_centered)

lmb.boston <- lmb(
  form,
  data = Boston_centered,
  pfamily = dNormal(
    mu = ps$mu,
    Sigma = ps$Sigma,
    dispersion = ps$dispersion
  )
)
summary(lmb.boston)

lmb.boston_v2 <- lmb(
  form,
  data = Boston_centered,
  pfamily = dNormal_Gamma(
    mu = ps$mu,
    Sigma_0 = ps$Sigma_0,
    shape = ps$shape,
    rate = ps$rate
  )
)
summary(lmb.boston_v2)

## Independent Normal-Gamma (OpenCL path when available)

  if (has_opencl()) {
    lmb.boston_v3 <- lmb(
      n = 1000L,
      form,
      data = Boston_centered,
      pfamily = dIndependent_Normal_Gamma(
        ps$mu,
        ps$Sigma,
        shape = ps$shape_ING,
        rate = ps$rate
      ),
      use_parallel = TRUE,
      use_opencl = TRUE,
      verbose = FALSE
    )
    summary(lmb.boston_v3)
  }

###############################################################################
## End of Boston_centered dataset example
###############################################################################

---

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

---

Case and Variable Names of Fitted Models

Description

Simple utilities returning (non-missing) case names, and (non-eliminated) variable names.

Usage

## S3 method for class 'glmb'
case.names(object, full = FALSE, ...)

## S3 method for class 'glmb'
variable.names(object, full = FALSE, ...)

Arguments

object	a fitted model object, typically of class "glmb".
full	logical; if TRUE, all names (including zero weights, ...) are returned.
...	further arguments passed to or from other methods.

Value

A character vector

See Also

glmb, glmbayes-package; rglmb, rlmb, lmb; case.names.

Examples

## ----dobson-------------------------------------------------------------------
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)

ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson())
## Call to glmb
glmb.D93 <- glmb(
  n = 1000,
  counts ~ outcome + treatment,
  family = poisson(),
  pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma)
)

## ----glmb vcov----------------------------------------------------------------
case.names(glmb.D93)

---

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

Examples

############################### Start of Cleveland dataset example ####################

data("Cleveland")
head(Cleveland)
summary(Cleveland)

# OpenCL-accelerated Bayesian logistic regression example
# This only runs if OpenCL is available

if (has_opencl()) {
  ps <- Prior_Setup(
    hd ~ age + sex + cp + trestbps + chol +
      fbs + restecg + thalach + exang + oldpeak + slope + ca + thal,
    family = binomial(logit),
    data = Cleveland
  )

  fit <- glmb(
    hd ~ age + sex + cp + trestbps + chol +
      fbs + restecg + thalach + exang + oldpeak + slope + ca + thal,
    family       = binomial(link = "logit"),
    pfamily      = dNormal(mu = ps$mu, Sigma = ps$Sigma),
    data         = Cleveland,
    n            = 1000,
    Gridtype     = 2,
    use_parallel = TRUE,
    use_opencl   = TRUE,
    verbose      = FALSE
  )
  summary(fit)
}

###############################################################################
## End of Cleveland dataset example
###############################################################################

---

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.

---

Credible Intervals for Model Parameters

Description

Computes Bayesian credible intervals for model parameters from posterior draws. Intervals are constructed as quantiles of the posterior distribution, following the standard approach for posterior summaries (Gelman et al. 2013).

Usage

## S3 method for class 'glmb'
confint(object, parm, level = 0.95, ...)

Arguments

object	a fitted model object of class "glmb". Typically the result of a call to glmb.
parm	a specification (not yet implemented) of which parameters are to be given credible sets, either a vector of numbers or a vector of names. If missing, all parameters are considered.
level	the credible interval required.
...	additional argument(s) for methods.

Value

A matrix (or vector) with columns giving lower and upper credible limits for each parameter. These will be labeled (1-level)/2 and 1-(1-level)/2 in \

References

Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB (2013). Bayesian Data Analysis, 3rd edition. CRC Press.

See Also

summary.glmb, vcov.glmb, glmb, glmbayes-package; rglmb, rlmb, lmb; confint for classical confidence intervals

Examples

## ----dobson-------------------------------------------------------------------
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)

ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson())
## Call to glmb
glmb.D93 <- glmb(
  n = 1000,
  counts ~ outcome + treatment,
  family = poisson(),
  pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma)
)
## ----glmb confint-------------------------------------------------------------
confint(glmb.D93)

---

Model Deviance

Description

Returns the deviance of a fitted Bayesian Generalized Linear Model

Usage

## S3 method for class 'rglmb'
deviance(object, ...)

Arguments

object	an object of class "rglmb", typically the result of a call to rglmb
...	further arguments to or from other methods

Value

A vector with the deviance extracted from the object.

See Also

rglmb; glmb, glmbayes-package; rlmb, lmb; deviance.

Examples

## ----dobson-------------------------------------------------------------------
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)

ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson())
## Call to glmb
glmb.D93 <- glmb(
  n = 1000,
  counts ~ outcome + treatment,
  family = poisson(),
  pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma)
)
## ----rglmb deviance-------------------------------------------------------------
deviance(glmb.D93)

---

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 (has_opencl).

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

Usage

diagnose_glmbayes()

has_opencl()

Details

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

Start with diagnose_glmbayes() for a single readable report; use 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.

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

Host / runtime checks (opencltools)

• detect_environment_and_gpus()

• detect_compute_runtimes()

• verify_opencl_runtime()

• check_runtime_env()

• 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, has_opencl, opencltools, glmb, rglmb.

---

Directional Tail Diagnostic

Description

Computes the directional tail probability based on posterior draws and prior mean, using whitening transformation and projection onto the direction of disagreement. This diagnostic identifies directional disagreement between posterior and prior, and is especially useful for visualizing rejection regions in whitened space. The whitening uses Mahalanobis distance (Mahalanobis 1936) in posterior-precision-scaled coordinates.

Usage

directional_tail(fit, mu0 = NULL)

## S3 method for class 'directional_tail'
print(x, ...)

Arguments

fit	A fitted model object of class 'glmb' or 'lmb'
mu0	An optional argument containing a reference vector relative to which the directional tail is computed. Defaults to the prior mean.
x	An object of class directional_tail
...	Additional arguments passed to or from other methods.

Details

Whitening is performed using the posterior precision matrix. The direction vector is computed as the mean shift in whitened space. Tail probability is the proportion of draws with negative projection onto this direction. For theory, interpretation, and relation to t/F statistics, see (Nygren 2025).

Value

An object of class 'directional_tail' containing:

mahalanobis_shift	Measures the standardized Mahalanobis distance between the posterior and prior means, using posterior precision for scaling. In the Gaussian case, this directly determines the directional tail probability via Phi(-||w||).
p_directional	Directional tail probability (proportion of draws in the direction of disagreement)
delta	Mean shift in whitened space
draws	List containing whitened draws, raw draws, and tail flags

References

Mahalanobis PC (1936). “On the generalized distance in statistics.” Proceedings of the National Institute of Sciences of India, 2(1), 49–55.

Nygren K (2025). “Chapter A04: Directional Tail Diagnostics for Prior-Posterior Disagreement.” Vignette in the glmbayes R package. R vignette name: Chapter-A04.

See Also

summary.glmb, anova.glmb

Examples

## Annette Dobson (1990) "An Introduction to Generalized Linear Models".
## Page 9: Plant Weight Data.
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)

dat  <- data.frame(weight, group)
dat2 <- rbind(dat)

lm.D9_null <- lm(weight ~ 1, data = dat2)
summary(lm.D9_null)

lm.D9_default <- lm(weight ~ group, data = dat2)
summary(lm.D9_default)

ps_null <- Prior_Setup(
  weight ~ group,
  family = gaussian(),
  pwt = 0.01,
  intercept_source = "null_model",
  effects_source   = "null_effects",
  data = dat2
)
mu_null     <- ps_null$mu
V_null      <- ps_null$Sigma
disp_ML_null <- ps_null$dispersion

lmb.D9_null <- lmb(
  weight ~ group,
  dNormal(mu_null, V_null, dispersion = disp_ML_null),
  data = dat2,
  n = 10000
)
summary(lmb.D9_null)

glmb.D9_default <- glmb(
  weight ~ group,
  family  = gaussian(),
  pfamily = dNormal(mu_null, V_null, dispersion = disp_ML_null),
  data    = dat2
)
summary(glmb.D9_default)

colMeans(residuals(lmb.D9_null))

glm.D9_default <- glm(weight ~ group, family = gaussian(), data = dat2)
summary(glm.D9_default)
disp_D9_default <- summary(glm.D9_default)$dispersion

solve(vcov(lmb.D9_null$lm))
t(lmb.D9_null$lm$x) %*% lmb.D9_null$lm$x / disp_D9_default
t(lmb.D9_null$lm$x) %*% lmb.D9_null$lm$x * (20 / 18)

tailprobs <- directional_tail(lmb.D9_null)
tailprobs

summary(lm.D9_null)
summary(lm.D9_null)$sigma^2

tailprobs$Prec_lik * summary(lm.D9_null)$sigma^2
t(lmb.D9_null$x) %*% lmb.D9_null$x
tailprobs$p_directional

##########################################

Z     <- tailprobs$draws$Z
flag  <- tailprobs$draws$is_tail
delta <- tailprobs$delta
w     <- tailprobs$delta

asp <- 1

## Plot posterior draws in whitened space
plot(
  Z,
  col = ifelse(flag, "red", "blue"),
  pch = 19,
  xlab = "Z1",
  ylab = "Z2",
  main = "Directional Tail Diagnostic"
)

abline(a = 0, b = -w[1] / w[2], col = "darkgreen", lty = 2)

## Add radius boundary centered at posterior mode (delta)
r <- sqrt(sum(delta^2))
symbols(
  delta[1], delta[2],
  circles = r,
  inches  = FALSE,
  add     = TRUE,
  lwd     = 2,
  fg      = "gray"
)

## Add prior mean at origin
points(0, 0, pch = 4, col = "black", lwd = 2)

## Add posterior mode at delta
points(delta[1], delta[2], pch = 3, col = "purple", lwd = 2)

## Add legend
legend(
  "topright",
  legend = c(
    "Tail draws", "Non-tail draws", "Direction vector",
    "Radius boundary", "Prior mean", "Posterior mode"
  ),
  col = c("red", "blue", "darkgreen", "gray", "black", "purple"),
  pch = c(19, 19, NA, NA, 4, 3),
  lty = c(NA, NA, 1, 1, NA, NA),
  lwd = c(NA, NA, 2, 2, 2, 2),
  bty = "n"
)

############################ Original Scales ###################

B       <- tailprobs$draws$B
flag    <- tailprobs$draws$is_tail
mu0     <- as.numeric(lmb.D9_null$Prior$mean)
mu_post <- colMeans(B)
x_range <- range(B[, 1])  # Intercept values
padding <- diff(x_range) * 0.1  # 10% margin

oldpar <- par(no.readonly = TRUE)
par(mar = c(5, 6, 4, 2))  # bottom, left, top, right

plot(
  B,
  col  = ifelse(flag, "red", "blue"),
  pch  = 19,
  xlab = "Intercept",
  ylab = "groupTrt",
  xlim = c(x_range[1] - padding, x_range[2] + padding),
  main = "Directional Tail Diagnostic (Raw Space)"
)

points(mu0[1], mu0[2], pch = 4, col = "black", cex = 1.5)       # Prior mean
points(mu_post[1], mu_post[2], pch = 3, col = "darkgreen", cex = 1.5)  # Posterior mean

legend(
  "topright",
  legend = c("Tail draws", "Non-tail draws", "Prior", "Posterior"),
  col    = c("red", "blue", "black", "darkgreen"),
  pch    = c(19, 19, 4, 3)
)

par(oldpar)

---

Extract Coefficients in Original Coding

Description

Extracts coefficients in terms of the original factor levels rather than the coded variables. The logic mirrors dummy.coef for lm objects (Chambers 1992; Venables and Ripley 2002).

Usage

## S3 method for class 'glmb'
dummy.coef(object, use.na = FALSE, ...)

## S3 method for class 'dummy_coef.glmb'
print(x, ...)

Arguments

object	a glmb model fit
use.na	logical flag for coefficients in a singular model. If use.na is true, undetermined coefficients will be missing; if false they will get one possible value.
x	object to be printed
...	arguments passed to or from other methods

Details

A fitted linear model has coefficients for the contrasts of the factor terms, usually one less in number than the number of levels. This function re-expresses the coefficients in the original coding; as the coefficients will have been fitted in the reduced basis, any implied constraints (e.g., zero sum for contr.helmert or contr.sum) will be respected. There will be little point in using dummy.coef for contr.treatment contrasts, as the missing coefficients are by definition zero.

Value

A list giving for each term the draws for the coefficients.

References

Chambers JM (1992). “Linear Models.” In Chambers JM, Hastie TJ (eds.), Statistical Models in S, chapter 4, 85–124. Wadsworth & Brooks/Cole, Pacific Grove, CA.

Venables W~N, Ripley B~D (2002). Modern Applied Statistics with S. Springer, New York.

See Also

summary.glmb, glmb, glmbayes-package; rglmb, rlmb, lmb; dummy.coef

Examples

set.seed(333)
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)

ps <- Prior_Setup(counts ~ outcome + treatment, family = poisson())
glmb.D93 <- glmb(
  n = 1000,
  counts ~ outcome + treatment,
  family = poisson(),
  pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma)
)
summary(glmb.D93)


myd <- dummy.coef(glmb.D93)
print(myd)

---

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, glmb. 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; glmb, 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
###############################################################################

---