1.1. Epistemic Posture and the Critical Audit of v07 §5

The extension is performed under three epistemic constraints documented in Charter §2.5 (Sub-phase 8.6):

1. Extension, not rewriting. v07’s Theorem 7A scalar statement and proof remain canonical in their domain (\(p = 1\), \(K = 1\)). The multivariate extension is a strictly stronger statement requiring strictly stronger hypotheses (information matrix positive definite rather than non-zero; weakly informative prior in \(\mathbb{R}^p\) rather than \(\mathbb{R}\); bounded \(C^2\) regularity of the log-marginal in a ball of \(\mathbb{R}^p\) rather than in an interval of \(\mathbb{R}\)).

2. Pedagogical hierarchy. Readers entering the EB-vs-FB comparison for the first time are referred to v07 (lower technical entry barrier; less notational overhead; identical conceptual content). This vignette is the technical extension for users whose model lives in \(p > 1\) and/or \(K > 1\).

3. Critical audit of v07 §5. Per Charter §2.5, the redaction of the multivariate Theorem 7A* is preceded by an explicit line-by-line audit of the scalar proof of v07 §5 intended to detect any implicit use of a property valid in \(\mathbb{R}\) but false in \(\mathbb{R}^p\) (e.g., real monotonicity; one-dimensional change-of-variables; integration by parts on \(\mathbb{R}\); sign of one-dimensional Fisher information taken as positivity rather than positive definiteness). The audit, summarized in §4.5 below, detected no error: v07 §5 transfers cleanly to \(\mathbb{R}^p\) under the strengthened hypotheses listed in §3.

If the audit had detected an implicit error, the redaction of this vignette would have been suspended and a parallel sub-sub-sub-phase 8.6.A.1 would have opened to reopen v07 explicitly. The fact that no error was detected confirms that v07’s “reference and argument” delegation to Petrone–Rousseau–Scricciolo (2014) and Rousseau–Szabo (2017) is faithful to those references: both treat hyperparameters on Polish spaces of which finite-dimensional \(\mathbb{R}^p\) is the routine specialization.

1.2. Reading Order

For a reader already familiar with v07: §§2–3 (notation and hypotheses), §4 (Theorem 7A*), §5 (Proposition 7B*), §§6–7 (Theorem 7C* and Proposition 7D*), §8 (recommendation table for the multivariate case), §9 (open questions including the numerical-anti-fragility question O5*-EBFB specific to multivariate Laplace), §10 (connections to the operational vignette vop07 and to the API gdpar_eb()).

For a reader unfamiliar with v07: start with v07 first; return here once the scalar single-slot case is conceptually clear.

2.1. Multivariate Population Reference and Heterogeneous Slots

Let \(\eta = (\theta_{\mathrm{ref}}, a, b, W)\) as in Block 4 §2.1, with the partition into the upper-level parameter \(\theta_{\mathrm{ref}}\) and the lower-level components \(\xi = (a, b, W)\) unchanged. The multivariate extension introduces three structural dimensions:

• \(p := \dim(\theta_{\mathrm{ref}}) \geq 1\) is the dimensionality of the population reference. Under \(p = 1\), \(\theta_{\mathrm{ref}} \in \mathbb{R}\) as in v07; under \(p \geq 2\), \(\theta_{\mathrm{ref}} \in \mathbb{R}^p\) is a vector whose components anchor different aspects of the hierarchy jointly.
• \(K \geq 1\) is the number of distributional slots that the family declares. For a single-slot family (\(K = 1\), e.g., Gaussian or Poisson), \(\theta_{\mathrm{ref}}\) anchors that one slot. For a multi-slot family (\(K \geq 2\), e.g., Zero-Inflated Poisson with \(K = 2\) slots indexing \((\mu, \pi)\); Tweedie with \(K = 3\) slots indexing \((\mu, \phi, p_{\mathrm{Tweedie}})\); heterogeneous Beta+Gamma with \(K = 2\) slots), \(\theta_{\mathrm{ref}}\) anchors each slot via a slot-indexed sub-parameter \(\theta_{\mathrm{ref}}^{(k)} \in \mathbb{R}^{p_k}\) with \(\sum_{k=1}^K p_k = p\) (additively in the no-overlap case). The framework’s canonical implementation (Sub-phase 8.6.B–D of Charter §3) takes the slot-indexed case as the only operational one in the current release. A theoretically natural alternative —a shared \(\theta_{\mathrm{ref}}\) vector with cross-slot coupling across the \(K\) slots— exists in the multivariate AMM formulation but is scoped for a future block; it is not exposed as a configuration in gdpar_eb() or anywhere else in gdpar 0.0.0.9001.
• \(J \geq 1\) is the number of groups when the anchored hierarchy admits a grouping covariate (e.g., \(J\) Argentine provinces; \(J\) ecological sites). For \(J = 1\) the hierarchy reduces to the single-group case of Blocks 1–6.

The full upper-level parameter tensor of an AMM specification is therefore \(\theta_{\mathrm{ref}} \in \mathbb{R}^{J \times K \times p}\), with \((p, K, J) = (1, 1, 1)\) recovering the simplest setting and \((p, K, J)\) all \(\geq 2\) requiring the full machinery of this vignette.

2.2. Marginal Likelihood in \(\mathbb{R}^p\) and per-Slot Tensor under \(K > 1\)

The marginal likelihood of \(\theta_{\mathrm{ref}}\) integrating out the lower-level parameters under their prior \(\pi_\xi\) is, as in v07, \[L_n^{\mathrm{marg}}(\theta_{\mathrm{ref}}) \;:=\; \int p(Y_{1:n} \mid \theta_{\mathrm{ref}}, \xi)\, \pi_\xi(\xi)\, d\xi,\] now read as a function \(L_n^{\mathrm{marg}}: \mathbb{R}^{J \times K \times p} \to \mathbb{R}_{\geq 0}\).

The Empirical Bayes estimator of \(\theta_{\mathrm{ref}}\) is \[\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}} \;=\; \mathop{\arg\max}_{\theta_{\mathrm{ref}} \in \Theta}\; L_n^{\mathrm{marg}}(\theta_{\mathrm{ref}}),\] with \(\Theta \subseteq \mathbb{R}^{J \times K \times p}\) the admissible domain (typically the Cartesian product of per-component intervals or all of \(\mathbb{R}^{J \times K \times p}\)).

