Core C++ Sampling Engine for 'glmbayes'

Core C++ engine for 'glmbayes': envelope-based iid linear and generalized linear model samplers, prior-family routing, and optional 'OpenCL' acceleration. Sampling for supported non-conjugate models uses accept-reject methods based on likelihood subgradients as in Nygren and Nygren (2006) . Intended as a developer backend for the 'glmbayes' formula interface; end users should use 'glmbayes' for modelling with interfaces analogous to 'lm' and 'glm'. Mixed-model engines are planned for a future release.

glmbayesCore

glmbayesCore is the compiled sampling engine that powers the glmbayes ecosystem. It holds the C++/OpenCL envelope samplers, the family-function infrastructure, and the R-level prior and simulation interfaces that downstream packages depend on. End users should install glmbayes rather than this package directly.

The relationship to the broader ecosystem parallels how StanHeaders / rstan serve as the compiled backbone for rstanarm: glmbayesCore is the infrastructure layer; glmbayes and the in-development lmebayes are the user-facing packages built on top of it.

Current staging note. This tree currently ships the iid GLM/LM envelope engine used by glmbayes. Mixed-model (LMM/GLMM / two-block) engines are still part of the long-term glmbayesCore API and are under active development in the temporary lmebayesCore fork (consumed by lmebayes). Think of lmebayesCore as a development holding package: features return here gradually once the iid backend is stable (CRAN / glmbayes re-import path).

Package Ecosystem

Target architecture (after mixed-model reintegration):

                ┌─────────────────────────────────────────┐
                │           End-user packages             │
                │   glmbayes  ·  lmebayes  ·  (others)    │
                └──────────────────┬──────────────────────┘
                                   │ Imports / LinkingTo
                ┌──────────────────▼──────────────────────┐
                │              glmbayesCore               │
                │  iid GLM/LM · LMM/GLMM · OpenCL         │
                │  pfamily · simfunctions · rglmb/rlmb    │
                └──────────────────┬──────────────────────┘
                                   │ Imports
                ┌──────────────────▼──────────────────────┐
                │   opencltools  ·  nmathopencl            │
                │   Rcpp · RcppArmadillo · RcppParallel   │
                └─────────────────────────────────────────┘

Temporary staging (today): glmbayes → glmbayesCore (iid); lmebayes → lmebayesCore (full fork including mixed-model stack). The fork collapses back into glmbayesCore as features are merged.

glmbayes adds the formula interface (glmb(), lmb()), MCMC diagnostics, and the full suite of S3 methods that mirror base-R's lm() / glm().

lmebayes (in development) extends the engine to linear / generalized linear mixed-effects models (lmerb(), glmerb()).

What Is Inside glmbayesCore

C++ sampling engine (src/)

The core is organized under the glmbayes:: namespace:

Export wrappers in export_wrappers.cpp and kernel_wrappers.cpp expose selected entry points to R via Rcpp.

OpenCL kernels (inst/cl/)

For systems with an OpenCL-capable device, envelope construction can be offloaded to the GPU. The inst/cl/ tree contains family/link f2/f3 kernels, an OpenCL port of R Mathlib probability functions, and shim headers. Kernel loading for exploration uses opencltools; runtime GPU assembly uses kernel_loader.cpp and kernel_runners.cpp.

R-level infrastructure (R/)

Architecture: How pfamilies Route to Simulation Functions

A pfamily object is a self-contained prior specification. Every constructor bundles the hyperparameters into a prior_list and embeds a simfun function pointer. When rglmb() draws samples, it calls pfamily$simfun(y, x, prior_list, family, ...) — there is no internal switch on prior type.

rglmb(y, x, pfamily = dNormal(...), family = poisson())
          │
          └─► pfamily$simfun  ──►  rNormal_reg()
                                       │
                              family == gaussian?
                              ├── Yes ──► conjugate multivariate normal draw
                              └── No  ──► envelope sampling (Nygren & Nygren, 2006)
                                              │
                                              └──► rNormalGLM (C++)

Prior_Setup() fits an auxiliary GLM and returns calibrated hyperparameters on the same scale as the design matrix.

Architecture: How Simulation Functions Route to C++ Samplers

rNormal_reg()

rNormal_reg(y, x, prior_list, family, ...)
       │
  family$family == "gaussian"?
  ├── Yes ──► direct MVN draw via backsolve / Cholesky
  └── No  ──► EnvelopeOrchestrator (R)
                   ├── EnvelopeBuild (C++)
                   └── rNormalGLM (C++)   [accept-reject; optional OpenCL envelope]

rindepNormalGamma_reg()

rindepNormalGamma_reg(y, x, prior_list, ...)
       │
       └──► rIndepNormalGammaReg (C++)
                   ├── EnvelopeBuild_Ind_Normal_Gamma per dispersion grid point
                   └── joint accept-reject over (beta, dispersion)

rGamma_reg()

rGamma_reg(y, x, prior_list, family, ...)
       │
  family$family == "gaussian"?
  ├── Yes ──► rGammaGaussian (C++)
  └── No  ──► rGammaGamma (C++)

Architecture: How rglmb() Orchestrates a Draw

rglmb() validates the family × pfamily combination and delegates sampling to the simfun embedded in the pfamily object. In glmbayes, glmb() and lmb() wrap rglmb() / rlmb() with formula parsing.

rglmb(y, x, family = poisson(), pfamily = dNormal(mu, Sigma), n = 1000)
  │
  ├─ 1. Resolve family
  ├─ 2. Unpack pfamily (okfamilies, plinks, prior_list, simfun)
  ├─ 3. Validate combination
  ├─ 4. outlist ← simfun(...)
  └─ 5. Post-process → class c("rglmb", "glmb", "glm", "lm")

