---
title: "Chapter A11: Implementation Companion for Independent Normal-Gamma"
author: "Kjell Nygren"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Chapter A11: Implementation Companion for Independent Normal-Gamma}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: REFERENCES.bib
reference-section-title: References
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(glmbayes)
```

# Chapter A11: Implementation Companion for Independent Normal-Gamma

## 1. Introduction

This vignette documents the implementation of the **independent Normal-Gamma extension** for Gaussian regression, focusing on the internal workflow of:

1. `rIndepNormalGammaReg` (non-`_std`) - the top-level director of traffic for the independent Normal-Gamma sampling path; and
2. `EnvelopeOrchestrator` - which composes the coefficient-envelope construction (`EnvelopeBuild`) with the dispersion-aware refinement (`EnvelopeDispersionBuild`).

The discussion then deep-dives into the **accept-reject sampler** used after the envelope is constructed:
`rIndepNormalGammaReg_std` (serial) / `rIndepNormalGammaReg_std_parallel` (parallel).

Theory justification and the underlying accept-reject / envelope construction are described in Chapter A07 (the independent Normal-Gamma extension theory) and in [@Nygren2006] for the base likelihood subgradient approach. Implementation of independent Normal--Gamma regression is discussed in [@GriffinBrown2010].

For the fixed-dispersion coefficient-envelope implementation details, see Chapter A05.

Prior construction and the Gaussian calibration of `shape`, `rate`, `dispersion`, and coefficient covariances returned by `Prior_Setup()` are derived in **Chapter A12** (<https://knygren.r-universe.dev/articles/glmbayes/Chapter-A12.html>).

## 2. Core Function #1: `rIndepNormalGammaReg` (non-`_std`) - Internal Workflow

The entry point in C++ is `src/rIndepNormalGammaReg.cpp`, function `rIndepNormalGammaReg(...)`.
At a high level, it performs the following stages:

1. **Dispersion centering / refinement** via `EnvelopeCentering()` to obtain:
   - `dispersion2` (the dispersion anchor used for the envelope), and
   - `RSS_Post2` (the expected posterior weighted RSS used to update the Gamma precision parameters).
2. **Posterior mode optimization** via `optim()` (BFGS) to obtain:
   - posterior mode `bstar`, and
   - Hessian `A1`.
3. **Model standardization** via `glmb_Standardize_Model()` to obtain standardized inputs for envelope building:
   `bstar2`, `A`, `x2_std`, `mu2_std`, `P2_std`, and transform matrices for back-transformation.
4. **Envelope construction** via `EnvelopeOrchestrator()` (Step D):
   - coefficient envelope via `EnvelopeBuild`,
   - dispersion-aware refinement via `EnvelopeDispersionBuild`, and
   - sorting of envelope components via `EnvelopeSort`.
5. **Sampling** by delegating to the standardized accept-reject sampler:
   - `rIndepNormalGammaReg_std` (serial) or
   - `rIndepNormalGammaReg_std_parallel` (parallel),
   depending on `use_parallel` and `n`.
6. **Back-transform** to the original parameterization and assemble outputs.

### 2.1 Stage A: `EnvelopeCentering()` for `dispersion2` and `RSS_Post2`

`rIndepNormalGammaReg` calls:

```r
centering_out = EnvelopeCentering(y, x, mu, P, offset, wt, shape, rate, Gridtype, verbose)
dispersion2 = centering_out$dispersion
RSS_Post2   = centering_out$RSS_post
```

In implementation terms, `EnvelopeCentering()`:
- initializes `dispersion2` from a weighted residual variance computed from `lm.wfit`,
- iteratively refines `dispersion2` using a fixed-point style loop, and
- updates the Gamma parameters using the expected posterior weighted RSS.

The key output for the dispersion-aware refinement is `RSS_Post2`, which drives the precision-rate update inside `EnvelopeOrchestrator`.

### 2.2 Stage B: Posterior mode (`bstar`) and Hessian (`A1`)

After centering, `rIndepNormalGammaReg` performs posterior mode optimization using R's `optim()` with:
- BFGS,
- supplied functions `f2` and `f3` (negative log-posterior and its gradient),
- `hessian = TRUE` to obtain a Hessian estimate `A1`.

The mode `bstar` and curvature `A1` are then used to standardize the model into a form suitable for envelope construction.

### 2.3 Stage C: Standardization (`glmb_Standardize_Model`)

`glmb_Standardize_Model()` transforms the model so that the prior mean is zero and the posterior curvature is represented in a standardized coordinate system.

In implementation terms, the code extracts:
- `bstar2` (standardized mode),
- `A` (standardized Hessian / precision contribution),
- `x2_std`, `mu2_std`, `P2_std` (standardized design and prior components),
- `L2Inv`, `L3Inv` (matrices used later to back-transform sampled coefficients).

For the full theoretical and implementation details of standardization, see Chapter A05.

### 2.4 Stage D: `EnvelopeOrchestrator()` - envelope + dispersion refinement

`rIndepNormalGammaReg` delegates the construction of the coefficient envelope and the dispersion-aware extension to `EnvelopeOrchestrator`.

It passes:
- standardized mode and curvature objects (`bstar2`, `A`),
- standardized design/prior terms (`x2_std`, `mu2_std`, `P2_std`),
- offsets/weights (`alpha`, `wt`),
- envelope grid controls (`n`, `Gridtype`, `n_envopt`),
- Gamma prior hyperparameters (`shape`, `rate`),
- and centering outputs (`RSS_Post2`, plus bounds / truncation controls).

The orchestrator returns:
- `Env` (sorted envelope object including the PLSD component),
- `gamma_list` (tilt parameters and dispersion bounds),
- `UB_list` (the various constants used by the sampler accept-reject test),
- and additional diagnostics including `low` / `upp`.

### 2.5 Stage E: Delegation to `rIndepNormalGammaReg_std` / `_parallel`

After envelope construction, the code selects:
- serial sampler `rIndepNormalGammaReg_std(...)`, or
- parallel sampler `rIndepNormalGammaReg_std_parallel(...)`,
based on `use_parallel` and `n`.

This is where the accept-reject loop actually executes and where `iters_out` is recorded.

### 2.6 Stage F: Back-transform and return values

The standardized samples `beta_out` are mapped back to the original parameterization using:
- the inverse standardization transforms (`L2Inv`, `L3Inv`),
- and the prior mean shift (`mu`).

The returned structure mirrors the expected R-level contract:
- samples (coefficients and dispersion),
- plus runtime / iteration diagnostics (`iters_out`, `weight_out`).

## 3. Core Function #2: `EnvelopeOrchestrator` - Composition of Envelope Steps

`src/EnvelopeOrchestrator.cpp` provides a single orchestrator that composes:

1. `EnvelopeBuild()` - fixed-dispersion coefficient envelope, and
2. `EnvelopeDispersionBuild()` - dispersion-aware refinement for the independent Normal-Gamma extension,
followed by:
3. `EnvelopeSort()` on the R side.

### 3.1 Orchestrator inputs and internal Gamma anchoring

Inside `EnvelopeOrchestrator`, the Gamma posterior parameters for the dispersion precision are computed from:
- the prior hyperparameters (`shape`, `rate`),
- and the centering output `RSS_Post2`.

The implementation computes:
- `shape2 = shape + n_w/2`,
- `rate3  = rate + RSS_Post2/2`,
- and the dispersion anchor:
  `d1_star = rate3 / (shape2 - 1)`.

This anchor is used to rescale weights (`wt2 = wt / d1_star`) when building the coefficient envelope.

### 3.2 Step 1: `EnvelopeBuild()` call (coefficient-envelope)

The orchestrator calls `EnvelopeBuild(...)` in C++ with:
- standardized mode (`bstar2`),
- standardized curvature / precision object (`A`),
- standardized design/prior pieces (`x2`, `mu2`, `P2`),
- offsets (`alpha`) and weight scaling (`wt2`),
- family/link fixed to Gaussian/identity.

Implementation detail worth noting: the orchestrator uses internal policy around grid size and sorting to avoid redundant work when followed by the dispersion-aware build.

### 3.3 Step 2 (Deep Dive): `EnvelopeDispersionBuild()` - dispersion-aware extension

This subsection focuses on `EnvelopeDispersionBuild()` as the **extension-specific** computational workhorse for dispersion.

#### 3.3.1 Inputs and the role of `RSS_Post2`

`EnvelopeDispersionBuild` receives:
- the coefficient-envelope object `Env` from `EnvelopeBuild`,
- prior hyperparameters (`shape`, `rate`),
- the standardized model pieces (`P`, `y`, `x2`, `alpha`, `mu`, `wt`),
- the number of observations,
- and critically:
  - `RSS_post` = the centering output `RSS_Post2`,
  - `RSS_ML` (currently `NA` in this path).

At a conceptual level, `RSS_post` determines where the dispersion-axis Gamma proposal is anchored and how the dispersion refinement changes the proposal shape/rate parameters.

#### 3.3.2 Dispersion interval construction (`low`, `upp`)

The implementation uses:
- `max_disp_perc` (default 0.99) as a central truncation mass level,
- and optional `disp_lower` / `disp_upper` overrides,
to produce the final dispersion truncation interval:
`[low, upp]` on the dispersion scale.

This interval becomes the domain over which all envelope components must dominate.

#### 3.3.3 Anchor dispersion and face slopes

The function computes an anchor dispersion (often the mode of the surrogate Gamma distribution on the dispersion scale) and uses it to evaluate:
- face slopes / gradient-related quantities used in the dispersion-axis refinement.

This anchor is used to extrapolate dispersion-dependent face constants to both ends of the dispersion interval.

#### 3.3.4 RSS/UB minimization and construction of `UB_list`

The dispersion-aware machinery constructs the constants required by the accept-reject sampler:
- global residual-minimum quantity `RSS_Min`,
- per-face `UB2min`,
- along with the various log/linear constants and diagnostic maxima used in the acceptance bound.

These populate `UB_list`, which is later consumed by `rIndepNormalGammaReg_std` to compute:
- `UB1`, `UB2`, `UB3A`, `UB3B`,
- and the final accept/reject `test` statistic.

#### 3.3.5 Gamma tilt parameters (`gamma_list`)

Finally, `EnvelopeDispersionBuild` updates the Gamma proposal parameters through the dispersion-axis correction logic.
It returns:
- `shape3` and `rate2`-type parameters for the precision proposal,
- plus the resulting dispersion bounds (`disp_lower`, `disp_upper`),
and a collection of diagnostics fields.

#### 3.3.6 Returned object summary

The main return object is a list with:
- `Env_out` (the updated envelope with dispersion-aware PLSD / mixture structure),
- `UB_list` (accept-reject constants and UB diagnostics),
- `gamma_list` (Gamma proposal tilt params),
- and additional diagnostics (e.g., anchor value, computed maxima, UB minima).

This entire structure is consumed by the sampler in `rIndepNormalGammaReg_std`.

### 3.4 Step 3: `EnvelopeSort()` (component reordering)

`EnvelopeOrchestrator` then calls the R function `EnvelopeSort(...)` to reorder components based on probability mass in a way that improves simulation runtime.

The sorted envelope and associated UB constants are carried forward into the accept-reject simulation.

## 4. Simulation Execution: `rIndepNormalGammaReg_std` - Serial Accept-Reject Loop

After `EnvelopeOrchestrator()` returns the sorted envelope together with the dispersion proposal parameters and the accept-reject constants, `rIndepNormalGammaReg_std` generates iid posterior draws using an envelope-based rejection sampler.

Conceptually, each accepted draw consists of two pieces:

1. a regression coefficient vector `beta` proposed from an envelope-restricted multivariate normal (face-wise), and
2. a dispersion/precision draw `dispersion` (via a truncated inverse-Gamma proposal).

The accept-reject decision is then made by checking that a bound constructed from `UB_list` and envelope-derived tangency quantities dominates the target density.

### 4.1 Precomputation (`Inv_f3_precompute_disp`)

Before the per-draw loop, the sampler builds a cache once:

```r
cache = Inv_f3_precompute_disp(cbars, y, x, mu, P, alpha, wt)
```

This cache is reused inside the accept-reject loop so that repeated evaluations at dispersion-dependent tangency points avoid rebuilding expensive matrices.

### 4.2 Proposal generation (per accepted draw)

For each draw index `i`, the sampler repeats a `while` loop until the proposal is accepted:

1. **Select an envelope face (component)**
   - Draw an index `J` from the discrete distribution defined by `PLSD`.
2. **Propose regression coefficients `beta` on the restricted interval**
   - For each coordinate `j = 1,...,p`, sample from the restricted normal using the envelope bounds `loglt[J, ]` and `logrt[J, ]`.
   - The restricted normal is centered using the face-specific anchoring values `-cbars(J, j)`.
3. **Propose dispersion**
   - Draw `dispersion` from the truncated proposal using:
     `rinvgamma_ct_safe(shape3, rate2, disp_upper, disp_lower)`.
4. **Update dispersion-dependent tangency quantities**
   - With the proposed `dispersion`, recompute the dispersion-adjusted tangency point(s) by calling:
     `Inv_f3_with_disp(cache, dispersion, cbars_small)`.
   - Evaluate the log-likelihood at both:
     - the updated tangency point (stored as `LL_New2`), and
     - the candidate `beta` (stored as `LL_Test`).
     The likelihood evaluations use weights adjusted as `wt2 = wt / dispersion`.

### 4.3 Accept-reject test statistic

With the proposal `(beta, dispersion)` and the chosen face `J`, the sampler computes the acceptance bound via four non-negative surplus terms:

- `UB1`: tangent-plane style surplus using `LL_New2` and the dispersion-adjusted anchoring (`cbars` and tangency points).
- `UB2`: RSS-based surplus derived from `rss_face_at_disp(...)` and the global/per-face RSS minima `RSS_Min` and `UB2min[J]`.
- `UB3A`: face-specific surplus term involving `lg_prob_factor[J]` and the function `g1_face_at_disp(...)`.
- `UB3B`: dispersion-axis correction term involving the log-tilt components and precomputed maxima.

Implementation-wise, it forms:

```text
test1 = LL_Test - UB1
test  = test1 - (UB2 + UB3A + UB3B)
test  = test - log_U2   (where U2 ~ Uniform(0,1))
accept if test >= 0
```

Equivalently, since `log_U2 = log(U2)`:

```text
accept iff  log(U2) <= (LL_Test - UB1) - (UB2 + UB3A + UB3B).
```

The implementation also enforces sign constraints expected from a valid dominating envelope:

- `test1 <= 0`,
- `UB2 >= 0` (with extra tolerance/diagnostics when `UB2` is numerically tiny),
- `UB3A >= 0` and `UB3B >= 0`.

If any sign constraint is violated, the code throws a runtime error with context to help diagnose invalid envelope constants.

### 4.4 Record outputs

Once accepted:

The fields below are the ones from the envelope/UB structures that directly drive the sampler:
- `PLSD`: selects the face index `J`;
- `loglt` / `logrt` and `cbars`: define the restricted normal proposal for `beta` (and the face-centering);
- `shape3`, `rate2`, `disp_lower`, `disp_upper`: define the dispersion proposal;
- `RSS_Min`, `UB2min`, `lg_prob_factor`, `lmc1`, `lmc2`, and related precomputed maxima: define `UB1`-`UB3B` in the acceptance test.

- `beta_out[i, ]` stores the accepted coefficients,
- `disp_out[i]` stores the accepted dispersion draw,
- `iters_out[i]` records how many rejection attempts were needed for that accepted draw.

## 5. Dataflow: From `EnvelopeOrchestrator` Outputs to the Sampler

This section is a "developer dataflow map" connecting `EnvelopeOrchestrator` / `EnvelopeDispersionBuild` outputs to the serial accept-reject sampler.

### 5.1 What the sampler reads from `Env`

The sampler reads from `Env`:
- `PLSD` (mixture component selection weights),
- `loglt` / `logrt` (restricted normal coordinate bounds),
- `logP` / `LLconst` and related constants used internally,
- `cbars` (face gradients used to center restricted normals).

### 5.2 What the sampler reads from `gamma_list` / `UB_list`

The sampler reads from:

- `gamma_list`:
  - `shape3` / `rate2` parameters,
  - `disp_lower` / `disp_upper` truncation parameters on the dispersion proposal.

- `UB_list`:
  - `RSS_Min`,
  - `UB2min`,
  - and the global constants used for computing `UB3A` and `UB3B`,
  - plus diagnostic maxima used in these corrections.

This is the key interface that ensures the dispersion-aware envelope dominates the target on the allowed truncation domain.

## 6. Practical Knobs and Diagnostics

### 6.1 Tuning parameters

Key arguments affecting the internal workflow:
- `n` (number of iid draws) - affects how often sorting / envelope sizing is amortized,
- `Gridtype` and `n_envopt` - controls envelope grid complexity,
- `use_parallel` and `use_opencl` - controls whether the sampling and/or envelope evaluation is parallelized / accelerated,
- `max_disp_perc` and optional `disp_lower` / `disp_upper` - controls truncation interval tightness.

### 6.2 Interpreting `iters_out`

The vector `iters_out` captures the number of accept-reject attempts per accepted draw.
Large values indicate a loose envelope (or expensive bounds), while values near 1 suggest a tight envelope and efficient acceptance.

### 6.3 Minimal examples

For runnable examples that exercise the envelope construction and sampler path, see:
- Chapter A07 (independent Normal-Gamma extension theory),
- Chapter A05 (fixed-dispersion envelope build),
- and the `inst/examples/` directory for end-to-end usage patterns.

## 7. Cross-References

- Chapter A05 for coefficient-envelope construction (`EnvelopeBuild`) and standardization (`glmb_Standardize_Model`).
- Chapter A07 for the independent Normal-Gamma extension theory and the accept-reject decomposition.
- Chapter A06 for dispersion accept-reject foundations in other Gaussian families (as relevant).