The Fisher information matrix of the marginal likelihood at the population truth \(\theta_{\mathrm{ref}}^*\) is \[I_{\theta\theta}^{\mathrm{marg}}(\theta_{\mathrm{ref}}^*) \;:=\; \mathbb{E}_{\eta_*}\!\left[ - \nabla^2_{\theta_{\mathrm{ref}}} \log L_n^{\mathrm{marg}}(\theta_{\mathrm{ref}})\big|_{\theta_{\mathrm{ref}}^*} \right] \;\in\; \mathbb{R}^{p \times p}.\] Under \(p = 1\), \(I_{\theta\theta}^{\mathrm{marg}}\) is a positive scalar (the Fisher information of v07 §4). Under \(p \geq 2\), \(I_{\theta\theta}^{\mathrm{marg}}\) is a \(p \times p\) symmetric matrix; the strictly stronger identifiability assumption (EB-MARG-ID)\(_p\) of §3 requires it to be positive definite (all eigenvalues strictly positive), not merely non-zero.

Under \(K > 1\) with slot-indexed sub-parameters \(\theta_{\mathrm{ref}}^{(k)} \in \mathbb{R}^{p_k}\), the Fisher information becomes a per-slot block tensor of marginal information matrices: \[I_{\theta\theta}^{\mathrm{marg}, K}(\theta_{\mathrm{ref}}^*) \;=\; \left\{ I_{\theta\theta}^{\mathrm{marg}, k}(\theta_{\mathrm{ref}}^{*,(k)}) \right\}_{k=1}^K, \qquad I_{\theta\theta}^{\mathrm{marg}, k} \;\in\; \mathbb{R}^{p_k \times p_k}.\] Under the factored prior \(\pi_\xi = \prod_{k=1}^K \pi_{\xi^{(k)}}\) of Theorem 7C* (§6), the block tensor is block-diagonal and each block \(I_{\theta\theta}^{\mathrm{marg}, k}\) is the Fisher information of the slot-\(k\) marginal likelihood. Under the coupled prior of §6, off-block cross terms appear and the full Fisher matrix \(I_{\theta\theta}^{\mathrm{marg}} \in \mathbb{R}^{p \times p}\) is not block-diagonal.

2.3. Sensitivity Jacobian of the Conditional Posterior

A novel object that has no counterpart in v07 (because v07’s scalar Proposition 7B absorbs it into a constant) is the sensitivity Jacobian of the lower-level conditional posterior mode with respect to the upper-level reference: \[J^\xi(\theta_{\mathrm{ref}}) \;:=\; \frac{\partial\, \xi^*(\theta_{\mathrm{ref}})}{\partial\, \theta_{\mathrm{ref}}} \;\in\; \mathbb{R}^{\dim(\xi) \times p},\] where \(\xi^*(\theta_{\mathrm{ref}}) := \arg\max_\xi \Pi_n^{\mathrm{cond}}(\xi \mid \theta_{\mathrm{ref}}, \cdot)\) is the lower-level conditional posterior mode given \(\theta_{\mathrm{ref}}\). By the implicit-function theorem applied to the conditional first-order condition \(\nabla_\xi \log p(Y_{1:n} \mid \theta_{\mathrm{ref}}, \xi^*) + \nabla_\xi \log \pi_\xi(\xi^*) = 0\), \[J^\xi(\theta_{\mathrm{ref}}) \;=\; -\left[ \nabla^2_{\xi\xi} \log \Pi_n^{\mathrm{cond}}(\xi^* \mid \theta_{\mathrm{ref}}, \cdot) \right]^{-1} \cdot \nabla^2_{\xi\theta_{\mathrm{ref}}} \log p(Y_{1:n} \mid \theta_{\mathrm{ref}}, \xi^*),\] which is well-defined under (HIER-COMPLEX)\(_p\) of §3. Under \(p = \dim(\xi) = 1\), \(J^\xi\) is a scalar absorbed into v07’s constant \(C_{g,\alpha}\) without explicit naming. Under \(p \geq 2\), the matrix \(J^\xi\) enters explicitly in the sandwich form of Proposition 7B* (§5).

2.4. Total Variation, Weak Metric, and Posterior Distance

The asymptotic comparison of \(\Pi_n^{\mathrm{EB}}\) and \(\Pi_n^{\mathrm{FB}}\) over \(\xi\) uses two metrics, paralleling v07 §2.2:

• Total variation \(d_{\mathrm{TV}}(\Pi_n^{\mathrm{EB}}, \Pi_n^{\mathrm{FB}}) = \sup_B |\Pi_n^{\mathrm{EB}}(\xi \in B) - \Pi_n^{\mathrm{FB}}(\xi \in B)|\) for finite-dim parametric \(\xi\) (Regime A of v07 §5).
• Weak metric over smooth functionals \(\sup_{g \in \mathcal{G}} |\mathbb{E}_{\Pi_n^{\mathrm{EB}}}[g(\xi)] - \mathbb{E}_{\Pi_n^{\mathrm{FB}}}[g(\xi)]|\) over a class \(\mathcal{G}\) of bounded \(L^2(\mu)\)-Lipschitz functionals for non-parametric \(\xi\) (Regime B of v07 §5).

The dimensionality of \(\theta_{\mathrm{ref}}\) (\(p\)) and the dimensionality of \(\xi\) are independent: the multivariate extension of this vignette concerns \(p\); the parametric-vs-non-parametric distinction of \(\xi\) governs the choice of metric. The four resulting combinations are tabulated in §4.

2.5. Notation Summary

3.1. (EB-MARG-ID)\(_p\) – Multivariate Marginal Identifiability

(EB-MARG-ID)\(_p\). \(L_n^{\mathrm{marg}}(\theta_{\mathrm{ref}})\) has a unique global maximum on the admissible domain \(\Theta \subseteq \mathbb{R}^p\) in the limit \(n \to \infty\), with the limiting Fisher information matrix \[I_{\theta\theta}^{\mathrm{marg}}(\theta_{\mathrm{ref}}^*) \;\in\; \mathbb{R}^{p \times p}\] positive definite: all \(p\) eigenvalues \(\lambda_1 \geq \cdots \geq \lambda_p\) of \(I_{\theta\theta}^{\mathrm{marg}}\) satisfy \(\lambda_p > 0\).

Compared to v07’s scalar (EB-MARG-ID), the requirement is strictly stronger: in \(\mathbb{R}\) “non-zero Fisher information” is equivalent to positivity, but in \(\mathbb{R}^p\) “non-singularity” (det \(\neq 0\)) does not imply positive definiteness in general — and only the latter ensures that the multivariate Bernstein–von Mises (BvM) limit of the marginal posterior is a proper \(p\)-variate Gaussian with covariance \(n^{-1} (I_{\theta\theta}^{\mathrm{marg}})^{-1}\). The smallest eigenvalue \(\lambda_p\) controls the effective sample size for the most poorly identified direction: when \(\lambda_p\) is small, \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\) has a large variance along the eigenvector of \(\lambda_p\), and the gap between EB and FB widens in that direction (Proposition 7D* (ii)*, §7).

