1.1 Why Envelopes?

Accept-reject sampling requires a proposal density \(q(\theta)\) and a constant \(a\) such that \(a \cdot q(\theta) \geq \pi(\theta \mid y)\) for all \(\theta\). The posterior \(\pi(\theta \mid y)\) is then sampled by: draw \(\theta \sim q\), accept with probability \(\pi(\theta \mid y) / (a \cdot q(\theta))\). The constant \(a\) controls efficiency: smaller \(a\) means fewer candidates per acceptance. The likelihood subgradient approach constructs tight envelopes, with \(a\) bounded above by \((2/\sqrt{\pi})^k\) for \(k\)-dimensional models when the likelihood is approximately normal (Nygren & Nygren 2006, Theorems 2 and 3). This ensures efficient iid sampling even in high dimensions.

1.2 High-Level Flow

The envelope-based pipeline follows this sequence:

1. Optimize - Find the posterior mode.
2. Standardize - Reparameterize into standard form (prior precision = identity, posterior precision diagonal at the mode).
3. Size - Choose the number of tangency points per dimension (EnvelopeSize, EnvelopeOpt).
4. Build - Construct the envelope (EnvelopeBuild).
5. Sort - Order envelope components for sampling efficiency (EnvelopeSort).
6. Sample - Generate candidates from the mixture envelope and accept/reject (rNormalGLM_std or C++ equivalents).

For models with unknown dispersion, an additional dispersion envelope step (EnvelopeDispersionBuild) extends the coefficient envelope to a joint envelope over \((\beta, \phi)\). This is orchestrated by EnvelopeOrchestrator for the independent Normal-Gamma case.

2.1 Posterior and Definition 2 (paper Section 2)

Let the posterior be \[ \pi(\theta \mid y) \propto f(y \mid \theta)\,\pi(\theta), \] where \(f(y \mid \theta)\) is the likelihood and \(\pi(\theta)\) is the prior density.

Definition 2. A density \(q(\cdot)\) is a generalized likelihood-subgradient density for \(\pi(\cdot\mid y)\) with prior \(\pi(\cdot)\) and likelihood \(f(y\mid\cdot)\) at \(\bar{\theta}\in\Theta\) if there exist a function \(g\) and a subgradient \(c(\bar{\theta})\) of \(-\log g\) at \(\bar{\theta}\) such that:

1. (a) \(g(\cdot)\) bounds \(f(y\mid\cdot)\) from above;
2. (b) \(\mathrm{MGF}(-c(\bar{\theta})) = \int_{\Theta} \exp(-c(\bar{\theta})^T\theta)\,\pi(\theta)\,d\theta\) is finite;
3. (c) for all \(\theta\in\Theta\), \[ q(\theta) = \frac{\exp\!\big(-c(\bar{\theta})^T\theta\big)\,\pi(\theta)} {\mathrm{MGF}\!\big(-c(\bar{\theta})\big)}. \]

Special cases (paper). If \(g\) is the likelihood, \(q\) is a likelihood-subgradient density; log-concave likelihoods (e.g., Poisson and logistic regression) are the main examples. If \(g\) is a constant upper bound on the likelihood, one obtains the Bayesian construction where the prior itself is viewed as a generalized likelihood-subgradient density.

Existence (paper Appendix). Likelihood-subgradient densities exist at every \(\bar{\theta}\) when the prior is a finite mixture of multivariate normals and the likelihood is log-concave. Mixture priors yield particularly simple algebraic forms (Claim 1 below).

2.2 Claim 1: finite mixture of normal priors

Suppose \(\pi(\theta) = \sum_{i=1}^{k} p_i\,\pi_i(\theta\mid \mu_i,\Sigma_i)\) and \(c(\bar{\theta})\) is a subgradient for \(-\log g\) at \(\bar{\theta}\). Then \[ \mathrm{MGF}\!\big(-c(\bar{\theta})\big) = \sum_{i=1}^{k} p_i\, \exp\!\Big( -c(\bar{\theta})^T\mu_i + \tfrac12\,c(\bar{\theta})^T\Sigma_i\,c(\bar{\theta}) \Big), \] and the generalized likelihood-subgradient density is again a mixture of multivariate normals: \[ q(\theta) = \sum_{i=1}^{k} \tilde{p}_i\,\pi_i(\theta\mid \tilde{\mu}_i,\Sigma_i), \] with \[ \tilde{p}_i = \frac{ p_i\,\exp\!\big(-c(\bar{\theta})^T\mu_i + \tfrac12\,c(\bar{\theta})^T\Sigma_i\,c(\bar{\theta})\big) }{ \mathrm{MGF}\!\big(-c(\bar{\theta})\big) }, \qquad \tilde{\mu}_i = \mu_i - \Sigma_i\,c(\bar{\theta}). \]