Adding a new prior family requires a new pfamily constructor and simulation function — not changes to rglmb() itself.

Function overview

Symbols below are exported from glmbayesCore today (iid path). End users typically load glmbayes (or lmebayes for mixed models). Mixed-model exports temporarily ship from lmebayesCore and will return here.

Shared with glmbayes (iid GLM / LM)

Retain as glmbayes re-exports

Phase out of glmbayes (stay in glmbayesCore)

Planned mixed-model API (temporary lmebayesCore; returns here)

These are part of the long-term glmbayesCore surface for lmebayes. They are not exported from this tree today.

Typical lmebayes workflow (via lmebayesCore for now): model_setup() → Prior_Setup_lmebayes() → pfamily_list(ps) → lmerb() / glmerb().

Developer Interface Levels

Level 1 — C++ (via LinkingTo)

#include "glmbayesCore/famfuncs.h"
#include "glmbayesCore/Envelopefuncs.h"
#include "glmbayesCore/simfuncs.h"
#include "glmbayesCore/R_interface.h"

Level 2 — R simulation functions

library(glmbayesCore)
fit <- rindepNormalGamma_reg(
  y = y, x = X, n = 2000,
  prior_list = dIndependent_Normal_Gamma(mu, Sigma, shape, rate)$prior_list,
  family = gaussian()
)

Level 3 — rglmb() / rlmb() with pfamily objects

ps  <- Prior_Setup(y, X, family = poisson())
fit <- rglmb(y = y, x = X, n = 1000,
             pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma),
             family  = poisson())

Installation

GitHub / R-Universe (recommended for developers):

install.packages("glmbayesCore",
                 repos = c("https://cloud.r-project.org",
                           "https://knygren.r-universe.dev"))

From source (required for OpenCL GPU support):

install.packages("glmbayesCore", type = "source",
                 repos = "https://knygren.r-universe.dev")

See Chapter 16 — Large models: GPU acceleration using OpenCL for system-level setup instructions.

Dependencies that must be installed first:

install.packages(c("Rcpp", "RcppArmadillo", "RcppParallel", "MASS", "Rdpack"))
install.packages(c("opencltools", "nmathopencl"),
                 repos = "https://knygren.r-universe.dev")

Extending glmbayesCore

Adding a new pfamily

See inst/ADDING_PFAMILY.md. In summary:

1. Write a constructor in pfamily.R that builds prior_list and sets simfun.
2. Implement or reuse a simulation function in simfunction.R.
3. If a new C++ sampler is needed, add it under src/, register via Rcpp::compileAttributes(), and expose it through export_wrappers.cpp.
4. For GPU support, add the corresponding f2/f3 OpenCL kernel under inst/cl/src/ and register it in kernel_loader.cpp.

Block Gibbs / mixed-model engines

The orchestrator pattern is intentionally generic: validate a model specification, unpack a routing object, call the embedded simfun, post-process. Mixed-effects drivers in lmebayes (lmerb(), glmerb()) build on two-block Gibbs engines that will return to this package. While those engines are staged in lmebayesCore, architecture notes (ergodicity, rGLMM_sweep / Block~1–Block~2 call chains, C++ migration plans) live there and will move back with the code.

Key References

• Nygren, K.N. and Nygren, L.M. (2006). Likelihood Subgradient Densities. Journal of the American Statistical Association, 101(475), 1144–1156. — The accept-reject envelope method at the heart of the non-Gaussian samplers.
• Lindley, D.V. and Smith, A.F.M. (1972). Bayes estimates for the linear model. Journal of the Royal Statistical Society B, 34, 1–41. — Conjugate Normal-Gamma foundations.
• Gelman, A. et al. (2013). Bayesian Data Analysis, 3rd ed. — Reference for prior specifications and dispersion modeling.

A complete bibliography is in inst/REFERENCES.bib.

Future plans

• Reintegrate mixed-model stack from temporary lmebayesCore: Gradually merge LMM/GLMM setup (model_setup, Prior_Setup_lmebayes, pfamily_list), matrix drivers (rlmerb / rglmerb), and two-block / block-ING engines back into glmbayesCore, then point lmebayes at this package again and retire lmebayesCore.
• Sweep-outer drivers and sweep_history on all two-block paths: Mixed-model sampling should use a sweep-outer loop (all chains complete inner sweep m, then m+1, …) on every route, for consistency with rGLMM_sweep() / rGLMM_reg_*. Each stored draw should attach sweep_history (class two_block_sweep_history) so print() and plot_sweep_history_diag() can diagnose inner-Gibbs convergence. Sweep-outer R drivers and ING pilot/main paths already capture history in the lmebayesCore fork; gaps remain on some Gaussian fixed-τ² / fixed-σ² routes that still use a chain-outer C++ driver without per-sweep cross-chain stats.
• C++ inner-chain loops and within-block parallel sampling: Per-sweep inner chain loops (Block~1 + Block~2 updates across replicate chains) should migrate from R orchestration into src/*.cpp drivers. Parallelism should be within-block, across chains at fixed inner sweep m (not parallel inner sweeps): Block~1 random-effect updates over replicate chains first (higher priority); Block~2 fixed-effect / hyperparameter updates over chains where safe (ideal follow-on). Use native threading (e.g. RcppParallel) so large n does not pay full R-loop overhead.
• OpenCL loader alignment with glmbayes: LinkingTo: opencltools, thin manifest-based loader.
• CRAN path for the iid engine: Submit the slim backend, point glmbayes at glmbayesCore, and strip duplicated backend code from glmbayes — then grow mixed-model features back into Core.

License

GPL-2. See the LICENSE file and inst/COPYRIGHTS for attribution of incorporated R Mathlib sources.