Under \(K > 1\) with slot-indexed \(\theta_{\mathrm{ref}}^{(k)} \in \mathbb{R}^{p_k}\), (EB-MARG-ID)\(_p\) generalizes componentwise: each block matrix \(I_{\theta\theta}^{\mathrm{marg}, k} \in \mathbb{R}^{p_k \times p_k}\) must be positive definite. Under the factored prior \(\pi_\xi = \prod_k \pi_{\xi^{(k)}}\) of §6, this is equivalent to (EB-MARG-ID)\(_p\) for \(L_n^{\mathrm{marg}}\) as a whole; under coupled prior, the cross-block terms must satisfy a strict-diagonal-dominance condition for the full matrix to remain positive definite.

3.2. (PRIOR-FB-WEAK)\(_p\) – Weak FB Prior in \(\mathbb{R}^p\)

(PRIOR-FB-WEAK)\(_p\). The FB prior \(\pi_\Theta\) on \(\theta_{\mathrm{ref}} \in \mathbb{R}^p\) has positive density at \(\theta_{\mathrm{ref}}^*\) and the FB posterior asymptotically dominates the prior componentwise, in the sense that \[\frac{\mathrm{tr}\,\mathrm{Cov}_n(\theta_{\mathrm{ref}}^*)}{\mathrm{tr}\,\mathrm{Cov}_\pi(\theta_{\mathrm{ref}})} \;\to\; 0 \qquad \text{as } n \to \infty.\] Equivalently, the maximum eigenvalue of \(\mathrm{Cov}_n(\theta_{\mathrm{ref}}^*) \cdot \mathrm{Cov}_\pi(\theta_{\mathrm{ref}})^{-1}\) tends to zero.

The trace-based formulation is the natural multivariate analog of v07’s scalar \(\mathrm{Var}_n / \mathrm{Var}_\pi \to 0\). Under \(p = 1\), it reduces exactly to v07’s (PRIOR-FB-WEAK). The eigenvalue-based equivalent formulation has the advantage of being invariant under linear reparametrization of \(\theta_{\mathrm{ref}}\): a reparametrization \(\theta_{\mathrm{ref}} \mapsto A \theta_{\mathrm{ref}}\) leaves the eigenvalues of the covariance ratio unchanged. The trace formulation is sufficient for the asymptotic theory of §§4–6 and is what the framework’s library reports as a diagnostic in gdpar_eb_fit$diagnostics$prior_dominance_ratio (Sub-phase 8.6.B, Charter §3.2).

The hypothesis is satisfied by typical multivariate priors (independent Gaussian per component with sufficient variance; multivariate Cauchy; multivariate Student-\(t\)). It fails when the prior is strongly informative on one component while the data informs that component weakly (Proposition 7D* (iii)*, §7).

3.3. (HIER-COMPLEX)\(_p\) – Multivariate Hierarchical Regularity

(HIER-COMPLEX)\(_p\). The log-marginal likelihood \(\log L_n^{\mathrm{marg}}(\theta_{\mathrm{ref}})\) is twice continuously differentiable on an open neighborhood \(B_\delta(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}) \subseteq \mathbb{R}^p\) of radius \(\delta > 0\), with Hessian bounded uniformly in \(n\). The number of upper-level hyperparameters \(p\) is bounded as \(n\) grows.

This combines v07’s scalar (HIER-COMPLEX) (bounded number of upper-level hyperparameters) with the new \(C^2\) regularity requirement that supports the multivariate Laplace approximation of step (i) in §11.1 of v07 — namely, that the marginal likelihood can be expanded to second order around \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\) in \(\mathbb{R}^p\) with the Hessian inverse providing the Gaussian curvature of the approximation. Under \(p = 1\), \(C^2\) regularity on an interval is the same condition required (implicitly) by v07 §5; the multivariate extension makes it explicit because it is the operative regularity used in the proof of Theorem 7A* (§4) and in the matrix-correction form of Proposition 7B* (§5).

The boundedness of \(p\) as \(n\) grows is the standard hierarchical complexity bound. If \(p\) itself grew with \(n\) (e.g., \(p = p_n \to \infty\)), the analysis would require non-parametric extensions and is deferred to Block 9.x (Open Question O3*-EBFB, §9). The framework’s release-gate (Charter §2.4) treats \(p\) as fixed (chosen by the analyst from the family specification) and \(n\) as the asymptotic index.

3.4. Discussion: Why the Asterisked Hypotheses are Strictly Stronger

The strengthenings are not stylistic conveniences: each is necessary to support the multivariate proof of §4. In particular, (EB-MARG-ID)\(_p\) as positive definiteness (not merely non-zero determinant) is the gateway condition for the multivariate BvM of the marginal posterior; without it, the posterior may concentrate as a degenerate Gaussian on a lower-dimensional subspace, breaking the equivalence with the conditional EB posterior.

4.1. Proof Sketch (Extension of v07’s Reference Argument)

The proof of v07 §5 (the “reference and argument”) generalizes to \(\mathbb{R}^p\) in four steps, each retracing the v07 argument with the strengthened hypotheses of §3.

Step 1 – FB posterior decomposition. As in v07, \[\Pi_n^{\mathrm{FB}}(\xi \mid \cdot) \;=\; \int_{\mathbb{R}^p} \Pi_n^{\mathrm{cond}}(\xi \mid \theta_{\mathrm{ref}}, \cdot)\; \Pi_n^{\mathrm{marg}}(\theta_{\mathrm{ref}} \mid \cdot)\; d\theta_{\mathrm{ref}}.\] The decomposition is algebraically valid for any \(p \geq 1\); the dimensionality enters only through the domain of integration.