2.3 Theorem 1: envelope dominance and tangency

Define \[ a(\bar{\theta}) = \frac{ g(\bar{\theta})\,\mathrm{MGF}\!\big(-c(\bar{\theta})\big) }{ f(y)\,\exp\!\big(-c(\bar{\theta})^T\bar{\theta}\big) }, \] and \[ h_{\bar{\theta}}(\theta) = \frac{ \exp\!\big(-c(\bar{\theta})^T\bar{\theta}\big)\,f(y\mid\theta) }{ \exp\!\big(-c(\bar{\theta})^T\theta\big)\,g(\bar{\theta}) }. \]

Then Theorem 1 states \[ a(\bar{\theta})\,q_{\bar{\theta}}(\theta) \ge a(\bar{\theta})\,h_{\bar{\theta}}(\theta)\,q_{\bar{\theta}}(\theta) = \pi(\theta\mid y), \] with \[ 0 \le h_{\bar{\theta}}(\theta) \le 1 \quad \forall \theta\in\Theta. \]

If \(f(y\mid\bar{\theta}) = g(\bar{\theta})\), then \(h_{\bar{\theta}}(\bar{\theta})=1\), so the envelope touches the posterior at tangency.

Operationally, \(a(\bar{\theta})\) governs acceptance efficiency: the acceptance probability is \(1/a(\bar{\theta})\).

2.4 Models in standard form

Following Nygren & Nygren (2006, Section 3.3), the envelope construction assumes the model has been reparameterized into standard form. In this form:

• The prior precision matrix is the identity.
• The posterior precision matrix, evaluated at the posterior mode \(\theta^\star\), is diagonal.

Example 2 in the paper illustrates the special case of a zero-mean normal prior with identity covariance. In this setting:

• The generalized likelihood-subgradient density at a tangency point \(\bar{\theta}\) has mean vector \(-c(\bar{\theta})\) and covariance \(I\).

• The normalizing integrals over restricted sets factorize across dimensions. Specifically, for a rectangular set \(A = \{\theta : l_L \le \theta \le l_U\}\), we have

  \[ \int_{\theta \in A} q^{\bar{\theta}}(\theta)\, d\theta \;=\; \prod_{r=1}^{p} \Big[ \Phi\!\big(l_{U,r} + c_{r}(\bar{\theta})\big) - \Phi\!\big(l_{L,r} + c_{r}(\bar{\theta})\big) \Big], \]

  and the truncated expectation in coordinate \(r\) is

  \[ \mathbb{E}_{\tilde{q}^{\bar{\theta}}}[\theta_r \mid \theta \in A] = -\,c_{r}(\bar{\theta}) + \frac{ \phi\!\big(l_{L,r} + c_{r}(\bar{\theta})\big) - \phi\!\big(l_{U,r} + c_{r}(\bar{\theta})\big) }{ \Phi\!\big(l_{U,r} + c_{r}(\bar{\theta})\big) - \Phi\!\big(l_{L,r} + c_{r}(\bar{\theta})\big) }. \]

These closed-form expressions explain why EnvelopeBuild evaluates logU, loglt, and logrt using univariate normal CDFs and densities, rather than numerical integration. The gradients (cbars) directly determine the shifted means of the restricted densities, and the separability across dimensions makes the grid-based construction computationally tractable.

Models with Zellner’s \(g\)-priors are essentially in standard form, since in the whitened design space both the prior and likelihood precisions are diagonal. Each dimension still needs to be scaled so that the prior precision is exactly the identity matrix. For other models, standard form can be achieved by reparameterization (e.g. via Cholesky of the posterior precision) or by shifting part of the prior quadratic form into the likelihood.