Step 2 – Marginal posterior concentration via multivariate BvM. Under (EB-MARG-ID)\(_p\) (positive definite Fisher) and (PRIOR-FB-WEAK)\(_p\) (data swamps prior), the multivariate Bernstein–von Mises theorem (van der Vaart 1998, Chapter 10; Castillo and Rousseau 2015 for smooth functionals; Rousseau and Szabo 2017 §3 for total variation quantification) gives \[d_{\mathrm{TV}}\!\left(\, \Pi_n^{\mathrm{marg}}(\theta_{\mathrm{ref}} \mid \cdot),\; \mathcal{N}\!\left(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}},\; n^{-1} \bigl(I_{\theta\theta}^{\mathrm{marg}}\bigr)^{-1}\right) \right) \;\xrightarrow{P_{\eta_*}}\; 0.\] The covariance \(n^{-1} (I_{\theta\theta}^{\mathrm{marg}})^{-1}\) is the inverse \(p \times p\) matrix Fisher info; its existence requires the positive definiteness of (EB-MARG-ID)\(_p\), not merely the non-singularity. In particular, the concentration rate is uniform across all directions in \(\mathbb{R}^p\), with the slowest direction governed by the smallest eigenvalue \(\lambda_p > 0\) of \(I_{\theta\theta}^{\mathrm{marg}}\). This step is the multivariate generalization of v07’s “the marginal posterior of \(\theta_{\mathrm{ref}}\) concentrates at \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\)”.

Step 3 – Identification of MML maximizer and posterior mode. Under (PRIOR-FB-WEAK)\(_p\), the gradient \(\nabla_{\theta_{\mathrm{ref}}} \log \pi_\Theta(\theta_{\mathrm{ref}}) / n \to 0\) componentwise (the prior’s logarithmic-derivative scaled by \(n^{-1}\) vanishes), so the maximizer of \(L_n^{\mathrm{marg}} \cdot \pi_\Theta\) (posterior mode) and the maximizer of \(L_n^{\mathrm{marg}}\) alone (\(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\)) differ by \(O_P(n^{-1})\). The two coincide to leading order, and the Gaussian centering in Step 2 can be taken at either.

Step 4 – Lipschitz/smoothness of \(\theta_{\mathrm{ref}} \mapsto \Pi_n^{\mathrm{cond}}(\xi \mid \theta_{\mathrm{ref}}, \cdot)\). Under (HIER-COMPLEX)\(_p\) (\(C^2\) in a ball), the conditional posterior depends on \(\theta_{\mathrm{ref}}\) smoothly: for any \(h \in \mathbb{R}^p\) small, \[\left\| \Pi_n^{\mathrm{cond}}(\xi \mid \widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}} + h, \cdot) - \Pi_n^{\mathrm{cond}}(\xi \mid \widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}, \cdot) \right\|_{\bullet} \;\leq\; L\, \|h\|,\] where \(\|\cdot\|_\bullet\) is \(d_{\mathrm{TV}}\) in Regime A and the weak-over-\(\mathcal{G}\) metric in Regime B, \(\|h\|\) is the Euclidean norm in \(\mathbb{R}^p\), and \(L > 0\) is the Lipschitz constant given by the maximum eigenvalue of the conditional Hessian (existence of \(L\) follows from (HIER-COMPLEX)\(_p\)). Combined with the \(O_P(n^{-1/2})\) concentration of Step 2, the difference \[\left| \int \Pi_n^{\mathrm{cond}}(\xi \mid \theta_{\mathrm{ref}}, \cdot)\, \Pi_n^{\mathrm{marg}}(\theta_{\mathrm{ref}} \mid \cdot)\, d\theta_{\mathrm{ref}} \;-\; \Pi_n^{\mathrm{cond}}(\xi \mid \widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}, \cdot) \right|\] in the relevant metric is \(O_P(L \cdot n^{-1/2}) \to 0\). This identifies the FB posterior with \(\Pi_n^{\mathrm{EB}}(\xi \mid \cdot) = \Pi_n^{\mathrm{cond}}(\xi \mid \widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}, \cdot)\) asymptotically. \(\square\)

The four steps reproduce v07’s “reference and argument” with \(\mathbb{R}\) replaced by \(\mathbb{R}^p\) throughout. Each step is the routine multivariate generalization established in the cited literature (Petrone–Rousseau–Scricciolo 2014; Rousseau–Szabo 2017; Castillo–Rousseau 2015; van der Vaart 1998).

4.2. Why Total Variation is Restricted to Regime A (Multivariate Refinement)

The v07 §5 caveat that TV is restricted to Regime A is independent of \(p\): it depends on the dimensionality of \(\xi\), not of \(\theta_{\mathrm{ref}}\). In non-parametric Bayesian inference, the posterior over an infinite-dim \(\xi\) does not converge in TV but in Hellinger or \(L^2(\mu)\) or in the weak-over-\(\mathcal{G}\) metric (van der Vaart and van Zanten 2008; Ghosal–van der Vaart 2017). The multivariate dimension of \(\theta_{\mathrm{ref}}\) does not change this conclusion: the only modification is the rate of concentration of the marginal posterior (which is \(n^{-1/2}\) in each of the \(p\) coordinates).

A subtlety specific to the multivariate setting: even in Regime A (finite-dim parametric \(\xi\)), TV convergence under \(p \geq 2\) requires the additional regularity that the Laplace approximation in \(\mathbb{R}^p\) is uniformly accurate over the relevant neighborhood. This is guaranteed by (HIER-COMPLEX)\(_p\) but is worth noting because the multivariate Laplace approximation can fail numerically (not theoretically) when \(I_{\theta\theta}^{\mathrm{marg}}\) has condition number \(\kappa = \lambda_1 / \lambda_p\) large — see §9 O5*-EBFB.

4.3. Practical Implication of Theorem 7A*

For large \(n\), weakly informative priors on \(\theta_{\mathrm{ref}}\), and well-conditioned \(I_{\theta\theta}^{\mathrm{marg}}\), EB and FB give essentially the same posterior over \(\xi\) in the metric appropriate to the dimensionality of \(\xi\). The dimensionality \(p\) of \(\theta_{\mathrm{ref}}\) and the number of slots \(K\) enter only through (i) the strength of the regularity hypotheses (which scale with \(p\) and \(K\)), (ii) the rate at which the smallest eigenvalue \(\lambda_p\) of \(I_{\theta\theta}^{\mathrm{marg}}\) governs identifiability, and (iii) the numerical difficulty of the multivariate Laplace step (§9 O5*-EBFB).

The choice between EB and FB therefore remains, to first order, primarily computational and methodological: EB is faster (Sub-phase 8.6.B’s cmdstanr::laplace() call typically completes in \(\sim 10\%\) of the FB sampling time on the same model); FB integrates over \(\theta_{\mathrm{ref}}\) honestly (does not require choosing \(\pi_\Theta\)). Higher-order coverage and finite-sample considerations (§§5–6) refine this picture.

4.4. Reduction to v07’s Theorem 7A under \((p, K) = (1, 1)\)

Under \(p = 1\) and \(K = 1\):

• (EB-MARG-ID)\(_1\) = (EB-MARG-ID) of v07 §4 (positivity of scalar Fisher information).
• (PRIOR-FB-WEAK)\(_1\) = (PRIOR-FB-WEAK) of v07 §4 (scalar variance ratio).
• (HIER-COMPLEX)\(_1\) = (HIER-COMPLEX) of v07 §4 augmented with explicit \(C^2\) on an interval (the implicit regularity of v07’s “reference and argument”).
• The four-cell table reduces to its two left-column cells: Regime A (TV) and Regime B (weak-over-\(\mathcal{G}\)).
• The proof sketch (§4.1) reduces to v07’s Steps 1–4 with \(\mathbb{R}\) in place of \(\mathbb{R}^p\) and scalars in place of matrices.

Theorem 7A* therefore strictly contains Theorem 7A of v07: every assertion of v07 §5 is an instance of Theorem 7A* under \((p, K) = (1, 1)\).

4.5. Audit Trail: v07 §5 Reviewed Line by Line

The redaction of this vignette was preceded by a line-by-line audit of v07 §5 intended to detect implicit use of any property valid in \(\mathbb{R}\) but false in \(\mathbb{R}^p\). The four audit-relevant assertions, in the order they appear in v07 §5, are:

1. “The argument decomposes the FB posterior as the average of the conditional posterior given \(\theta_{\mathrm{ref}}\) over the marginal posterior of \(\theta_{\mathrm{ref}}\).” Algebraic identity (mixture decomposition) valid in any dimension. No scalar-specific property used.
2. “Under (PRIOR-FB-WEAK), the marginal posterior of \(\theta_{\mathrm{ref}}\) concentrates at \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\).” Multivariate BvM in \(\mathbb{R}^p\) requires positive definite \(I_{\theta\theta}^{\mathrm{marg}}\); the scalar v07 version requires positivity. The strengthening is explicit in (EB-MARG-ID)\(_p\) of §3. No implicit scalar-specific assumption.
3. “The conditional posterior at the concentration point coincides with the EB posterior.” Definition of \(\Pi_n^{\mathrm{EB}}\); valid in any dimension.
4. “The difference between FB and EB is thus the difference between integrating with respect to a concentrated distribution and evaluating at its mode, which goes to zero in the appropriate metric for the dimension of \(\xi\).” Requires Lipschitz/smoothness of \(\theta_{\mathrm{ref}} \mapsto \Pi_n^{\mathrm{cond}}(\xi \mid \theta_{\mathrm{ref}}, \cdot)\) in the appropriate metric. In \(\mathbb{R}\) this is \(|h|\)-Lipschitz; in \(\mathbb{R}^p\) it is \(\|h\|\)-Lipschitz with Euclidean norm. Norm equivalence in finite dimension ensures transfer. No scalar-specific argument.

Audit conclusion. v07 §5 transfers cleanly to \(\mathbb{R}^p\) under the strictly stronger hypotheses of §3. No implicit error detected. The “reference and argument” delegation of v07 §5 to Petrone–Rousseau–Scricciolo (2014) and Rousseau–Szabo (2017) is faithful: both treat hyperparameters on Polish spaces, of which finite-dimensional \(\mathbb{R}^p\) is the routine specialization.

A minor pedagogical opacity in v07 §6 (not §5) is addressed in §5 of this vignette: v07’s scalar formula \(C_{g,\alpha} \approx (g'(\xi^*))^2 / I_{\theta\theta}^{\mathrm{marg}} \cdot \kappa(\alpha)\) implicitly absorbs a sensitivity factor \(\partial \xi^* / \partial \theta_{\mathrm{ref}}\) into the leading coefficient. The multivariate Proposition 7B* makes this Jacobian (called \(J^\xi\) in §2.3) explicit as part of a sandwich quadratic form, recovering the scalar v07 §6 formula as the \(p = \dim(\xi) = 1\) specialization.

5.1. Derivation Sketch

The under-coverage arises from the fact that EB conditions on a point estimate \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\) and ignores the variability in this estimator, whereas FB integrates over the posterior of \(\theta_{\mathrm{ref}}\). Under the multivariate BvM of Step 2 of §4.1, \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\) has approximate covariance \(n^{-1} (I_{\theta\theta}^{\mathrm{marg}})^{-1}\). The EB-side variance of the plug-in estimator \(g(\xi^*(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}))\), by the multivariate delta method, includes the propagated contribution \[\mathrm{Var}_{\mathrm{EB-FB}}\!\bigl[ g(\xi^*) \bigr] \;\approx\; \bigl[\nabla_\xi g(\xi^*)^\top \cdot J^\xi \cdot \mathrm{Cov}(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}) \cdot (J^\xi)^\top \cdot \nabla_\xi g(\xi^*)\bigr] \;\approx\; \frac{1}{n} \nabla_\xi g(\xi^*)^\top J^\xi (I_{\theta\theta}^{\mathrm{marg}})^{-1} (J^\xi)^\top \nabla_\xi g(\xi^*).\] The EB credible interval, computed from \(\Pi_n^{\mathrm{EB}}(\xi \mid \widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}})\) alone, omits this contribution; FB’s posterior includes it via the marginalization over \(\theta_{\mathrm{ref}}\). The coverage deficit at level \(1 - \alpha\) is therefore proportional to the trace of the omitted covariance contribution scaled by \(\kappa(\alpha)\), which is exactly \(\mathrm{tr}(C_{g,\alpha}^*) / n\). The Edgeworth-type expansion conditions of Bickel–Ghosh (1990) generalized to \(\mathbb{R}^p\) ensure the \(o(n^{-1})\) remainder.

5.2. Reduction to v07’s Scalar Proposition 7B under \(p = \dim(\xi) = 1\)

When \(p = 1\) and \(\dim(\xi) = 1\):