2.5 Construction of restricted subgradient densities

Following Remark 5 in Nygren & Nygren (2006), each unrestricted likelihood-subgradient density \(q_{\bar{\theta}}(\cdot)\) can be restricted to a subset \(A \subset \Theta\). The restricted density is defined as

\[ \tilde{q}_{\bar{\theta}}(\theta) = \frac{ q_{\bar{\theta}}(\theta)\,\mathbf{1}_{\{\theta \in A\}} } { \int_{\theta' \in A} q_{\bar{\theta}}(\theta')\, d\theta' }, \] and the corresponding constant is

\[ \tilde{a}(\theta_{\text{bar}}) = a(\theta_{\text{bar}})\, \int_{\theta \in A} q_{\theta_{\text{bar}}}(\theta)\, d\theta, \] where \(a(\theta_{\text{bar}})\) is the global normalizing constant from Theorem 1. For every \(\theta \in A\), the identity

\[ \tilde{a}(\theta_{\text{bar}})\, h_{\theta_{\text{bar}}}(\theta)\, \tilde{q}_{\theta_{\text{bar}}}(\theta) = \pi(\theta \mid y) \] holds, ensuring that the restricted densities reproduce the posterior when combined.

The constant \(a(\bar{\theta})\) is defined in Theorem 1 as

\[ a(\bar{\theta}) = \frac{g(\bar{\theta})\,\mathrm{MGF}\!\big(-c(\bar{\theta})\big)} {f(y)\,\exp\!\big(-c(\bar{\theta})^{T}\bar{\theta}\big)}, \] where \(g(\bar{\theta})\) is the reference density at the tangency point, \(c(\bar{\theta})\) is the subgradient of the log-likelihood, and \(\mathrm{MGF}(-c(\bar{\theta}))\) is the moment-generating function of the prior evaluated at \(-c(\bar{\theta})\).

The envelope function \(h_{\bar{\theta}}(\theta)\) is defined as

\[ h_{\bar{\theta}}(\theta) = \frac{\exp\!\big(-c(\bar{\theta})^{T}\bar{\theta}\big)\,f(y\mid\theta)} {\exp\!\big(-c(\bar{\theta})^{T}\theta\big)\,g(\bar{\theta})}, \] and satisfies

\[ 0 \le h_{\bar{\theta}}(\theta) \le 1 \quad \forall\,\theta \in \Theta, \qquad h_{\bar{\theta}}(\bar{\theta}) = 1 \quad \text{if } f(y\mid\bar{\theta}) = g(\bar{\theta}). \]

In the standardized model (zero-mean normal prior with identity covariance and diagonal posterior precision at the mode), the restricted integral \(\int_{A} q_{\bar{\theta}}(\theta)\, d\theta\) factorizes across dimensions and can be evaluated in closed form using normal CDFs:

\[ \int_{\theta \in A} q_{\bar{\theta}}(\theta)\, d\theta = \prod_{r=1}^{p} \Big[ \Phi\!\big(l_{U,r} + c_{r}(\bar{\theta})\big) - \Phi\!\big(l_{L,r} + c_{r}(\bar{\theta})\big) \Big]. \]

In standardized models, where the prior is \(\mathcal{N}(0, I)\), the moment-generating function simplifies to

\[ \mathrm{MGF}\big(-c(\bar{\theta})\big) = \exp\left(\tfrac{1}{2} c(\bar{\theta})^{T} c(\bar{\theta})\right), \] so the constant \(a(\bar{\theta})\) from Theorem 1 becomes

\[ a(\bar{\theta}) = \frac{g(\bar{\theta})\,\exp\left(\tfrac{1}{2} c(\bar{\theta})^{T} c(\bar{\theta})\right)} {f(y)\,\exp\!\big(-c(\bar{\theta})^{T} \bar{\theta}\big)}. \]

This is why EnvelopeBuild computes and stores logU, loglt, and logrt using univariate normal CDF evaluations. The constants \(\tilde{a}(\bar{\theta})\) are then obtained by scaling the global constant with these integrals, and the mixture weights (PLSD) are normalized accordingly. In practice:

• EnvelopeSet_Grid_C2_pointwise evaluates restricted densities at each grid point.
• LLconst stores the log of the restricted integrals.
• EnvelopeSet_LogP_C2 computes \(\tilde{a}(\bar{\theta})\) and normalizes mixture weights.

2.6 Mixture construction and tractable probabilities

Claim 2 in Nygren & Nygren (2006) shows that the posterior density \(\pi(\theta \mid y)\) can be expressed as a mixture of restricted likelihood-subgradient densities. Let \(A_1, \dots, A_m\) be a partition of the parameter space \(\Theta\), and define

\[ \tilde{q}_{\bar{\theta}}(\theta) = \sum_{i=1}^{m} \tilde{p}_i \, q^{\bar{\theta}}_{A_i}(\theta), \qquad \tilde{p}_i = \frac{ \tilde{a}_i }{ \sum_{j=1}^{k} \tilde{a}_j }. \]

In Remark 6 of Nygren & Nygren (2006), the mixture weights \(\tilde{p}_i\) for each restricted likelihood-subgradient density \(q^{\bar{\theta}_i}_{A_i}\) are defined as

\[ \tilde{p}_i = \frac{ g(\bar{\theta}_i)\,\mathrm{MGF}\!\left(-c(\bar{\theta}_i)\right)\, \dfrac{ \int_{\theta \in A_i} q^{{\bar{\theta}_i}}(\theta)\, d\theta }{ \exp\!\left(-c(\bar{\theta}_i)^{\top}\bar{\theta}_i\right) } }{ \displaystyle\sum_{j=1}^{k} g(\bar{\theta}_j)\,\mathrm{MGF}\!\left(-c(\bar{\theta}_j)\right)\, \dfrac{ \int_{\theta \in A_j} q^{{\bar{\theta}_j}}(\theta)\, d\theta }{ \exp\!\left(-c(\bar{\theta}_j)^{\top}\bar{\theta}_j\right) } }. \]

This expression reflects the full normalization of the mixture, where each \(\tilde{p}_i\) is proportional to the restricted constant \(\tilde{a}_i(\bar{\theta}_i)\) from Remark 5, and the denominator sums over all such constants across the partition. The resulting mixture \(\tilde{q}^{\bar{\theta}}(\theta)\) is a valid approximation to the posterior density \(\pi(\theta \mid y)\).

Remark 6 emphasizes that these mixture weights are tractable to compute in standardized models. When the prior is \(\mathcal{N}(0, I)\) and the posterior precision is diagonal at the mode, each integral \(\int_{A_i} q^{\bar{\theta}}(\theta)\, d\theta\) factorizes across dimensions and can be evaluated using normal CDFs:

\[ \int_{\theta\in A_i} q^{{\bar{\theta}}}(\theta)\, d\theta = \displaystyle\prod_{r=1}^{p} \Bigl[ \Phi\!\bigl(l_{U,r}^{(i)} + c_{r}(\bar{\theta})\bigr) - \Phi\!\bigl(l_{L,r}^{(i)} + c_{r}(\bar{\theta})\bigr) \Bigr]. \]

where \(l_{L}^{(i)}\) and \(l_{U}^{(i)}\) are the bounds defining region \(A_i\).

This tractability is central to the envelope construction. It allows EnvelopeBuild to compute:

• LLconst, which stores the log of each restricted integral \(\log \int_{A_i} q^{\bar{\theta}}(\theta)\, d\theta\).
• PLSD, which stores the normalized mixture weights \(\tilde{p}_i\).

These quantities are used to construct the combined density \(\tilde{q}_{\bar{\theta}}(\theta)\) and to evaluate the envelope approximation to the posterior. Because all components are normalized and tractable, the mixture is both valid and computationally efficient.

2.7 Log-scale properties of the envelope function

The envelope function \(h_{\bar{\theta}}(\theta)\) is defined in Theorem 1 as

\[ h_{\bar{\theta}}(\theta) = \frac{\exp\!\big(-c(\bar{\theta})^{T} \bar{\theta}\big)\,f(y \mid \theta)} {\exp\!\big(-c(\bar{\theta})^{T} \theta\big)\,g(\bar{\theta})}. \]

When \(g(\bar{\theta}) = f(y \mid \bar{\theta})\), this simplifies to

\[ h_{\bar{\theta}}(\theta) = \exp\!\big( c(\bar{\theta})^{T}(\theta - \bar{\theta}) \big) \cdot \frac{f(y \mid \theta)}{f(y \mid \bar{\theta})}. \]

Taking logarithms yields

\[ \log h_{\bar{\theta}}(\theta) = c(\bar{\theta})^{T}(\theta - \bar{\theta}) + \log f(y \mid \theta) - \log f(y \mid \bar{\theta}), \]

which is tractable as long as the log-likelihood \(\log f(y \mid \theta)\) is. In particular, if the log-likelihood is concave or piecewise affine, then \(\log h_{\bar{\theta}}(\theta)\) inherits that structure and can be efficiently evaluated across grid regions.

This tractability is central to the envelope construction: it allows EnvelopeBuild to evaluate the envelope function pointwise using EnvelopeSet_Grid_C2_pointwise, and ensures that the resulting approximation remains bounded between 0 and 1. At the tangency point, we recover

\[ h_{\bar{\theta}}(\bar{\theta}) = 1, \] confirming that the envelope touches the posterior density exactly.

The key inequality that ensures envelope dominance follows from the subgradient inequality for concave functions. If \(\log f(y \mid \theta)\) is concave and \(c(\bar{\theta})\) is a subgradient at \(\bar{\theta}\), then

\[ \log f(y \mid \theta) \le \log f(y \mid \bar{\theta}) + c(\bar{\theta})^{T}(\theta - \bar{\theta}), \]

which implies

\[ \log h_{\bar{\theta}}(\theta) \le 0, \qquad h_{\bar{\theta}}(\theta) \le 1. \]

This inequality guarantees that the envelope function dominates the posterior density pointwise, as required by Theorem 1. Equality holds at the tangency point \(\theta = \bar{\theta}\), where

\[ h_{\bar{\theta}}(\bar{\theta}) = 1. \]

2.8 Standard form: Definition 3 and Remarks 11-15 (paper Section 3.3)

Remark 11. If the log-likelihood is concave and twice continuously differentiable, a Cholesky reparameterization from the posterior precision at the unique mode achieves diagonal posterior Hessian at the mode.

Remark 12. For any positive definite prior precision \(P\), there exists diagonal \(D\succ 0\) such that \(P-D\succ 0\).

Remark 13. With \(P\) and \(D\) as in Remark 12, shifting \((P-D)\) from the prior quadratic form into the likelihood leaves the posterior unchanged: the same \(\pi(\theta\mid y)\) arises from prior precision \(D\) and an augmented log-likelihood.

Remark 14. If the prior covariance is diagonal and the posterior Hessian at the mode is diagonal, the model can be reparameterized into standard form.

Remark 15. Any model with Normal prior and log-concave, twice differentiable log-likelihood can be reparameterized into standard form.

2.10 Remarks 7-8: sampling restricted normals (paper)

Remark 7. Sampling from restricted Normal laws can be done by inverse transform (Fishman 1999).

Remark 8. When intervals lie in extreme Gaussian tails, accurate evaluation of normal CDFs (or their logs) with uniformly small relative error matters; standard approximation theory and tables are discussed in (Hart et al. 1968); normal distribution properties relevant to tail behavior appear in (Bryc 2002).

2.11 Theorem 2: univariate Normal data, three intervals, and \(\tilde{a}\to 2/\sqrt{\pi}\)

Empirically (paper), a single optimally placed likelihood-subgradient component, and two-interval partitions, can become loose as sample size grows; the three-interval construction does not exhibit the same deterioration.

Let \(\theta^\ast\) be the posterior mode. Define (univariate case) \[ \omega = \frac{ \sqrt{2} - \exp\!\big(-1.20491 - 0.7321\,(0.5 - \partial^2\log f(\theta^\ast\mid y)/\partial\theta^2)\big) }{ 1 - \partial^2\log f(\theta^\ast\mid y)/\partial\theta^2 }, \] \[ \ell_1=\theta^\ast-0.5\,\omega,\quad \ell_2=\theta^\ast+0.5\,\omega, \] regions \(A_1=(-\infty,\ell_1)\), \(A_2=[\ell_1,\ell_2]\), \(A_3=(\ell_2,\infty)\), and tangencies \(\bar{\theta}_1=\theta^\ast-\omega\), \(\bar{\theta}_2=\theta^\ast\), \(\bar{\theta}_3=\theta^\ast+\omega\).

For \(N\) i.i.d. Normal\((\theta,1)\) observations, \(-\partial^2\log f(\theta^\ast\mid y)/\partial\theta^2 = N\).

Theorem 2. With \(\tilde{a}^\ast(N)\) the combined constant at sample size \(N\), \[ \lim_{N\to\infty} \tilde{a}^\ast(N) = \frac{2}{\sqrt{\pi}}. \] The paper gives a proof sketch using symmetry \(\tilde{a}_1(\bar{\theta}_1)=\tilde{a}_3(\bar{\theta}_3)\) and appendix claims.

2.12 Multivariate partition in standard form (paper; Remark 16)

Let \(\theta^\ast\) be the unique posterior mode. For each coordinate \(i\), define \[ \omega_i = \frac{ \sqrt{2} - \exp\!\big(-1.20491 - 0.7321\,(0.5 - \partial^2\log f(\theta^\ast\mid y)/\partial\theta_i^2)\big) }{ 1 - \partial^2\log f(\theta^\ast\mid y)/\partial\theta_i^2 }, \] then \(\ell_{i,1}=\theta^\ast_i-0.5\,\omega_i\), \(\ell_{i,2}=\theta^\ast_i+0.5\,\omega_i\), and \[ A_{i,1}=(-\infty,\ell_{i,1}),\quad A_{i,2}=[\ell_{i,1},\ell_{i,2}],\quad A_{i,3}=(\ell_{i,2},\infty). \] Let \(j=(j_1,\ldots,j_p)\in\{1,2,3\}^p\) and \(A^\ast_j=\prod_{i=1}^p A_{i,j_i}\). The \(A^\ast_j\) partition \(\mathbb{R}^p\).

With \(C_{j1}=\{i:j_i=1\}\), \(C_{j2}=\{i:j_i=2\}\), \(C_{j3}=\{i:j_i=3\}\), tangency \(\bar{\theta}_j\) has coordinates \(\theta^\ast_i-\omega_i\), \(\theta^\ast_i\), or \(\theta^\ast_i+\omega_i\) according as \(i\) lies in \(C_{j1}\), \(C_{j2}\), or \(C_{j3}\).

Remark 16. The resulting mixture generalized likelihood-subgradient density is a mixture of restricted multivariate Normal laws for which a direct sampling procedure exists.

2.13 Theorem 3 (paper)

For \(k\)-dimensional Normal regression in standard form, the three-interval construction yields \[ \tilde{a}\le \left(\frac{2}{\sqrt{\pi}}\right)^k. \]

For other log-concave models, posterior asymptotic normality suggests analogous practical tightness as the posterior concentrates.

3.1 Where Envelopes Are Used

For conjugate cases (e.g., Gaussian likelihood with Normal or Normal-Gamma prior), no envelope is constructed; direct sampling is used.

3.2 Pipeline Order

The following functions are invoked in sequence during envelope-based sampling:

1. glmb_Standardize_Model - Transforms the model into standard form. Required before EnvelopeBuild.

2. EnvelopeSize - Determines the grid structure (number of tangency points per dimension). Calls EnvelopeOpt when Gridtype = 2.

3. EnvelopeBuild - Constructs the envelope: tangency points, gradients, mixture weights, and sampling constants. This is the main envelope construction function.

4. EnvelopeSort - Reorders envelope components by probability mass to improve sampling efficiency.

5. EnvelopeSetGrid, EnvelopeSetLogP - Low-level helpers that compute log-densities and constants for the envelope. Typically called internally by EnvelopeBuild or the C++ routines.

6. rNormalGLM_std - Samples from the posterior using the constructed envelope via accept-reject. Used for non-Gaussian GLMs.

7. EnvelopeDispersionBuild - For Gaussian regression with independent Normal-Gamma priors, refines the coefficient envelope to a joint envelope over \((\beta, \phi)\). Called by EnvelopeOrchestrator.

8. EnvelopeOrchestrator - High-level coordinator that runs EnvelopeBuild, EnvelopeDispersionBuild, and EnvelopeSort for the independent Normal-Gamma case. Used inside the C++ sampler.