• \(J^\xi \in \mathbb{R}^{1 \times 1}\) is the scalar sensitivity \(\partial \xi^* / \partial \theta_{\mathrm{ref}}\).
• \(\nabla_\xi g(\xi^*) \in \mathbb{R}^{1 \times 1}\) is the scalar derivative \(g'(\xi^*)\).
• \((I_{\theta\theta}^{\mathrm{marg}})^{-1} \in \mathbb{R}^{1 \times 1}\) is the scalar inverse Fisher information \(1 / I_{\theta\theta}^{\mathrm{marg}}\).
• \(C_{g,\alpha}^*\) degenerates to a scalar: \(C_{g,\alpha}^* = \kappa(\alpha) \cdot (g'(\xi^*))^2 \cdot (J^\xi)^2 / I_{\theta\theta}^{\mathrm{marg}}\).

The scalar v07 §6 formula \(C_{g,\alpha} \approx (g'(\xi^*))^2 / I_{\theta\theta}^{\mathrm{marg}} \cdot \kappa(\alpha)\) absorbs \((J^\xi)^2\) into the proportionality constant. Under the convention that \(J^\xi = 1\) for a one-dimensional smooth identity \(\xi \mapsto \xi\) (Carlin–Gelfand 1990 §3), the formulas coincide exactly. In the multivariate extension this absorption is no longer harmless: \(J^\xi\) depends explicitly on the hierarchy’s coupling between \(\theta_{\mathrm{ref}}\) and \(\xi\), and a value \(\|J^\xi\| \gg 1\) amplifies the under-coverage proportionally.

5.3. Effective Sensitivity Diagonalization (Practical Form)

For a practical formula independent of explicit Jacobian computation, diagonalize the sandwich: \[\mathrm{tr}\!\bigl(C_{g,\alpha}^*\bigr) \;=\; \kappa(\alpha) \cdot \nabla_\xi g(\xi^*)^\top \cdot J^\xi (I_{\theta\theta}^{\mathrm{marg}})^{-1} (J^\xi)^\top \cdot \nabla_\xi g(\xi^*).\] Let \(\Sigma_{\mathrm{EB-FB}} := J^\xi (I_{\theta\theta}^{\mathrm{marg}})^{-1} (J^\xi)^\top \in \mathbb{R}^{\dim(\xi) \times \dim(\xi)}\) be the EB-vs-FB covariance gap matrix, computable post-hoc from the EB fit (via finite differences on \(\theta_{\mathrm{ref}} \mapsto \xi^*(\theta_{\mathrm{ref}})\) for \(J^\xi\) and from the Laplace Hessian for \(I_{\theta\theta}^{\mathrm{marg}}\)). Then \[\mathrm{tr}\!\bigl(C_{g,\alpha}^*\bigr) \;=\; \kappa(\alpha) \cdot \nabla_\xi g(\xi^*)^\top \cdot \Sigma_{\mathrm{EB-FB}} \cdot \nabla_\xi g(\xi^*),\] and the EB credible interval for \(g(\xi^*)\) is corrected by inflating its width by a factor approximately \(\sqrt{1 + \mathrm{tr}(C_{g,\alpha}^*)/(n - q)}\) to recover nominal coverage, where \(q\) is the effective degree-of-freedom adjustment (typically the dimension of \(\xi\) relevant to the conditional posterior of \(g(\xi)\)).

This is the form implemented by eb_correction = TRUE (the default in gdpar_eb(), per Charter §2.6 of Sub-phase 8.6).

5.4. Operational Note and the eb_correction Argument

Per Charter §2.6 (Sub-phase 8.6) and per the operative description of v07 §11.1, the framework’s library applies the Proposition 7B* correction by default (eb_correction = TRUE):

• Under \(p = 1\), the scalar correction of v07 §6 is applied.
• Under \(p \geq 2\), the matrix correction \(C_{g,\alpha}^* \in \mathbb{R}^{p \times p}\) of this section is applied via its diagonalized form (§5.3).
• The user can set eb_correction = FALSE to retrieve nominal (un-corrected) intervals, in which case the library emits a gdpar_diagnostic_warning advising of the expected \(O(n^{-1})\) under-cover (per Sub-phase 8.6.B’s implementation).

Caveat: the correction is approximate. As in v07 §6, the correction depends on standard Edgeworth-type expansion conditions and on smoothness of \(g\) at \(\xi^*\). Under non-smooth functionals, near-singular \(I_{\theta\theta}^{\mathrm{marg}}\), or distributions for which Edgeworth expansions fail (e.g., heavy-tailed Student-\(t\) slots; mixture families with near-modes-of-non-identifiability), the correction’s accuracy degrades. The framework recommends FB whenever the standard expansion conditions are in doubt; the EB correction is offered as a partial correction for cases where the user has chosen EB for computational reasons.

6.1. Reduction to v07’s Theorem 7C under \(K = 1\)

Under \(K = 1\): Case A collapses to v07’s Theorem 7C with \(K\) replaced by \(J\) (in v07’s notation, the symbol \(K\) stood for the number of exchangeable units within the single slot; in this vignette, \(K\) stands for the number of distributional slots and \(J\) stands for the number of groups within each slot). The relabeling does not change content: the gap \(B_{J, 1}\) tends to zero as the number of within-slot exchangeable units grows, reproducing v07’s \(B_K \to 0\) as \(K \to \infty\) verbatim. Case B is vacuous under \(K = 1\) (no cross-slot coupling exists).

6.2. Two Regimes of Practical Use (Multi-Slot Refinement)

• \(J\) large within slots, \(K\) moderate: per-slot EB risk approaches per-slot FB risk; EB is the practical choice slot-by-slot. The aggregate bound is dominated by the per-slot gaps.
• \(J\) small, \(K\) small: the gap \(B_{J, k}\) is non-negligible per slot, and under coupled prior the residual \(\mathcal{B}_{\mathrm{coupling}}\) adds an explicit cross-slot penalty. FB is preferred to capture the full covariance.
• \(K\) large, \(J\) large: per-slot bound tightens uniformly; if the prior is factored (Case A), EB is asymptotically efficient. If the prior is coupled (Case B), the cross-slot interaction term decays at \(J^{-1}\) but adds a constant overhead at any fixed \(J\), advising caution.

6.3. Caveat (Squared-Error Loss; Coverage Picture Unchanged)

As in v07 §7, Theorem 7C* is for squared-error loss on point estimates. For interval coverage, Proposition 7B* (§5) shows that EB systematically under-covers regardless of \(K\) or \(J\); the choice between EB and FB therefore depends on whether the user prioritizes point estimate accuracy under squared-error loss (EB acceptable for large \(J\), possibly large \(K\) under Case A) or coverage calibration (FB preferred uniformly).

6.4. Reference

This is the standard compound decision result of Robbins (1956) extended to \(K \geq 1\) slots. The per-slot bound of Case A is the routine slotwise application of Brown (2008) and Efron (2019, Theorem 1.1). The cross-slot residual of Case B is derived from the second-order Taylor expansion of the coupled prior’s log-density at the truth and is, to the author’s knowledge, not previously stated in this form in the literature; it is canonized here as the AMM specialization of the multi-hyperparameter compound decision setting.

7.1. Reduction to v07’s Proposition 7D under \((p, K) = (1, 1)\)

Under \(p = 1\) and \(K = 1\), the four asterisked conditions reduce to v07’s Proposition 7D conditions (i)–(iv) with the following correspondences:

The renumbering (v07’s (iii) becomes v07b’s (i)*) reflects the ordering of practical importance in \(\mathbb{R}^p\): multimodality is the most frequent failure mode of multivariate EB and is treated first.

9.1. (O1*-EBFB) Adaptive Choice between EB and FB Based on Data (Multivariate)

v07’s O1-EBFB is the data-driven adaptive selection between EB and FB based on a diagnostic. The multivariate extension requires a diagnostic that handles the multi-dimensional Fisher information, the per-slot block structure under \(K > 1\), and the inter-init dispersion of multi-start Laplace. A candidate diagnostic is the ratio of EB-side coverage gap \(\mathrm{tr}(C_{g,\alpha}^*)/n\) (computable post-hoc from gdpar_eb_fit) to the FB-side variance of \(g(\xi^*)\) obtained from a short pilot FB run; the framework will release a candidate diagnostic in Block 9.x as gdpar_eb_diagnostic() after sufficient empirical calibration. Partial work: Petrone–Rousseau–Scricciolo (2014, §7) discuss adaptive empirical Bayes in the abstract; the AMM specialization in \(\mathbb{R}^p\) with heterogeneous slots is open.

9.2. (O2*-EBFB) Higher-Order Coverage Correction for Multivariate EB

v07’s O2-EBFB is the search for higher-order corrections (Edgeworth-type, bootstrap-based) beyond the first-order Proposition 7B correction. The multivariate extension requires corrections that handle (i) the sandwich matrix form of \(C_{g,\alpha}^*\), (ii) the per-slot tensor under \(K > 1\), and (iii) the propagated uncertainty from the multi-start Laplace inter-init dispersion. Multivariate Edgeworth expansions are technically more delicate (Bhattacharya and Ghosh 1978; Hall 1992 Chapter 5) and bootstrap calibration in \(\mathbb{R}^p\) requires more replications to control simulation noise. The framework’s current implementation uses the first-order matrix correction (§5); higher-order corrections are deferred to Block 9.x.

9.3. (O3*-EBFB) Behavior under Model Misspecification (Multivariate)

v07’s O3-EBFB is the question of which of EB or FB degrades more gracefully under misspecification. Under multivariate hyperparameter and heterogeneous slots, the question further refines: does misspecification of the per-slot prior \(\pi_{\xi^{(k)}}\) in one slot affect the EB estimator of \(\theta_{\mathrm{ref}}^{(k)}\) in only that slot (a slot-local insulation property of Case A of Theorem 7C*) or across slots through the coupled prior of Case B? Partial empirical observation from the AMM regression tests (tests/testthat/test-eb_misspecification.R to be written in Sub-phase 8.6.D) will inform; full theoretical treatment is deferred.

9.4. (O4*-EBFB) Computational Competitiveness of FB at Moderate \(n\) under \(K \cdot p\) Large

v07’s O4-EBFB is the question of when modern HMC narrows the computational gap enough to make FB the default even for very large \(n\). The multivariate extension asks: how does the FB sampling time scale with the joint dimensionality \(K \cdot p\)? Empirical observation from Bloque 7 benches (5 species × 4 NE-USA sub-regions, single-slot \(K = 1\), \(p = 1\)): the framework’s FB sampling took \(\sim 6600\) s vs. mgcv’s \(\sim 0.7\) s; under \(K \cdot p \geq 4\) the gap may grow super-linearly. The threshold at which FB becomes infeasible per \(K \cdot p\) is problem-dependent and remains open.

9.5. (O5*-EBFB) Numerical Anti-Fragility of Multivariate Laplace under Non-Conjugacy (New)

The most important new open question raised by the multivariate extension is the formal characterization of when cmdstanr::laplace() fails numerically in \(\mathbb{R}^p\) under non-conjugate family + prior combinations. The failure modes are:

• Singular or ill-conditioned Hessian of the marginal log-likelihood at the candidate \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\). The condition number \(\kappa(H) = \lambda_1(H) / \lambda_p(H)\) controls the relative error of the inverse Hessian used in the Laplace approximation; when \(\kappa(H) \gtrsim 10^{10}\), the Gaussian approximation around \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\) is unreliable and the EB credible intervals are invalid.
• Non-positive-definite Hessian when the L-BFGS optimization converges to a saddle-point or to a local maximum of \(L_n^{\mathrm{marg}}\) that is not a true global maximum.
• Multiple disconnected modes when \(L_n^{\mathrm{marg}}\) has more than one local maximum in \(\mathbb{R}^p\) (recall §7 Proposition 7D* (i)*: multimodality is generic for \(p \geq 3\) under non-conjugacy).
• Slow gradient decay along ridge directions when one or more eigenvalues of the marginal Hessian are near zero, causing L-BFGS to oscillate rather than converge.

The framework’s anti-fragility strategy (canonized in Charter §2.8 of Sub-phase 8.6) addresses these four failure modes with four operational components:

1. Preventive detection of ill-conditioned Hessian. Before accepting the output of cmdstanr::laplace(), the library computes the condition number \(\kappa(H)\) of the marginal Hessian at the candidate \(\widehat{\theta}_{\mathrm{ref}}^{\,\mathrm{EB}}\). If \(\kappa(H) > \kappa_{\mathrm{threshold}}\) (default \(10^{10}\), configurable via laplace_control$kappa_threshold), the EB fit aborts with class gdpar_eb_numerical_error and an explicit diagnostic message identifying the most poorly conditioned direction (via the eigenvector of \(\lambda_p\)).
2. Adaptive Levenberg-Marquardt ridge. If \(|\det(H)| < \epsilon_{\mathrm{LM}}\) (Hessian effectively singular) or \(H\) is not positive-definite (a true Newton step would diverge), the library substitutes \(H \leftarrow H + \lambda I\) with an adaptive ridge \(\lambda\) chosen via the standard L-M heuristic (start with \(\lambda_0 = \mathrm{tr}(H) \cdot 10^{-3}\); geometrically increase by 10 until \(H + \lambda I\) is positive-definite with condition number below threshold; report the used \(\lambda\) in gdpar_eb_fit$diagnostics$lm_perturbation).
3. Multiple random restarts. The L-BFGS optimization is restarted from \(M\) independent random inits within a bounding box derived from the prior \(\pi_\Theta\) (default \(M = 5\); configurable via laplace_control$multi_start_M$). The library retains the init that maximizes the marginal log-likelihood and reports the inter-init dispersion (coefficient of variation of marginal log-likelihood at each converged init) in gdpar_eb_fit$diagnostics$multi_start_dispersion. A coefficient of variation \(> 0.05\) across inits signals likely multimodality and the library issues a gdpar_diagnostic_warning advising the user to consider FB.
4. Documented fallback. If all \(M\) inits fail (all produce singular Hessian or all diverge), gdpar_eb() aborts with class gdpar_unsupported_feature_error and a diagnostic message recommending the user to use gdpar() (FB) instead. The library never returns a silently invalid estimator.

The formal characterization of the conditions under which cmdstanr::laplace() fails (independent of the empirical heuristics above) is open. A candidate framework is to express failure as a function of two intrinsic quantities of the AMM specification: (a) the condition number \(\kappa(I_{\theta\theta}^{\mathrm{marg}})\) of the population Fisher matrix, and (b) the Morse index spread of \(L_n^{\mathrm{marg}}\) in \(\mathbb{R}^p\) (the variance of the number of negative eigenvalues of the Hessian across the interior of the admissible domain \(\Theta\)). A precise theorem of the form “if \(\kappa \leq \kappa_*(p)\) and Morse spread \(\leq \mu_*(p)\), then Laplace succeeds with probability \(\geq 1 - \delta_n(p)\) where \(\delta_n(p) = O(n^{-1/p})\)” is the open form; the constants \(\kappa_*(p)\), \(\mu_*(p)\), and \(\delta_n(p)\) are unknown for general non-conjugate AMM specifications and require systematic empirical-and-theoretical investigation. This is the primary open theoretical question that motivates the anti-fragility strategy as a working substitute pending formal resolution.

11.1. The Stein Operator and the Stein Kernel

For a target distribution \(p\) with score function \(s_p(x) = \nabla_x \log p(x)\), the Stein operator \(\mathcal{A}_p\) acting on a vector-valued function \(f\) is \[\mathcal{A}_p f(x) \;=\; \langle s_p(x), f(x) \rangle + \nabla \cdot f(x).\] Stein’s identity asserts \(\mathbb{E}_{x \sim p}[\mathcal{A}_p f(x)] = 0\) for all \(f\) in a sufficiently smooth class. Given a base reproducing kernel \(k(x, y)\), the Stein kernel is \[k_p(x, y) \;=\; \langle s_p(x), s_p(y) \rangle k(x, y) + \langle s_p(x), \nabla_y k(x, y) \rangle + \langle s_p(y), \nabla_x k(x, y) \rangle + \mathrm{tr}\!\left(\nabla_x \nabla_y^{\top} k(x, y)\right),\] and the KSD of a sample \(\{x_i\}_{i=1}^n \sim Q\) from a target \(p\) is the V-statistic \[\mathrm{KSD}(Q, p) \;=\; \sqrt{\max\!\left\{ 0,\; \tfrac{1}{n^2} \sum_{i, j} k_p(x_i, x_j) \right\}}.\] A value close to zero indicates \(Q\) is well-aligned with \(p\); positive values indicate joint distributional deviation.

11.2. Empirical Gaussian Target (B9.4 iteration) and Future Extension

The KSD requires the score of the target. In this iteration the target is the empirical Gaussian Laplace approximation of the FB posterior over the common \(\xi\)-variables: \(\hat\mu = \mathbb{E}^{\mathrm{FB}}[\xi]\), \(\hat\Sigma = \mathrm{Cov}^{\mathrm{FB}}(\xi)\), so the Gaussian score is the linear closed-form \[s(x) \;=\; -\hat\Sigma^{-1}(x - \hat\mu).\] The implementation evaluates \(\mathrm{KSD}(\Pi_n^{\mathrm{EB}}, \mathcal{N}(\hat\mu, \hat\Sigma))\) in closed form (no auxiliary MCMC required for the target score). This is sufficient to detect EB-vs-FB deviations of joint first- and second-order moment structure. The full-KSD variant against the actual FB target via the Stan model’s grad_log_prob() method (a one-sample KSD of EB samples against the true non-Gaussian FB posterior) is a documented extension reserved for Block 9.x.

11.3. Base Kernel and Bandwidth

The default base kernel is the inverse multi-quadric (IMQ) of Gorham-Mackey \[k(x, y) \;=\; \left(h + \|x - y\|^2\right)^{\beta}, \qquad \beta \in (-1, 0), \quad \beta_{\mathrm{default}} = -\tfrac{1}{2},\] which is the canonical choice for distinguishing distributions whose densities differ in tails: Gorham-Mackey Theorem 8 establishes that the IMQ kernel detects mismatch with a dimension-independent rate for log-concave targets, a property the RBF kernel \(k(x, y) = \exp(-\|x - y\|^2 / (2 h))\) does not enjoy. The bandwidth \(h\) (in squared units of \(x\)) defaults to the median heuristic: \(h = \mathrm{median}\!\left( \{ \|x_i - x_j\|^2 \}_{i < j} \right)\) over the FB sample, which is dimension-adaptive. A fixed bandwidth is also accepted as a configuration option.

11.4. ESS-Weighted Variant

When MCMC chains are short relative to the integrated autocorrelation time of \(\xi\), the V-statistic \(\tfrac{1}{n^2}\sum_{i,j} k_p(x_i, x_j)\) understates the asymptotic KSD because successive draws are not independent. The ess_weighted = TRUE variant addresses this by thinning both EB and FB samples to \(\min(\widehat{\mathrm{ESS}}_{\mathrm{EB}}, \widehat{\mathrm{ESS}}_{\mathrm{FB}})\) rows (per-variable basic ESS via posterior::ess_basic), yielding an unbiased V-statistic on the thinned samples. The thinning seed is user-configurable for reproducibility.

11.5. Operational Reading: gdpar_compare_eb_fb and gdpar_ksd_joint are Complementary

The two diagnostics measure orthogonal aspects of the EB-vs-FB discrepancy:

• gdpar_compare_eb_fb() returns a marginal-by-parameter TV table: \(\widehat{\mathrm{TV}}(\Pi_n^{\mathrm{EB}, k}, \Pi_n^{\mathrm{FB}, k})\) for each scalar parameter \(k\) separately. This detects deviations in univariate marginals but is blind to joint dependence (e.g., two EB and FB posteriors with identical marginals but different correlation structure give the same marginal TV).

• gdpar_ksd_joint() returns a single joint KSD scalar that integrates over the full joint \(\xi\)-vector. This detects deviations in the joint dependence structure (correlations, copula geometry) that the marginal TV misses, at the cost of a coarser per-parameter resolution.

The recommended operational reading is to compute both: a small marginal TV plus a large joint KSD is the signature of an EB approximation that gets the per-coordinate spread right but distorts the cross-coordinate correlations. The converse (large marginal TV, small joint KSD) is uncommon in practice but indicates per-coordinate location/scale deviations that nevertheless preserve the joint dependence skeleton.