First public release of gdpar. It implements the General Dynamic Parameter
framework, in which each unit's parameter is decomposed as
theta_i = theta_ref + Delta(x_i, theta_ref) around a population reference, with
Delta following the canonical Additive--Multiplicative--Modulated (AMM) form
a(x) + b(x) * theta + W(theta) x. Estimation is via Path 1 (hierarchical
Bayesian inference through 'Stan' / 'cmdstanr'); Paths 2 (varying-coefficient)
and 3 (hypernetwork) are specified at reference grade and queued for future
versions.
Packaging note (pre-CRAN correction): the fitting entry points now run all
input validation before the optional cmdstanr dependency is required, which
is checked only immediately before it is used. The package and its full test
suite therefore work on machines where the Suggests package cmdstanr is not
installed. This is a packaging-only change with no effect on results when
cmdstanr is available.
This release consolidates the development cycles documented in detail below (Sessions 8.3--8.6, the Block 8.4 audit, and the Block 9 re-validation):
K >= 1 slots. Gaussian, Student-t, Gamma,
Beta, Tweedie, Poisson, negative-binomial, zero-inflated and hurdle families,
plus heterogeneous per-slot families, with B-spline W bases and arbitrary
covariate dimension p.K > 1 and p > 1 regime.gdpar does
not model the dependence; only its inference is made robust to it.Validation. Block 9 closed adversarial internal re-validation, a synthetic
recovery benchmark against mgcv / brms / INLA / rstanarm with known ground truth
(3200 cells), and an organic eBird benchmark (80 cells). gdpar leads on the
distributional structure it models (heavy-tail quantiles, zero-inflated count
means, heteroscedastic scale), is statistically indistinguishable from the
strongest competitor (mgcv) on generic predictive accuracy, is robust where it
makes no modelling claim (e.g. autocorrelation), and is the most computationally
expensive method in the roster. Reports under inst/benchmarks/results/.
Detailed per-session development notes follow.
gdpar_fit), not only a scalar Empirical-Bayes
fit (gdpar_eb_fit). gdpar_dependence_diagnostic(),
gdpar_dependence_robust(), gdpar_spatial_dependence_diagnostic() and
gdpar_spatial_dependence_robust() are unchanged in name and signature; they
simply work on either path now, closing the EB/FB asymmetry left open at the
start of Axis 2. The K > 1 / p > 1 paths remain deferred on both EB and FB.theta_ref, a_coef, b_coef, W_raw, the last on its
raw scale, for parity with the Empirical-Bayes extractor; the theta_ref
hyperparameters are excluded). robust_se is the block-bootstrap SD of the
per-refit posterior means and se_ratio = robust_se / model_se remains a
like-for-like (SD-vs-SD) ratio. Each refit re-runs the full HMC and is
markedly more expensive than an Empirical-Bayes refit (the opt-in cost message
says so).refit_diagnostics field on the
robust result: max_rhat, min_ess_bulk, n_divergent_refits,
n_high_rhat_refits). Refits are never excluded or down-weighted (that would
be a non-random, biasing screen); a single informational note is emitted only
when a refit's R-hat clearly exceeds 1.05. This applies to both paths.robust_se;
under an informative prior se_ratio < 1 is benign prior regularization, not
an overstated model SE; and the Empirical-Bayes theta_ref point estimate (the
Laplace mode) and the full-Bayes one (the posterior mean) are different
estimands that coincide only asymptotically. A bagged / widened posterior
(BayesBag) is recorded as a deferred lateral, not adopted. The full-Bayes
design was cross-lineage-reviewed; the review's concrete proposals (a hard
minimum-ESS abort, an estimator-replacing bag) were audited and not adopted
in favour of the more robust informational accounting, while its documentation
findings (the se_ratio < 1 and mode-vs-mean subtleties) were adopted.gdpar_dependence_robust() and gdpar_spatial_dependence_robust() gain an
opt-in data-driven block size, block_length = "auto" / block_size = "auto". The default behaviour is unchanged (the rate-only n^(1/3) /
n^(1/4)); "auto" replaces the rate's constant with a value selected from
the fitted residuals, with no extra Stan refit, and the rate as the fallback.
Both functions also gain residual_type / randomize_seed, used only to feed
the "auto" selector. The chosen value and method are reported
(block_length / block_size and the new block_length_method /
block_size_method).b_opt = (2 ghat^2 / D)^(1/3) n^(1/3) with the
overlapping moving/circular-block constant D = (4/3) spec^2; longer
dependence yields a longer block, white noise yields unit blocks, and a
degenerate series falls back to the rate. A p << n caveat is documented.g, cheap (no-refit)
spatial block resamples give the bootstrap variance of the design-weighted
residual functionals (the coefficient's influence directions); g minimises
an empirical mean-squared error — squared bias anchored at the largest blocks
(least biased) plus a 1/n_tiles sampling-variance term (Lahiri 2003). A
decorrelating cross-lineage review supplied the MSE skeleton; two of its
choices were corrected after an audit and empirical validation (the bias
anchor, and the variance term), because the proposed forms would have made the
selector either anticonservative or non-adaptive. A single isotropic g is
used (anisotropy is a documented, deferred limitation).gdpar_spatial_dependence_diagnostic() (new, exported) — Moran's I for the
residuals of a scalar Empirical-Bayes fit. The spatial sibling of
gdpar_dependence_diagnostic(): it builds a row-standardized spatial weight
matrix (default k-nearest-neighbour, the data-driven k = max(4, min(round( log n), n - 1)); or a distance band; or a user-supplied W carrying domain
knowledge), hand-rolls Moran's I in base R (no spdep / sf dependency),
and reports a significance test — a two-sided permutation test by default
(robust to non-normal Dunn-Smyth residuals and asymmetric weights), or the
cheaper analytic Cliff-Ord normal approximation (which warns under an
asymmetric W). Guards cover isolated locations (zero-weight rows warn and
return NA), duplicate coordinates (permitted), small n (hard/soft
warnings), and the lon/lat (project-first) and model-misspecification caveats
are documented. gdpar does not model the spatial dependence; this only
makes the violation visible.
gdpar_spatial_dependence_robust() (new, exported) — spatial
block-bootstrap-by-refit standard errors. Refits the model on B spatial
block-bootstrap resamples and reports bootstrap SEs and percentile intervals
alongside the model-based ones, in the working-independence + robust-variance
spirit of Liang & Zeger (1986) (point estimates unchanged). Default scheme is
non-overlapping tiled blocks with a randomized grid origin per
replicate (Politis-Romano-Lahiri), plus an opt-in overlapping moving
scheme. The default block side is g = max(2, round(n^(1/4))) cells per axis,
the d = 2 case of the variance-MSE-optimal rate M ~ n^(d/(d+2)) points per
block, which reduces exactly to the temporal n^(1/3) block length at
d = 1 (derivation and the registered cross-lineage dissent toward n^(1/6)
are in ?gdpar_spatial_dependence_robust). Collinear coordinates abort; a
single-cell collapse warns.
Internal refactor (regression-gated). The temporal gdpar_dependence_ robust() and the new spatial function share one extracted engine
(.gdpar_dependence_robust_engine) and one residual extractor
(.gdpar_dependence_residuals); the temporal path is verified bit-identical
pre/post refactor by a fixed-seed regression gate (a frozen pre-refactor copy
vs. the live engine-backed function, same robust-SE table) plus a Stan-free
golden that locks the engine's RNG-consumption order. No change to the
temporal API or output.
gdpar_geom_laplace() (new, exported) — the Laplace approximation as a
first-class capability. Given a gdpar_geom_target, it climbs to the
posterior mode (L-BFGS-B warm start + a modified-Newton polish when an exact
Hessian is exposed, reading the same target and gradient the sampler uses),
forms the precision M = -Hessian (the observed information, exact cmdstan
Hessian or finite differences), reports its covariance M^{-1}, optional iid
draws mode + M^{-1/2} z, and — crucially — a fidelity diagnostic of the
Gaussian against the true posterior over the same draws: the
importance-sampling ESS, the PSIS Pareto-k (when loo is installed), and
the mean/max log-density drop against its Gaussian expectation d/2. These are
distilled into a single scalar label "good" / "poor" / "very_poor" so the
approximation is never mistaken for exact MCMC. The un-floored condition
number is reported and a saddle (non-positive-definite curvature) is flagged
loudly. This promotes the validated extractor core of RG.7
(inst/benchmarks/scripts/rg7_laplace_elpd.R) into R/.
gdpar_geom_orchestrate() and gdpar_geom_fit() gain laplace_fallback = FALSE (opt-in) and laplace_draws = 0L. When opted in and a run ends in a
certified limit (the genuinely non-Gaussian canyon the geometry ladder
cannot sample — RG.7's eBird Tweedie count), the orchestrator attaches the
Laplace approximation ($laplace) and relabels the status
"certified_limit_laplace": the sampling limit still stands and a labelled,
fidelity-diagnosed Laplace posterior is provided — exactly the mgcv/REML and
INLA/Laplace competitor-parity regime. On the out-of-scope path the fallback is
attached only when the curvature at the mode is genuinely positive-definite
(status "out_of_scope_laplace"); otherwise the certificate is left untouched
(no overreach). The default FALSE leaves the output bit-identical (no
$laplace, the original status set), so the goldens stay bit-identical.
gdpar_geom_bridge() (new, exported) — the durable, path-agnostic core.
Turns an already-fitted gdpar object into the inputs the geometry-adaptive
controller gdpar_geom_orchestrate() consumes, without touching gdpar().
It reads the fit's compiled cmdstan model and Stan data, exposes the
standalone log_prob / grad_log_prob / hessian methods, derives the
unconstrained dimension and a posterior-mean warm-start, and returns a
(target, geom_target, fisher, reference) tuple. Because the diagnostic needs
a re-samplable model and a CmdStanMCMC object cannot be re-sampled, the
bridge recompiles one from the fit's own Stan source ($code(), with the
standalone methods; cmdstanr's content-hash cache makes it a cache hit), while
the engine target reuses the fitted object directly. The Hessian compilation is
best-effort: if higher-order autodiff will not compile (models with custom
densities, the Tweedie lpdf being the case RG.7 targets), it falls back to
gradient-only methods, which still serve the Euclidean, dense and
sub-Riemannian (expected-Fisher) levels. This is the tool RG.7 points at the
real Tweedie count of benchmark 9.2.O.
gdpar_geom_fit() (new, exported) — the one-call ergonomic entry. A
sister of gdpar() (not an internal branch) for the K-individual path: it
builds and compiles the model through the shared seam .gdpar_K_build(),
then runs the orchestrator instead of the default NUTS fit, returning a rich
gdpar_geom_fit object. Print methods for both new classes are provided.
.gdpar_K_build() (internal) extracted from .gdpar_K(). A
behaviour-preserving refactor of the model-building phase (everything up to but
not including cs_model$sample): .gdpar_K() and gdpar_geom_fit() now share
a single source (no duplication, no throwaway computation). The default
branch (compile_model_methods = FALSE) issues the byte-identical
cmdstanr::cmdstan_model(stan_path) call used before, so the K-path goldens
stay bit-identical (verified: the K=2 Gaussian golden re-fits with
max|delta| = 0); gdpar_geom_fit() passes compile_model_methods = TRUE to
expose the standalone methods.
Pedagogical vignette vop08_geometric_robustness.Rmd (new). An abundant,
two-level narrative of the whole Block RG (RG.0 → RG.6): the conceptual bridge
to the user's geometry manuscript (with an organic critical reading of its §12
hierarchy and its §14 "future statistical line", and of the ORPHEUS §16.3
demonstrated-vs-conjectured honesty), the diagnostic with size-invariant
signals, the five sampler levels, the orchestrator with its armour and its
certified limit, Option A, and the integration. The evaluated chunks are
cheap (suite metadata, the budget/criteria/thresholds, and the geometry engine
on pure-R closure targets — no Stan); the heavy cmdstan-backed runs are shown
with eval = FALSE and reproduced end to end by the gated script
inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1).
Opt-in throughout: the default gdpar() fit path is untouched, so the goldens
remain bit-identical (canonized as D97, session B9.31).
Compaction of the per-slot additive scale sigma_a_k to the slots that
carry free a coefficients (J_a_free > 0). A distributional / multivariate
slot whose a() component contributes no free coefficient — either because the
slot declares no a() at all (use_a_k == 0, e.g. phi ~ 1, p ~ 1) or
because it declares an intercept-only a() (a = ~ 1, fully absorbed into the
anchor theta_ref) — used to declare a sampled-but-unused half-prior scale
sigma_a_k[k]. That scale is a flat direction (a non-identified nuisance
dimension) that the no-U-turn sampler explores pointlessly and that, combined
with other pathologies, produced the rhat = Inf first diagnosed in session
B9.21. The Stan templates now compute, in transformed data, the count
n_sigma_a of slots carrying free a coefficients and an index map
sigma_a_idx (the compacted position per slot, 0 if none), declare
sigma_a_k with length n_sigma_a, and index it as sigma_a_k[sigma_a_idx[k]]
wherever the additive coefficients are reconstructed or given their hierarchical
prior. Touched: the three K-individual canonical pieces
(amm_canonical_distrib_K.stan, amm_canonical_eb_conditional_K.stan,
amm_canonical_eb_marginal_K.stan) via the shared codegen generate_a_blocks_K,
the three K×p multivariate pieces (amm_canonical_pmulti_KxP.stan,
amm_eb_conditional_KxP.stan, amm_eb_marginal_KxP.stan) inline, and the Path C
K×p random-init helper so the init vector matches the compacted length.
Methodological decision (session B9.30, the cornerstone rule of maximum multidimensional robustness; canonized as D96):
The criterion is J_a_free > 0, not use_a_k == 1. The mathematically
correct definition of when sigma_a_k exists is "there is at least one free
a coefficient for it to scale". This single criterion removes the flat
direction in both sub-cases (no a() and intercept-only a()), so it is
strictly more robust than a use_a_k-based rule, which would miss the
intercept-only case. Normalising the inconsistency upstream in amm_spec
(forcing use_a_k = 0 for an intercept-only a()) was rejected as an
over-correction: "an a() was declared" and "an a() has free coefficients"
are legitimately distinct facts, and the truth that governs sigma_a_k lives
at the Stan level, where J_a_free > 0 states it exactly.
Bit-identity and the one re-bootstrapped golden. When every slot carries
free a coefficients (the case of every golden whose formula puts a(x) on
each slot) n_sigma_a == K and sigma_a_idx is the identity, so the generated
model is mathematically identical and the frozen golden draws stay bit-exact
(verified by re-fit-and-compare on the K, K×p and EB paths). The sole
exception is the K×p student_t golden g_KxP_004, whose nu slot declares
an intercept-only a() (a = ~ 1) and therefore carried a flat
sigma_a_k[3]; Option A removes it. That golden was re-bootstrapped, and the
change was validated as a faithful marginalisation: sigma_a_k[3] appears
only in its own prior (never in the likelihood or any other parameter's
prior), so it is a-priori independent of the rest and its removal leaves the
joint over the surviving parameters mathematically invariant; the re-fit shows
exactly one fewer parameter and no surviving marginal mean disagreeing beyond
Monte-Carlo error. A regression guard test fixes that a slot with an
intercept-only a() yields a sigma_a_k of length < K.
Geometry-adaptive orchestrator (gdpar_geom_orchestrate()). The closed loop
of the Block RG charter: an opt-in controller that diagnoses the geometry of a
posterior (reusing gdpar_geometry_diagnostic()), classifies the pathology,
selects a level of the sampler hierarchy (euclidean_diagonal →
euclidean_dense → riemannian → relativistic → sub_riemannian), samples
with that level (reusing gdpar_geom_hmc() and the RG.2/RG.3/RG.4 metrics),
re-diagnoses, and either escalates the level or emits a certified limit. It
is standalone and does not touch the package's fit path, so the default branch
is bit-identical (goldens intact). Returns an object of class
gdpar_geom_orchestration (status one of "resolved", "certified_limit",
"out_of_scope"), carrying the per-round decision/sampling ledger, the winning
metric, the best result so far, the budget spent and a reproducibility block.
Methodological decisions of this sub-phase (session B9.29, four hinges consulted upfront and decided by the user under the cornerstone rule of maximum multidimensional robustness):
Level selection — the combined map. The diagnostic's transparent,
calibrated rule-based classifier is the primary selector; a continuous
proximity score over the size-invariant signals is a second, independent
estimator. When they agree and confidence is high, the discrete level is used
directly; when they conflict or confidence is low, the controller starts at
the lower of the two candidate levels and lets the escalation climb — robust
at the funnel/heavy-tail border the RG.1.c calibration found mutually
confusable. A user level_map / entry_level overrides this.
Escalation — the closed-loop, re-diagnosis-guided controller, fully armoured. On a failed level the controller re-diagnoses with a fresh pilot and lets the new signature pick the next level (a re-diagnosis-guided upward jump that may skip ladder rungs), under a monotone ratchet (never below a level already tried), a visited-state memo (a provably acyclic walk), a global round cap, per-level caps, a no-progress (stall) detector, hierarchical budget accounting with cost-aware admission and cheap-probe-then-full tiering (the successive-halving spirit; Li et al. 2018), a coarse step-size search per level (so a level is not failed for a mis-set step), a per-fit wall-time watchdog, graceful degradation (the best result so far is always returned), deterministic seeding (the adaptive trajectory is bit-reproducible, so a re-run retraces it — the resumability guarantee, with an optional atomic ledger checkpoint), a multi-signal success gate with hysteresis, a frozen-level final sampling phase (preserving the exactness / ergodicity of the returned draws; cf. the diminishing-adaptation condition of Roberts & Rosenthal 2007), and the capability-subsumption guarantee (level L+1 strictly generalises level L, so the monotone ladder fallback is a proven monotone-improvement backstop when an adaptive jump misfires). Every level is Metropolis-exact, so a wrong jump wastes only a bounded amount of budget, never correctness.
The certified limit (gdpar_geom_certificate). When the budget is
exhausted without success the controller returns a first-class, reproducible,
falsifiable certificate: the demonstrated evidence in three rigour layers
(algebraic = the geometry; statistical = the per-level sampler diagnostics;
numerical = the budget and fits) plus a conjectured prescription — the
smallest change predicted to break the limit, tagged as a conjecture in the
sense of the demonstrated/conjectured honesty convention and made falsifiable
by a re-run test. A diagnosis pointing outside the geometry ladder
(multimodality, a boundary, a flat direction) short-circuits to a certificate
naming the proper remedy (tempering, reparametrisation, Option A), without
overreaching by pretending a geometry metric fixes it.
Scope — standalone with a rich return. The controller never touches
gdpar(); the return carries everything a future integration (RG.6) will need
(the winning metric, the ledger, the reproducibility block) without freezing a
contract that the integration might have to redesign. The statistical-layer
validation exercises a real compiled cmdstan model end to end.
The tunable factories gdpar_geom_orchestrate_budget() (budget and stopping
rule) and gdpar_geom_orchestrate_criteria() (the multi-signal success gate)
expose the controller's data so it can be re-calibrated. Anchors: Betancourt
(E-BFMI, sampling pathologies); Girolami & Calderhead 2011; Roberts & Rosenthal
2007 (diminishing adaptation / containment); Li et al. 2018 (Hyperband /
successive halving); Xu et al. 2008 (SATzilla, runtime-predicted algorithm
selection with a presolver fallback). gdpar integrates these established
patterns; it does not claim to invent them.
Finsler / relativistic metric (gdpar_geom_metric_relativistic()).
The level-4 geometry of the Block RG hierarchy and the remedy for heavy tails
and directional anisotropy (the G3_heavy_tails target, not the count). A
Gaussian kinetic energy gives the unbounded velocity M^{-1} p, so a large
momentum drawn in a heavy tail produces an arbitrarily large position step
that a fixed-step integrator cannot follow (overshoot, divergences, a step
size forced to the stiffest region). The relativistic kinetic energy caps the
velocity at a finite speed c, taming the tails and the ill-conditioning
while staying exact. Coupled to the position-dependent Riemannian metric
M(theta) of RG.3 (the chosen, maximally robust form), the kinetic energy is
the relativistic energy of a particle of rest mass m on the statistical
manifold,
K = c sqrt(p' M(theta)^{-1} p + m^2 c^2) + 0.5 log det M(theta).
Methodological decisions of this sub-phase (Lu, Perrone, Hasenclever, Teh and Vollmer 2017 for the relativistic kinetic; Livingstone, Faulkner and Roberts 2019 for the kinetic-energy choice under heavy tails; Girolami and Calderhead 2011 for the generalised implicit leapfrog; Randers 1941 for the asymmetric metric), with an organic-critical reading of the user's geometry document (section 12.3):
Bounded velocity by construction. The velocity
grad_p K = c M^{-1} p / sqrt(p' M^{-1} p + m^2 c^2) has M-norm strictly
below c for every momentum — the property that domps the tails. The
0.5 log det M term is the same normaliser as the Gaussian Riemannian
kinetic, so the theta-marginal of the joint is exactly the posterior: the
kinetic energy is a preconditioner, not part of the target, and the
Metropolis correction with the exact density keeps the sampler exact for any
speed and rest_mass (they govern only efficiency).
Strict generalisation of the Riemannian level. As c -> infinity the
kinetic energy reduces to the Gaussian Riemannian kinetic of RG.3 (the
non-relativistic limit); larger speed mixes faster in the bulk with less
tail-taming, smaller speed caps sooner and is more robust.
Finsler structure, with the asymmetric Randers piece deferred on purpose.
A relativistic kinetic energy is the Legendre dual of a Finsler norm on
velocities — a norm not induced by an inner product, the speed c being the
Finsler unit ball. The asymmetric Randers extension F = sqrt(g(v,v)) + beta(v) (an irreversible drift) is not included: it makes the kinetic
energy odd in p and would break the reversibility the Metropolis
correction relies on. It models irreversible dynamics — the document's
ontological interest, the cost of raising complexity differing from lowering
it — a different goal than exact sampling of a fixed posterior. The even,
bounded relativistic kinetic realises exactly the part of the Finsler
insight that serves robust, exact sampling.
Non-separable Hamiltonian -> dedicated integrator. Because K depends on
both theta (through M) and p, the Hamiltonian is non-separable: the
explicit leapfrog does not apply, and the existing implicit leapfrog
hardwires the Gaussian velocity. A dedicated generalised implicit leapfrog
(three reversible, volume-preserving sub-steps with the relativistic
velocity) is carried in the metric's integrator slot, run by
gdpar_geom_hmc() in place of the default leapfrog; the default branch stays
bit-identical. The momentum is refreshed from the exact relativistic
momentum law (an inverse-CDF radial sampler under p = L q).
Opt-in throughout: the default branch is bit-identical, no fit path is touched
and the golden regression is intact. gdpar_geom_hmc() now uses a metric's
own kinetic energy when it carries one (and rejects, rather than crashes on, a
proposal that drives the kinetic or target undefined). Validated in three
layers — algebraic (kinetic gradients vs finite differences, the M-norm
velocity bound, the Gaussian limit, the dedicated integrator's reversibility
and energy conservation, the momentum law); statistical ungated (a heavy-tailed
Student-t recovered); statistical gated (a correlated multivariate-t where the
Euclidean kinetic overshoots, and a real cmdstan-backed correlated
multivariate-t with a position-dependent SoftAbs mass matching the no-U-turn
sampler). Decision D94.
Sub-Riemannian metric and integrator (gdpar_geom_metric_subriemannian()).
The level-5 geometry of the Block RG hierarchy and the remedy for a
quasi-deterministic posterior (the eBird count / tweedie case), where the
typical set contracts onto a lower-dimensional manifold: the expected Fisher
information grows without bound along the stiff "wall" directions while the
"floor" directions still carry genuine variation. A standard or even
Riemannian sampler is then forced to a vanishing step by the walls; this
geometry equips only a distribution of accessible directions with motion and
glides along the floor instead.
Methodological decisions of this sub-phase (Montgomery 2002; Shahbaba, Lan, Johnson and Neal 2014 for the split):
gdpar_geom_fisher_simulator()) is
eigendecomposed at a reference; a continuous verticality filter
w_i = sigma((log lambda_i - log tau) / softness) blends each direction
between floor (small lambda) and wall (large lambda) with no hard cut, so
a borderline direction is split smoothly and reversibility is never broken by
a discontinuous classification.A = U diag(w_i lambda_i) U^T defines a fixed reference quadratic whose
harmonic flow is integrated in closed form (a symplectic rotation per mode);
the stiff walls cost no step size while the floor follows the free-drift
limit. A Strang splitting (half residual kick, exact flow, half kick) keeps
the scheme symplectic and time-reversible. The default threshold tau is the
floor scale, which caps the leapfrog residual at the floor so the step is
governed by the floor and never the walls; suggested_epsilon exposes that
step. The Metropolis correction with the exact density keeps the sampler
exact however coarse the Gaussian treatment of the walls is — only
efficiency, never correctness, depends on it.Validated in three layers: algebraic (the exact flow conserves the reference
quadratic and reverses to machine precision; the Strang trajectory is
reversible; exact on a pure quadratic at any step; the filter is monotone and
caps the residual); statistical (a mild canyon recovered through
gdpar_geom_hmc); and gated heavier runs — the G4 canyon at n = 1000 (the
sub-Riemannian sampler accepts 0.90 with zero divergences and recovers both
scales where the Euclidean sampler diverges on every proposal at the same
step) and a real cmdstan-backed near-deterministic Poisson count (the
simulation-based Fisher matches the analytic information to 2%, and the
sampler matches NUTS in mean and standard deviation at a step that wrecks the
Euclidean sampler) — the structural rehearsal for the eBird tweedie of RG.7.
Hardened proposal evaluation in gdpar_geom_hmc(). A proposal that drives
the target undefined (a cmdstan log-density that throws on a non-finite
unconstrained value, an overflowing gradient) is now caught and counted as
divergent rather than crashing the run. Stable proposals never trigger the
handler, so the default branch stays bit-identical.
Simulation-based estimator of the expected Fisher information
(gdpar_geom_fisher_simulator()). Completes the learned Riemannian metric
for models with no closed-form Fisher: the expected Fisher
I(theta) = E_y[ s s^T ] is estimated by the average outer product of the
log-likelihood score over data sets simulated from the model at theta,
positive semi-definite by construction and unbiased. The estimate is a
deterministic, reproducible function of theta (the RNG is reseeded from a
position key and restored). The SoftAbs mean of the Gaussian-process surrogate
acts as a structural control variate, so few replicates per site suffice. It
plugs into the fisher slot of gdpar_geom_metric_gp_fisher(). Validated
against closed forms: a constant Fisher (a normal location model) and a
position-dependent one (a Poisson log-rate model), to a few percent.
Generative targets. gdpar_geom_target() gains optional simulate(theta)
and score(theta, y) slots (the per-data-set log-likelihood score whose outer
product is the Fisher), distinct from the fixed-data grad_log_prob. Default
branch unchanged.
Adaptive Riemannian HMC with online active learning
(gdpar_geom_rmhmc_adaptive()). The full novelty-driven loop: the learned
metric is held fixed within each trajectory (preserving reversibility) and
re-learned only between rounds from a reservoir that grows where the
surrogate's epistemic novelty is high, with a decreasing (Robbins--Monro-style)
number of admitted sites per round so the metric sequence settles; a final
sampling phase uses the frozen metric. The sampler is exact in every phase
regardless of the metric (the metric is a preconditioner, not part of the
target), so no delayed acceptance is needed; only the efficiency is heuristic
and is measured (E-BFMI, acceptance, novelty trace), not asserted (ORPHEUS
section 16.3). A gated end-to-end demonstration runs the full pipeline over a
real cmdstan-backed Poisson count model — the structural rehearsal for the
eBird Tweedie of RG.7.
gdpar_geom_metric_gp_fisher()). The general realisation of the natural
Rao–Amari Riemannian metric where no closed form exists, building
M(theta) = L(theta) L(theta)^T from a Gaussian-process surrogate of the
log-Cholesky factor. Following an organic-critical reading of ORPHEUS-PIMC-A
(the new sections 8 and 16), the surrogate's mean function is the SoftAbs
curvature of Capa 1 and the process learns only the smooth residual to
the expected Fisher at the reservoir sites. A single object unifies three
ORPHEUS components without their drawbacks: (i) the metric degrades
continuously to the always-available SoftAbs far from the reservoir (the
kernel decays, the posterior mean returns to its SoftAbs mean) — no hard
metric switch that would break reversibility; (ii) the predictive variance is
the principled epistemic-novelty detector (novelty(theta)); (iii)
positive-definiteness, smoothness and an analytically differentiable
dmass(theta) hold by construction. The sampler stays exact for any
surrogate quality (the metric is a preconditioner, not part of the target;
the Metropolis correction with the exact density is the corrector — no delayed
acceptance is needed). Importance weighting (weights, 1/Q) and the
closed-form-or-simulated Fisher slot (fisher) are exposed for the general
case; gdpar_geom_reservoir() collects phase-one sites from a warmup run.Var(v)).gdpar_geom_metric_riemannian()). The first level-3 metric of the Block RG
hierarchy, extending the pluggable interface of the RG.2 engine with
position_dependent = TRUE. Two curvature sources: the expected Fisher
information (curvature = "fisher", the natural metric of the statistical
manifold, positive-definite by construction and supplied as a function), the
primary maximally robust choice; and the SoftAbs regularisation of the
observed Hessian (curvature = "softabs", Betancourt 2013, eigenvalues mapped
to lambda * coth(alpha * lambda)), fully general and used as the cold-start
and extrapolation fallback. The metric exposes the spatial derivatives
dmass(theta) (the SoftAbs derivative via the Daleckii–Krein formula; the
Fisher derivative analytic or finite-differenced).gdpar_geom_hmc() now records the energy trace and the E-BFMI
(energy Bayesian fraction of missing information) and routes position-dependent
metrics through this integrator automatically.Var(v), where the Euclidean metric
stalls; cmdstan-gated: Stan's $hessian matches the analytic Hessian of the
suite targets and yields the same SoftAbs metric. The Riemannian sampler stays
exact for any metric (the metric is a preconditioner, not part of the target;
the Metropolis correction with the exact log-density is the corrector — no
delayed acceptance is needed for the metric). Opt-in: the default fit path is
untouched, goldens stay bit-identical.gdpar_geom_hmc(),
gdpar_geom_target(), gdpar_geom_metric_euclidean()). A pure-R
Hamiltonian integrator that delegates the log-density, gradient and Hessian to
a compiled backend — a cmdstan model built with
compile_model_methods = TRUE (exposing $log_prob / $grad_log_prob /
$hessian) or an R closure such as the dual targets of
gdpar_geometry_suite(). This is the canonical decision-A scaffolding (D88)
for the Block RG geometry hierarchy: pluggable metric / kinetic-energy /
symplectic-integrator interfaces with the Euclidean level delivered now
(constant metric, standard leapfrog, static HMC with a Metropolis
correction). The Riemannian (Fisher / SoftAbs, RG.3), Finsler / relativistic
and sub-Riemannian levels extend the same interfaces. Validated in three
rigour layers: a deterministic check against the exact Hamiltonian flow of
a Gaussian (energy conservation, reversibility, second-order accuracy), a
statistical check that HMC recovers the easy targets' moments, and a
cmdstan-gated algebraic check that the R closure gradient matches Stan's
$grad_log_prob and that a cmdstan-backed target drives the same trajectory.
Opt-in and standalone: the default fit path is untouched, goldens stay
bit-identical.inst/benchmarks/results/block_rg_calibration.md: six of the eight pathology
classes are recovered at 0.93–1.00, while funnel and heavy tails remain
mutually confusable (recall ~0.6–0.7) and the deepest funnels and heaviest
tails are the hardest cases — gdpar's diagnostic can and does miss them.gdpar_geometry_thresholds() defaults. Eight thresholds
moved to catch the signals as the short diagnostic pilots actually attenuate
them (e.g. condition_high 50→12, heavy_kurtosis_high 3→1.8,
boundary_prox_high 0.10→0.02, nslope_grows 0.30→0.80). The calibration is
against the synthetic suite, not a claim of optimality on real posteriors; the
helper stays exposed for re-calibration.gdpar_geometry_diagnostic(). Probes the
geometry of a posterior with cheap Hamiltonian pilots and classifies the
sampling pathology, localises the culprit parameter(s), estimates the
difficulty-vs-n behaviour, recommends the geometry level that remedies it,
and estimates cost. It is the diagnostic half of the Block RG capability
(geometric robustness of sampling), motivated by the eBird Tweedie count
coordinate that the no-U-turn sampler cannot traverse (session B9.21). The
default fit path is untouched: this is a standalone, opt-in tool, so the
goldens stay bit-identical.gdpar_geometry_suite(): eight synthetic geometries of known
difficulty. A falsifiable calibration backbone covering the full pathology
taxonomy (isotropic control, anisotropic, Neal's funnel, heavy tails,
quasi-deterministic canyon, multimodal, boundary-pegged, flat direction).
Each target is dual (a Stan program plus an R log-density closure with an
analytic gradient on the unconstrained scale), so the same geometry can be
exercised by the cmdstan pilots and, later, by the R-native geometric engine
of RG.2. Analytic gradients are cross-checked against finite differences to
better than 1e-9.gdpar_geometry_thresholds() exposes the classifier thresholds as
tunable data for the calibration sweep.print methods). No existing fit path,
golden, or vignette changed.FAIL 0 | WARN 8 | SKIP 0 | PASS 76): the
statistical fit-fuzz invariances (INV-C: permutation invariance and
Lipschitz continuity of the EB location estimate), the 12-family
outcome-outlier stress (STR-C: every canonical family returns finite
estimates or aborts with a canonical error, never a silent NaN/Inf), and
the full bit-exact compare-path.compute_diagnostics()
computes max(rhat, na.rm = TRUE) over a vector that can be entirely
NA under degenerate outlier-laden data, yielding -Inf plus a warning
for the reported max R-hat. The point estimate stays finite (the contract
holds); the all-NA guard is a documented candidate fix, confined to the
gated stress path, and does not affect the default suite or R CMD check.inst/benchmarks/results/block9_internal.md §11 for the full report.gdpar() and gdpar_eb() now reject Inf / -Inf (as well as
NA / NaN) outcomes pre-fit with a gdpar_input_error. Previously an
infinite outcome slipped past the is.na() guard (is.na(Inf) is
FALSE) and surfaced late as an opaque gdpar_eb_numerical_error. The
broadened, type-aware predicate (!is.finite() for numeric outcomes,
is.na() otherwise) preserves behaviour on non-numeric outcomes.gdpar_dependence_diagnostic() (D75). Quantifies serial
dependence in the residuals of a scalar Empirical-Bayes fit (lag-1
autocorrelation, Durbin-Watson, Ljung-Box). It turns the iid-violation
risk from an invisible theoretical hazard into a measured quantity and
gates the remedy below. Scalar EB path only this session; the spatial
analogue (Moran's I) and the full-Bayes path are deferred.gdpar_dependence_robust() (D75). Dependence-robust standard
errors and percentile intervals via a temporal moving / circular block
bootstrap that refits the model on contiguous-block resamples of the
data. This is the working-independence + robust-variance stance of
Liang & Zeger (1986): the point estimates are unchanged (consistent
under a correct mean structure, not efficient), only the reported
uncertainty is made robust to temporal dependence. The default block
length follows the n^(1/3) rate (Künsch 1989); the Politis-White (2004)
data-driven constant is deferred.test-block9_dependence.R. Stan-free resampler algebra
and API guards run by default (including a bit-exact resampler golden,
D76); end-to-end fit-based behaviour (AR(1) detection, robust-SE table,
seed-determinism) is gated by GDPAR_RUN_BLOCK9_DEP_FITS. The
HMC-dependent bootstrap SEs are not frozen bit-exact (HMC is not
reproducible across Stan toolchains).test-block9_stress.R Section A gains STR-A4 (an infinite outcome
aborts pre-fit with gdpar_input_error, the D74 contract).test-block9_adversarial.R (Layer 1; 11 test_that).
Structural adversarial fuzz over the canonical architecture: aliased
grouping / collinear designs (C4, C4-bis), duplicate
prior_canonical_kind D-ID violations, structurally degenerate priors,
and latent HOM stratification. All aborts are pre-fit
(gdpar_input_error / gdpar_identifiability_error) or passed = FALSE verdicts; no Stan compilation. Documents two validation
boundaries: prior content (negative scale) is a Stan-side reject, and
HOM is a fit-time (not pre-fit) condition.test-block9_invariances.R (Layer 2; 6 test_that, 4
default + 2 gated). Deterministic equivariances of the
identifiability eigenstructure: row permutation and covariate
rescaling / sign flip leave the verdict and conditioning invariant
(exact up to floating-point summation order). Statistical fit-based
invariances (permutation invariance, noise-injection continuity) are
gated by GDPAR_RUN_BLOCK9_FIT_FUZZ (nocturnal B9.8 sub-unit).test-block9_stress.R (Layer 3; 6 test_that, 5
default + 1 gated). NA / missing robustness (NA in AMM covariate, NA
and NaN in outcome all abort with gdpar_input_error) and
extreme-but-finite covariate outliers (99.9% / 0.1% quantiles, 1e6
scale) keep the conditioning finite. The 12-family outcome-outlier
stress sweep is gated by GDPAR_RUN_BLOCK9_STRESS_FITS (B9.8).golden_fb_KxP_gaussian_polynomial_K2_p2, seed 91001) was
reproduced in isolation and matched its frozen RDS bit-exactly
(400 × 815 draws). No drift: the cluster-K4 Stan unification preserved
bit-exactness by construction, so the conditional full reseed (C.ii) is
not triggered.inst/benchmarks/results/block9_internal.md
opened (criterion 1 of the Block 9 closure; updated across B9.7-B9.9)..gdpar_fb_KxP_fit() (R/stan_codegen.R) que compone
.build_amm_design_KxP() + el assembler unificado
.assemble_stan_data_KxP(path = "FB") (M.iv B9.6) + cmdstanr
sampling sobre la piece amm_canonical_pmulti_KxP.stan (J.iv.A
B9.5). Bypassa el guard publico .gdpar_guard_multiparam_multivariate
(que sigue activo para gdpar() user-facing); init conservador
init = 0.5 para evitar rejecting initial values con las familias
Path B (beta / gamma / student_t cerca de los bordes del link).R/eb.R): la signature de
.assemble_stan_data_KxP() ahora acepta path = c("EB", "FB") y
cp_W = FALSE. El branch EB preserva byte-identico el comportamiento
de Sub-fase 8.6.D (use_W = 0 hardcoded por D39 + stan_id en {1, 3}
por D40'); el branch FB lifta ambas restricciones y emite los
campos W esperados por la KxP piece (X, dim_W, d,
W_per_kj_dim, W_type_id, W_n_knots_full, W_knots_full,
W_degree). Backward-compatible por default path = "EB".data-raw/bootstrap_fb_goldens_KxP.R
reemplaza los 4 fit closures stub (D69 canonizada en B9.5) por
closures reales que invocan .gdpar_fb_KxP_fit() con iter_warmup = 200, iter_sampling = 200, chains = 2, sobre datos sinteticos
pequenos (n = 80-100) por seed canonical:
91001 gaussian K=2 p=2 / 91002 beta K=2 p=2 / 91003 gamma K=2 p=3 /
91004 student_t K=3 p=2. Persiste un snapshot por golden con
draws_matrix + stan_data + canonical metadata bajo
tests/testthat/data/golden_fb_KxP_<family>_polynomial_K<K>_p<p>.rds
y appendea 4 filas al manifest CSV con sub_phase = "9.3.d" +
status = "fitable".test_that blocks en
tests/testthat/test-stan_codegen_canonical.R (frozen md5 de los
6 EB pieces relocadas + CANONICAL_HELPERS placeholder sanity +
legacy top-level removidos + KxP EB intactos + dispatcher render
sin leftover placeholders + assembler EB-side preservado +
assembler FB-side acepta Path B family set + driver guards).
Section 4 nueva en tests/testthat/test-fb_goldens_9_3_d_KxP.R:
fitable RDS structure + draws_matrix has theta_ref_kp + stan_data
carries FB extension fields (M.iv).Cascade L heredada cerrada operativamente (canonizada en B9.5
bajo decision L.iv.A para amm_distrib_K, ejecutada en B9.6 bajo
O.a para los 6 EB-side templates restantes con helpers inline). Los
6 EB-side templates con inline bspline_basis_eval +
apply_W_basis_diff fueron relocados a
inst/stan/_canonical_pieces/ con el placeholder
// {{CANONICAL_HELPERS}} reutilizando la infraestructura G.iv de
B9.4:
amm_canonical_eb_marginal.stan (md5 7a279bf087a0fd0a95d35a490814673f, 6442 bytes)amm_canonical_eb_marginal_multi.stan (md5 9250144e95df5fd409cb5cc9df3f3730, 13661 bytes)amm_canonical_eb_marginal_K.stan (md5 dca0b8282e4c617d259e8f17d1f95be5, 36709 bytes)amm_canonical_eb_conditional.stan (md5 8321f9d7d99d0c702045a2886cd1f98c, 6984 bytes)amm_canonical_eb_conditional_multi.stan (md5 186b6ac9f9810584ea21522974fe641d, 13878 bytes)amm_canonical_eb_conditional_K.stan (md5 169274eca2a64e2a0d5a060e2d0c3033, 36250 bytes)Bit-exact Stan semantics preservadas token-wise (6/6); solo
difieren comments del helpers block, unified al inyectar la piece
canonica amm_canonical_helpers.stan via el dispatcher.
Callers actualizados: generate_stan_code_multi() +
generate_stan_code_K() extendidos con dos switch entries cada
uno traduciendo el legacy nombre al canonical (mirror del patron
L.iv.A B9.5). .gdpar_eb_render_template() extendido con
translation + helpers substitution (mirror del dispatcher canonical).
6 legacy top-level eliminados ([[feedback_depuracion_obsoletos]]):
amm_eb_marginal{,_multi,_K}.stan + amm_eb_conditional{,_multi,_K}.stan.
Los 2 KxP EB templates (amm_eb_marginal_KxP.stan +
amm_eb_conditional_KxP.stan) permanecen en inst/stan/ root
intactos: NO tienen helpers inline porque la canonization 8.6.D
hardcodea use_W = 0 per D39.
Cluster K4 cerrado completamente (FB + EB cascade + fit harness): 4/4 deudas sustantivamente cerradas + L cascade propagada + D69 operativamente cerrada.
New canonical piece
inst/stan/_canonical_pieces/amm_canonical_pmulti_KxP.stan
(~740 lines): Full-Bayes multi-parametric multivariate template for
the K >= 2 K-individual slots crossed with p >= 2 AMM coordinates
regime, the bit-exact FB counterpart of the EB Path C templates
canonized in Sub-fase 8.6.D. Mirrors amm_eb_marginal_KxP.stan
extended with the globally shared W modulating component and the
full Path B family set
{1 gaussian, 3 neg_binomial_2, 5 beta, 6 gamma, 7 lognormal_loc_scale, 8 student_t, 9 tweedie, 10 zip, 11 zinb, 12 hurdle_poisson, 13 hurdle_neg_binomial_2}. Canonized under
decision I.iv lateral (DESIGN_9_3_D_PATH_C.md APERTURA cerrada
B9.4) + J.iv.A piece arquitectonica (sub-decision 3.2 dedicated
piece with canonical helpers reuse via // {{CANONICAL_HELPERS}}
placeholder consumed by the dispatcher per the G.iv pattern of
B9.4).
Dispatcher extension: .gdpar_emit_canonical_stan()
(R/stan_codegen.R) now accepts
p_class = "pmulti_KxP" (plus "distrib_K" from the L.iv.A
cascade below). The KxP branch consumes the canonical priors +
{{THETA_REF_PRIOR_BLOCK}} (per-family anchor block via the new
helper .gdpar_build_theta_ref_prior_block_KxP(); Tweedie K=3
carries the bounded p slot per coord via the iterated slice form)
{{W_SCALE}} / {{W_PRIOR}} for the CP / NCP toggle. The KxP
piece emits NCP-only a/b blocks inline (no
DATA_CP_A_PER_K_DECL / TP_A_BLOCK / MODEL_A_BLOCK overrides
for B9.5 atomic scope; per-slot per-coord CP / mixed variants are
deuda for B9.6+).Thin wrapper: new internal generate_stan_code_KxP_FB()
preserves the canonical source-of-truth pattern of B9.3; the
user-facing guard .gdpar_guard_multiparam_multivariate still
aborts gdpar() on K > 1 + p > 1, deferring the public surface
lift to a future sub-phase.
Reserved canonical seeds for the bootstrap goldens roster
(DESIGN_9_3_D_PATH_C.md §4): 91001..91004 cover the K.c roster
minimum (gaussian K=2 p=2, beta K=2 p=2, gamma K=2 p=3,
student_t K=3 p=2). 91005..91099 stay reserved for the
expanded roster diferido a B9.6.
Bootstrap script data-raw/bootstrap_fb_goldens_KxP.R
(env-var gated on GDPAR_BOOTSTRAP_FB_GOLDENS +
GDPAR_BOOTSTRAP_FB_GOLDENS_KxP, mirroring the 8.6.D pattern):
records 4 metadata-only manifest rows for the reserved seeds. The
actual cmdstanr sampling is deuda D69 canonized in B9.5 +
deferred to B9.6 under [[feedback_token_economy_session_planning]];
the B9.5 cycle closed the codegen + dispatcher + piece architecture
with bit-exact tests over the canonical Stan source (frozen md5
invariants).
New roster test
tests/testthat/test-fb_goldens_9_3_d_KxP.R: 8 codegen tests +
3 RDS-presence tests (skip_if RDS not present + skip_if manifest
not present, the canonical 8.6.D convention). Always-on codegen
tests verify the family-specific dispatch + the seed range
reservation + the bootstrap script env-gating + the determinism
of the wrapper.
Test extensions:
tests/testthat/test-stan_codegen_canonical.R adds frozen md5
invariants for the new KxP piece + the relocated distrib_K piece +
the dispatcher OUTPUT for gaussian / tweedie + helpers presence
sanity + cmdstanr parse-check (skip_if cmdstan unavailable).
Caveat heredado Path B logit-strict bajo Path C registered in the new piece header and in DESIGN_9_3_D_PATH_C.md §5: the K.c B9.5 roster (gaussian + beta + gamma + student_t) avoids logit-strict puro; roster expansion to {10 ZIP, 12 hurdle_poisson} (link logit in pi) reserved for B9.6 with skip_if granular if emergent numerical instability appears.
inst/stan/_canonical_pieces/amm_canonical_distrib_K.stan
(851 lines): byte-relocated copy of the legacy
inst/stan/amm_distrib_K.stan with the inline
bspline_basis_eval + apply_W_basis_diff definitions replaced
by the // {{CANONICAL_HELPERS}} placeholder consumed by the
dispatcher per the G.iv pattern of B9.4. The Tweedie family
helpers (tweedie_log_W_series, tweedie_log_f_series,
tweedie_log_f_saddlepoint, tweedie_lpdf, tweedie_rng) and
apply_inv_link_by_id stay inline because they are K-specific
(not shared across templates)..gdpar_emit_canonical_stan() accepts
p_class = "distrib_K" with the full K-style placeholder set
(THETA_REF_PRIOR_BLOCK + DATA_CP_A_PER_K_DECL + TP_A_BLOCK +
MODEL_A_BLOCK + canonical priors + W toggles). The
template_name override supports the EB-side K templates
(amm_eb_marginal_K.stan, amm_eb_conditional_K.stan per D34 of
Sub-fase 8.6.C); those EB templates retain inline helper copies
pending the EB cascade in B9.6 (deuda L heredada bajo L.iv.A
symmetry).generate_stan_code_K() refactor: now a thin wrapper that
delegates to the canonical dispatcher with p_class = "distrib_K",
preserving signatures and tests bit-for-bit. The legacy
template_name = "amm_distrib_K.stan" default is translated
internally to "amm_canonical_distrib_K.stan"; the EB-side K
templates pass through unchanged.inst/stan/amm_distrib_K.stan removed per
[[feedback_depuracion_obsoletos]].generate_stan_code_K(prior, family = ...) are preserved bit-exact
across L.iv.A (only the unified helper-block comments differ vs the
legacy template). The 14 K=2 / K=3 goldens in
tests/testthat/data/golden_K*.rds remain valid by construction
(they fit data, not source hashes); cmdstanr cache recompiles
once on first call post-B9.5 and is then warm.amm_eb_marginal*.stan, amm_eb_conditional*.stan for
single / multi / K variants) retain inline copies pending B9.6
cascade (canonized as deuda colateral L heredada in
CHARTER_BLOQUE_9.md §2.1).gdpar_ksd_joint(eb_fit, fb_fit, ...)
(R/ksd_joint.R): operationalizes the open question documented in
the Roxygen of gdpar_compare_eb_fb() that the marginal
total-variation distance is only a coarse proxy and the joint
posterior discrepancy deserves a density-free spectral metric.
Returns an object of class gdpar_ksd_joint with the KSD
V-statistic, kernel/bandwidth configuration, target empirical
Gaussian (mean + covariance), and call. Canonized under decision
H.iv lateral: IMQ base kernel default (Gorham-Mackey 2017,
dim-independent rate for log-concave targets) + median heuristic
bandwidth + ESS-weighted variant via posterior::ess_basic. RBF
base kernel is provided as a textbook alternative.grad_log_prob() is
documented as a Block 9.x extension.print.gdpar_ksd_joint, summary.gdpar_ksd_joint,
print.summary.gdpar_ksd_joint..gdpar_ksd_safe_min_ess
was updated to use the canonical posterior::summarise_draws() API
(was calling posterior::ess_basic() directly on a draws_matrix,
which is method-dispatched ambiguously across posterior versions and
could return either a numeric scalar or a tibble). Resolves the two
ESS-weighted-related errors observed in the full-suite run.v07b_eb_multivariate.Rmd ("Joint Kernel Stein Discrepancy: the
gdpar_ksd_joint() Helper") with subsections on the Stein
operator, the empirical Gaussian target, the base kernel and
bandwidth, the ESS-weighted variant, and the complementarity
between gdpar_compare_eb_fb (marginal TV) and gdpar_ksd_joint
(joint KSD). Existing Section 11 (Summary) renumbered to Section
12.tests/testthat/test-ksd_joint.R): 14 test_that blocks
/ 45 assertions covering input validation (5), happy path and
class invariants (3), statistical layer (KSD small under matched
Gaussians + KSD detects shifted Gaussians) (2), ESS-weighted
variant (2), S3 methods (2). Mock fits avoid Stan dependence.inst/stan/_canonical_pieces/amm_canonical_helpers.stan
(new file): dedicated Stan source FRAGMENT containing the two
canonical helpers bspline_basis_eval and apply_W_basis_diff.
Closes deuda f of 9.3.a (helpers duplication across the K=1 FB
pieces) under canonized decision G.iv lateral (R-side textual
substitution). The piece has NO surrounding functions { } block;
it is intentionally not standalone-parseable by cmdstanr::stanc
and is inserted by .gdpar_emit_canonical_stan() at the
// {{CANONICAL_HELPERS}} placeholder inside the body piece's
functions { } block.amm_canonical_p1.stan (13878 → 10722
bytes, md5 f899811e... → 730c93f5...) and
amm_canonical_pmulti.stan (16028 → 13764 bytes, md5
112a4e68... → 01fc4b04...) now contain only the
// {{CANONICAL_HELPERS}} placeholder inside their functions { }
block; helper definitions removed inline. Total bytes saved across
the K=1 FB pieces: 2240 (~7.5%); future EB-side cascade in B9.5+
will reuse the same helpers piece with additional savings..gdpar_emit_canonical_stan() reads the
helpers piece and substitutes the placeholder via gsub() before
applying the other prior / parametrization placeholders. The gate
is grepl("// {{CANONICAL_HELPERS}}", src, fixed = TRUE), so
templates that do not contain the placeholder are unaffected
(forward-compatible with the EB cascade and with any future piece
that ships its own helpers).8930c724... for NCP and 1e6d7d28... for MLE):
the helpers in the dedicated piece are verbatim from the
pre-G.iv amm_canonical_p1.stan. Output md5 changes for the
pmulti path (cf182996... → 2ede60f6..., +892 bytes) due to
unified helper comments (pre-G.iv amm_canonical_pmulti.stan
carried a brief 4-line block referencing the then-open deuda f;
it is replaced by the detailed 16-line algorithmic comment block).
Stan semantics preserved bit-exact (only comments differ);
goldens unaffected (they fit data, not source hashes);
cmdstanr::cmdstan_model() cache hash recompiles once on first
call post-B9.4.tests/testthat/test-stan_codegen_canonical.R): 5 new
test_that blocks (frozen md5 of helpers piece + frozen output
md5 invariants for p1 NCP / p1 MLE / pmulti default + sanity
check that helpers are inserted + defensive smoke). The 2 existing
frozen-piece-md5 tests updated with the new post-G.iv hashes..gdpar_emit_canonical_stan() (new internal
helper, R/stan_codegen.R): single source-of-truth dispatcher for
the (p, K=1) FB Stan templates. Replaces the duplicated
template-read plus substitute logic that previously lived inline in
generate_stan_code() (p = 1) and generate_stan_code_multi()
(p >= 1). Reads a canonical piece selected by spec$p_class and
applies the appropriate placeholder substitution. The dispatcher
emits the same substituted Stan strings as the legacy paths for the
same inputs.inst/stan/_canonical_pieces/ (new
directory): the legacy inst/stan/amm_main.stan and
inst/stan/amm_distrib_multi.stan are renamed to
amm_canonical_p1.stan and amm_canonical_pmulti.stan respectively
and moved into inst/stan/_canonical_pieces/. The relocation is
byte-identical (md5sum preserved); bit-exactness of the substituted
Stan source is preserved by construction; no golden re-bootstrap
required.inst/stan/amm_main.stan
and inst/stan/amm_distrib_multi.stan are deleted from the package
root. Backward-compat for callers of generate_stan_code_multi()
that pass template_name = "amm_distrib_multi.stan" is preserved
via wrapper-side translation to the canonical piece name; other
template_name values (the EB-side multi templates pending the
cascade in B9.5+: amm_eb_marginal_multi.stan,
amm_eb_conditional_multi.stan) are served from
inst/stan/ root unchanged.generate_stan_code() and
generate_stan_code_multi() become thin wrappers (~10 lines each)
that construct a canonical codegen spec and delegate to
.gdpar_emit_canonical_stan(). Signatures preserved bit-for-bit;
no caller migration required across the package. The K-individual
generator generate_stan_code_K() (K >= 2, p = 1, multi-parametric
via amm_distrib_K.stan) is out of scope for 9.3.a and is left
untouched.tests/testthat/test-stan_codegen_canonical.R,
9 tests): frozen md5 hashes of the two canonical pieces;
determinism of the dispatcher across calls (both p1 and pmulti
paths); equivalence of the public wrappers with the direct
dispatcher call; input validation (invalid p_class, malformed
cp_a_per_k); cmdstanr parse-check of the substituted output for
both pieces (skipped when cmdstan is unavailable).tests/testthat/test-stan_codegen_multi.R): the
legacy existence check for amm_distrib_multi.stan is canonized to
verify the relocation (canonical pieces exist; legacy top-level
templates are removed)._canonical_pieces/), A3.1 (snapshot frozen + md5 hash dual
verification), B4 (FB cluster first, EB cascade in B9.5+), B5
(deletion during B9.3 once tests pass), B6 (9.3.c intercalable in
B9.4), C7.2 (deuda f helpers dedup deferred to B9.4 under the
token-economy principle), D8 (rename without lifecycle since
functions are internal), E9 (Roxygen + this NEWS entry; vignette
section in v07_path_1.Rmd), F10 (closure criterion under rigor
de tres capas), F11 (atomic sub-unit pacing per
feedback_token_economy_session_planning) canonized concurrently.bspline_basis_eval and apply_W_basis_diff
duplicated in pieces) remains open as planned for B9.4 under the
token-economy session-planning principle; #include mechanism to
factor the helpers into _canonical_pieces/amm_canonical_helpers.stanfunctions
will be canonized in the next session.CHARTER_SUBFASE_8_6.md, component 2): the single-step ridge of
Sub-phase 8.6.B is extended into an adaptive geometric loop bounded
by laplace_control$ridge_max_iter (default 10) and multiplied by
laplace_control$ridge_grow_factor (default 10.0) at each iteration
until the post-ridge condition number is at or below
laplace_control$kappa_threshold or the loop is exhausted.|det(cov)| < epsilon_lm trigger (new condition): the helper
also triggers when the determinant of the empirical Laplace
covariance is strictly below laplace_control$epsilon_lm
(default sqrt(.Machine$double.eps), approximately 1.5e-8) even
when the matrix is technically positive-definite, matching the
canonical Section 2.8 wording \emph{|det(H)| < epsilon_LM} of the
ridge trigger..gdpar_eb_lm_perturb() (new internal helper, R/eb.R):
factored out from the per-Path inline ridge of Sub-phase 8.6.B/C/D,
shared by both .gdpar_eb_maximize_marginal() (Paths A/B) and
.gdpar_eb_maximize_marginal_KxP() (Path C, per-slot). Returns a
list with cov_perturbed, lambda_used, n_iter, kappa_post
and status (one of "not_needed", "converged", "exhausted").diagnostics_numerical fields: lm_n_iter, lm_status,
kappa_post_ridge (Paths A/B); lm_lambda_per_slot,
lm_n_iter_per_slot, lm_status_per_slot (Path C). The pre-9.3.b
lm_perturbation and kappa_per_slot fields are preserved.gdpar_eb_numerical_error messages include the L-M
status, lambda used and number of iterations for both Paths A/B and
Path C, enabling downstream code to disambiguate "exhausted" from
"converged but still above threshold".tests/testthat/test-eb_numerical_robustness.R exercising the new
helper across well-conditioned, rank-deficient, small-determinant,
exhausted-by-budget and 1x1 covariances; plus 4 new resolver tests
for the three new laplace_control defaults and their validators.gdpar_compare_eb_fb(eb_fit, fb_fit, level = 0.95, tv_bins = 30L)
(new function in R/compare_eb_fb.R): canonizes the operational
EB-vs-FB comparison promised in v07 Section 11.1. Takes a
gdpar_eb_fit and a gdpar_fit fitted on the same dataset and
returns a gdpar_eb_fb_comparison object with three tables:
theta_diff_table: per-anchor cell comparison of
$\widehat\theta_{\text{ref}}^{\text{EB}}$ vs the FB posterior mean
$E_{\text{FB}}[\theta_{\text{ref}}]$, with the difference and the
difference normalized by the FB standard error. Operationally
verifies Theorem 7A / 7A* (first-order asymptotic equivalence).tv_table: marginal empirical total variation distance per common
$\xi$ parameter via histogram plug-in. Theorem 7A predicts marginal
TV $\to 0$ in probability as $n \to \infty$.coverage_table: per anchor cell, EB credible-interval width
(with inflation applied when eb_correction = TRUE), FB
credible-interval width, and the ratio
width_eb / width_fb. Operationally verifies the $O(n^{-1})$
under-cover prediction of Proposition 7B / 7B* matricial /
7B* tensorial.print.gdpar_eb_fb_comparison,
summary.gdpar_eb_fb_comparison,
print.summary.gdpar_eb_fb_comparison in
R/compare_eb_fb_methods.R. The summary aggregates the TV table
and the coverage table into quartile / min / median / max / mean
summaries with the full per-anchor diff table.vignettes/vop07_eb_workflow.Rmd (new operational vignette,
~400 lines): walks through gdpar_eb() end to end across the four
path regimes, the numerical diagnostics ($\kappa(H)$, LM ridge,
multi-start dispersion, kappa_per_slot and slot dispositions under
Path C), the EB-vs-FB comparison via gdpar_compare_eb_fb(), and
the troubleshooting recipes for ill-conditioned Hessians and
conditional-HMC instability under logit-strict links (Path B
caveat heredado de 8.6.C + Path C caveat heredado de 8.6.D).
Chunks default to eval = FALSE because each fit compiles Stan
models and takes minutes; re-enable per-chunk or globally to
reproduce the runs locally.vignettes/v07_eb_vs_fb.Rmd Section 11.1 rewritten in
descriptive voice apuntando a gdpar_eb() con las cuatro
variantes de path Stan templates y a gdpar_compare_eb_fb(); se
elimina el contrato textual prospectivo eb_for_theta_ref = TRUE
(nunca implementado per Charter Decision 2.2). Section 13
(«Connections to Subsequent Blocks») extendida con dos viñetas
apuntando a v07b_eb_multivariate (extensión teórica multivariada)
y vop07_eb_workflow (operativa).tests/testthat/test-compare_eb_fb.R (new): unit tests for the
comparator with 4 input-validation guards (unconditional) + 1
end-to-end smoke gated by GDPAR_RUN_STAN_SMOKE_EB=1 + 4
internal-helper unit tests for .gdpar_eb_fb_tv_table /
.gdpar_eb_fb_extract_xi_draws.NAMESPACE: exports gdpar_compare_eb_fb + 3 S3method
registrations regenerated via devtools::document().man/: gdpar_compare_eb_fb.Rd +
print.gdpar_eb_fb_comparison.Rd +
summary.gdpar_eb_fb_comparison.Rd +
print.summary.gdpar_eb_fb_comparison.Rd regenerated.This closure aggregates Sub-fase 8.6 entera (A theoretical extension
CHARTER_SUBFASE_8_6.md; for the closure summary across the five
sub-sub-fases see HANDOFF_SUBFASE_8_6_CIERRE.md.gdpar_eb() is now operative.
The path inherits two new Stan templates canonized in Sub-phase 8.6.D
under decision D36 = (alpha) (Session 13a 2026-05-25):
inst/stan/amm_eb_marginal_KxP.stan and
inst/stan/amm_eb_conditional_KxP.stan. The canonical Path C
architecture is the composition of Path A multivariate ($p > 1$,
Theorem 7A* of v07b) with Path B multi-parametric ($K > 1$, Theorem
7C* of v07b) under decision D38'' = (h) (Session 13b 2026-05-25):
outcome y[n, p] matrix-column shared across the K slots, with
per-slot per-coord linear predictors eta_kp[k, i, j] for the
canonical Path B K-slot multi-parametric family..gdpar_eb_correction_tensor() returns a 3D array
$C^{\text{tensor}}{g, \alpha}[k, , ] = \kappa(\alpha) \cdot
\Sigma^{\text{marg}}{\theta_{\text{ref}, k}}$
retaining all cross-coordinate terms within each slot. The
block-diagonal structure across slots is the canonical factorization
assumption $\pi_\xi = \prod_k \pi_{\xi_k}$. Reduces algebraically to
the Path A matrix correction at $K = 1$ and to the scalar Path B
correction at $p = 1$. The correction_tensor_constant /
correction_tensor_dispositions slots of gdpar_eb_fit report the
tensor and per-slot dispositions ("ok", "non_finite",
"non_psd", "missing", "disabled").stan_id %in% c(1, 3) per slot (Gaussian K=2,
Negative Binomial K=2). The remaining Path B set
${5, 6, 7, 8, 9, 10, 11, 12, 13}$ (Beta / Gamma /
Lognormal-loc-scale / Student-t / Tweedie / ZIP / ZINB /
Hurdle-Poisson / Hurdle-NB) is deferred to a later iteration of 8.6.D
with the explicit numerical caveat of Section 6.1 of the opening
handoff (HMC condicional bajo plug-in EB cerca del borde de soporte
logit/log links + warmup corto). The .gdpar_eb_check_stan_id_for_path
guard rejects deferred families with a clear deferral message.use_W = 0 enforced in
Path C first iteration; the modulating component W is registered as
a follow-on debt against Block 9.x because its canonical
factorization across the $K \times p$ cross structure has no clean
canonization in v07b yet. The orchestrator upstream rejects input
with W != NULL on any slot with gdpar_unsupported_feature_error.a / b coefs are
per-slot per-coord ragged in Path C. Each (slot k, coord j) pair
has its own J_a_per_kp[k, j] and J_b_per_kp[k, j]; the design
matrices Z_a_kp[k, j] and Z_b_kp[k, j] are padded matrices
matrix[n, J_a_max] for maximum statistical generality without
coord-wise restrictions within a slot.theta_ref_kp
$[J_{\text{groups}} \times K \times p]$ anchors followed by
block-diagonal extraction per slot post-hoc to feed
.gdpar_eb_correction_tensor(). Minimum compile count (1 vs K
separate Stan compiles for the per-slot alternative); preserves the
EB-Path-B-Compound canonical estimator (Petrone-Rousseau-Scricciolo
2014). The block-diagonal extraction averages over groups when
$J > 1$ to match the per-slot canonical aggregation of v07b
Section 5.1..gdpar_eb_run_KxP() (R/eb.R,
~350 lines) drives the full pipeline: input validation +
homogeneous-p enforcement, design building via
.build_amm_design_KxP(), stan_data assembly via
.assemble_stan_data_KxP(), conservative path-aware init via
.gdpar_eb_make_random_init_KxP() (Path C-specific because the
default cmdstanr unconstrained random init explodes the K x p joint
via exp(eta_kp[2]) overflow), Laplace via
.gdpar_eb_maximize_marginal_KxP(), conditional HMC sampling, and
tensor correction application. Returns a gdpar_eb_fit with
path = "eb_KxP" and new slots theta_ref_kp_hat
(3D array $[J, K, p]$), theta_ref_kp_se, theta_ref_kp_cov_per_slot
(list of K matrices $p \times p$), correction_tensor_constant,
correction_tensor_dispositions.print.gdpar_eb_fit and summary.gdpar_eb_fit extended with a
Path C branch: report per-slot per-coord theta_ref_kp_hat and SE
for group 1, per-slot kappa diagnostics, and per-slot dispositions.
The summary tabulates the corrected credible intervals per (group, slot, coord) cell with the per-slot inflation factor derived from
the diagonal of the correction tensor block.R/golden_helpers.R: .gdpar_golden_roster_8_6_D() canonizes 8
configurations (4 fitable Gaussian / NB at $p \in {2, 3}$ + 4
metadata-only guard Beta / Gamma / Student-t / ZIP at $p = 2$).
.gdpar_golden_eb_rds_name_8_6_D(family, regime) produces the
canonical filename
golden_eb_<family>_pathC_<regime>.rds so the 8.6.B / 8.6.C / 8.6.D
on-disk namespaces never collide.tests/testthat/test-eb_smoke_by_family.R: 8 new smokes appended
(smokes 35-42). 4 Path C guard smokes (39-42) run unconditionally
(verify gdpar_unsupported_feature_error on Beta / Gamma /
Student-t / ZIP under $K, p > 1$). 4 Path C fitable smokes (35-38;
Gaussian K=2 and NB K=2 crossed with $p \in {2, 3}$) are gated by
GDPAR_RUN_STAN_SMOKE_EB=1.tests/testthat/test-eb_goldens_8_6_D.R: new structural test
mirroring test-eb_goldens_8_6_C.R. Validates roster coverage,
Path C invariants ($K > 1 \wedge p > 1$ in every entry), filename
convention, correction tensor shape under valid PSD inputs and under
injected non-PSD slots, .gdpar_eb_check_stan_id_for_path Path C
acceptance (${1, 3}$) and Path B preservation
(${1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13}$ at $p = 1$), guard class
canonicity, and manifest CSV cross-validation.data-raw/bootstrap_eb_goldens_D.R: new bootstrap script gated
by GDPAR_BOOTSTRAP_EB_GOLDENS=1 + GDPAR_BOOTSTRAP_EB_GOLDENS_D=1.
Bootstraps the 4 fitable Path C goldens with seed 20260526L and 4
metadata-only guard goldens; appends 8 rows to
tests/testthat/data/golden_manifest.csv with sub_phase = "8.6.D".gdpar_eb() are now operative. The two paths inherit four new Stan
templates canonized in Sub-phase 8.6.C under decision D34:
inst/stan/amm_eb_marginal_multi.stan and
inst/stan/amm_eb_conditional_multi.stan (Path A; bit-identical bodies
to amm_distrib_multi.stan modulo the theta_ref_data move + EB
banner) and inst/stan/amm_eb_marginal_K.stan and
inst/stan/amm_eb_conditional_K.stan (Path B; bit-identical bodies to
amm_distrib_K.stan modulo the theta_ref_k_data move + EB banner).
The combined regime $K > 1 \wedge p > 1$ remains guarded with
gdpar_unsupported_feature_error and is queued for Sub-phase 8.6.D..gdpar_eb_generate_stan_marginal() /
.gdpar_eb_generate_stan_conditional(): selects one of the four EB
templates based on the resolved $(K, p)$ pair; the K = 1 + p = 1 leaf
preserves the byte-identical 8.6.B template selection.
.gdpar_eb_check_stan_id_for_path(family, K, p) enforces the
per-path supported stan_id set: Path A supports
${1, 2, 3, 4}$ (inherits amm_distrib_multi.stan), Path B supports
${1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13}$ (inherits
amm_distrib_K.stan's full distributional-regression dispatch).gdpar_eb(): the function now
accepts the same three forms as gdpar() for declaring multi-slot
models: (i) formula as gdpar_formula_set (e.g.,
gdpar_bf(y ~ a(x), phi ~ a(z))); (ii) amm as a named list of
amm_spec; (iii) formula with a()/b()/W() wrappers in its
RHS. When the resolved K is > 1, the new internal orchestrator
.gdpar_eb_run_K() runs Path B end-to-end..gdpar_eb_correction_matrix()): returns
$C^*{g, \alpha} = \kappa(\alpha) \cdot \Sigma^{\text{marg}}{\theta_{\text{ref}}}$
as a $p \times p$ matrix; reduces algebraically to the scalar
Proposition 7B of v07 Section 6 at $p = 1$. The correction_applied
/ eb_correction_constant slot of gdpar_eb_fit now carries either
the scalar form (Sub-phase 8.6.B contract preserved bit-exactly) or
the matrix form (Path A under D34).R/golden_helpers.R: .gdpar_golden_roster_8_6_C() canonizes
17 configurations (10 Path A + 7 Path B; 13 fitable + 4 guard).
.gdpar_golden_eb_rds_name_8_6_C(family, path, regime) produces the
canonical filename
golden_eb_<family>_path<A|B>_<regime>.rds so the 8.6.B and 8.6.C
on-disk namespaces never collide.tests/testthat/test-eb_smoke_by_family.R: 17 new smokes appended
(smokes 18-34). 4 Path A guard smokes run unconditionally (verify
gdpar_unsupported_feature_error on Beta/Gamma/Student-t/Tweedie at
$K = 1, p = 2$). 6 Path A fitable + 7 Path B fitable smokes are
gated by GDPAR_RUN_STAN_SMOKE_EB=1.tests/testthat/test-eb_goldens_8_6_C.R: new structural test
mirroring test-eb_goldens_8_6_B.R. Validates roster coverage,
filename convention, Path A K=1+p>1 / Path B K>1+p=1 invariants,
guard class canonicity, and manifest CSV cross-validation.data-raw/bootstrap_eb_goldens_C.R: new bootstrap script gated
by GDPAR_BOOTSTRAP_EB_GOLDENS=1 + GDPAR_BOOTSTRAP_EB_GOLDENS_C=1.
Bootstraps the 13 fitable goldens with seed 20260525L and 4
metadata-only guard goldens; appends 17 rows to
tests/testthat/data/golden_manifest.csv with sub_phase = "8.6.C".R/eb.R necessary for the multivariate paths to run
end-to-end. (i) Path-aware init dispatch in
.gdpar_eb_maximize_marginal() (Path A / Path B delegate to
cmdstanr's default unconstrained-space random init; K = 1 + p = 1
preserves the explicit .gdpar_eb_make_random_init() bit-exactly).
(ii) Path-aware extraction of theta_ref from Laplace draws:
detects theta_ref_k[j,k] (Path B) and theta_ref[j,k] (Path A)
variable name conventions. (iii) Reshape of theta_ref_data to a
$J_{\text{groups}} \times p$ matrix in the K = 1 pipeline before
the conditional HMC step (the Stan conditional template declares
array[J_groups] vector[p] theta_ref_data and cmdstanr's
auto-packing requires the R-side matrix shape under $p > 1$). All
three are minimum-invasive patches at the R-Stan binding layer;
none touches the canonized Stan templates or the Charter Section 2
decisions.iter_warmup >= 1000, tighter sigma_a_k
prior, or alternative anchor selection. Documented as deferred
numerical caveat for Sub-phase 8.6.D Hessian-fragility expansion.[ FAIL 0 | WARN 0 | SKIP 82 | PASS 3021 ] (+103 PASS / +17
SKIP), of which Sub-phase 8.6.C contributes +63 PASS / +17 SKIP on
top of the +40 PASS canonized by 8.6.C pasos 1-4 (Session 11).gdpar_eb(formula, family, amm, data, ..., eb_correction, laplace_control):
new exported function that fits an AMM canonical model via Empirical Bayes
(Type II ML for theta_ref + conditional HMC for
xi = (a, b, W, sigma_*, phi) given the plug-in estimate
$\widehat\theta_{\mathrm{ref}}^{,\mathrm{EB}}$). Sub-phase 8.6.B is the
base regime $K = 1$, $p = 1$; multivariate $p > 1$ (8.6.C), multi-slot
$K > 1$ (8.6.C), and mixture/hurdle families with min_K > 1 (8.6.C) are
rejected by the input guard with gdpar_unsupported_feature_error.
Decisions 2.1, 2.2, 2.6, 2.7, 2.8 of CHARTER_SUBFASE_8_6.md.inst/stan/amm_eb_marginal.stan and
inst/stan/amm_eb_conditional.stan. The marginal template is the
cmdstanr::laplace() target of Step (i); the conditional template moves
theta_ref from the parameters block to the data block (as
theta_ref_data) and removes its prior so that Step (iii) HMC samples
the conditional posterior of xi only. Both templates declare the
dimension fields K_slots and p_dim so that Sub-phase 8.6.C can
habilitate $K > 1$ and $p > 1$ without rewriting the data block.
Decisions 2.3, 2.7 of CHARTER_SUBFASE_8_6.md.laplace_control$kappa_threshold, default 1e10; aborts with
gdpar_eb_numerical_error when violated);
(ii) adaptive Levenberg-Marquardt ridge on the marginal covariance
when its eigen-decomposition is singular or non-PD
(laplace_control$ridge_init, default 1e-6);
(iii) multi-start optimize() with laplace_control$multi_start_M
independent random inits (default 5), retaining the init with the
highest joint MAP log-density and reporting the dispersion across
inits in gdpar_eb_fit$diagnostics_numerical$multi_start_dispersion;
(iv) documented fallback that aborts with
gdpar_unsupported_feature_error recommending FB when every init
fails.eb_correction = TRUE): the conditional credible intervals reported
by summary.gdpar_eb_fit are inflated by a scalar factor derived
from the marginal variance of theta_ref; setting
eb_correction = FALSE issues a gdpar_diagnostic_warning advising
of the expected $O(n^{-1})$ under-cover. Decision 2.6 of
CHARTER_SUBFASE_8_6.md.print.gdpar_eb_fit, summary.gdpar_eb_fit,
print.summary.gdpar_eb_fit, coef.gdpar_eb_fit,
predict.gdpar_eb_fit (in-sample only in 8.6.B; newdata rejected
with gdpar_unsupported_feature_error until 8.6.C).test-eb_numerical_robustness.R suite are listed as TODO in
HANDOFF_SUBFASE_8_6_B_PARCIAL.md §5 and are deferred to the next
clean session that closes Sub-phase 8.6.B.gdpar_compare_meta_learners(bridge, methods, newdata, data, seed):
new exported function that benchmarks a gdpar_causal_bridge against a
user-supplied set of external meta-learner adapters on a common
evaluation grid. Reports per-method cate_mean (and native CIs when the
adapter exposes one) plus three concordance matrices (RMSE, Pearson,
mean absolute discrepancy) over every ordered pair of methods.
Descriptive by construction: no tests of hypothesis, no claims of
algorithmic equivalence across methods of different inferential origin.
Decisions 2.1, 2.3, 2.5, 2.6 of CHARTER_SUBFASE_8_5_B.md.gdpar_meta_learner_adapter(name, fit_predict_fun, predict_fun, ...):
constructor of the pluggable adapter contract. The mandatory
fit_predict_fun carries the full fit-and-predict cycle on
(X, Y, T, X_newdata, level, seed_run); the optional predict_fun
re-evaluates the meta-learner on a fresh grid by reusing the fitted
state returned by fit_predict_fun, which lets
predict.gdpar_meta_learner_comparison avoid a refit when supported.
is_gdpar_meta_learner_adapter() and print.gdpar_meta_learner_adapter()
round out the interface.gdpar_adapter_grf() wraps grf::causal_forest R-side with native CIs
from predict(..., estimate.variance = TRUE) and the normal
approximation; gdpar_adapter_econml() wraps
econml.dml.CausalForestDML Python-side via reticulate, with native
CIs from effect_interval(). Both expose fit_predict_fun and
predict_fun; the Python reference inside the EconML adapter's state
becomes invalid after saveRDS round-trip and the predict_fun aborts
cleanly with gdpar_unsupported_feature_error in that case. Decision
2.1 of CHARTER_SUBFASE_8_5_B.md.print.gdpar_meta_learner_comparison for a concise
summary including the three concordance matrices;
summary.gdpar_meta_learner_comparison for a structured object with
long-format pairwise metrics, per-method ATE (with CI bounds derived
from per-observation native CIs when available), and timing;
predict.gdpar_meta_learner_comparison to re-evaluate the CATE on a
new grid for every method, reusing cached state via the adapter's
predict_fun and emitting a gdpar_diagnostic_warning when any
adapter would require a full refit.K > 1
or p > 1 are rejected with gdpar_unsupported_feature_error;
multi-output external adapters are queued for Block 9. Decision 3.8
of CHARTER_SUBFASE_8_5_B.md.Suggests: grf (R-side reference adapter) and reticulate
(Python bridge for EconML). econml itself is a Python module the
user installs in their active Python environment (e.g.
reticulate::py_install("econml")); the package does not install
Python dependencies automatically. Decision 4.6 of
CHARTER_SUBFASE_8_5_B.md.vignette("v08c_meta_learner_comparison") canonizes the comparator
contract, the concordance criterion, and the limits of cross-method
comparison; vignette("vop06_meta_learner_comparison") is the
step-by-step operational recipe (grf CRAN-valid; EconML with chunks
gated by Python availability). Decision 2.2 of
CHARTER_SUBFASE_8_5_B.md.gdpar_causal_bridge(fit_treat, fit_ctrl, newdata, type, level): new
exported function that estimates the conditional average treatment effect
(CATE) from a pair of independent gdpar_fit objects fitted to the two
arms of a treatment / control design. Implements the T-learner
meta-learner of Kuenzel et al. (2019) on the AMM-side: each arm is fitted
independently and the per-observation CATE is the difference of the two
predictive distributions evaluated on a common evaluation set. Decision
2.1 of CHARTER_SUBFASE_8_5_A.md.K > 1), K, p, AMM level, modulating basis type, anchor value,
and the covariate column structure of every AMM component. Mismatches
abort with gdpar_unsupported_feature_error. Hierarchical fits
(stan_data$use_groups == 1L) are out of scope and likewise abort; the
abort uses the same condition class as gdpar_bvm_check so user code
that handles unsupported-feature errors covers both helpers.print.gdpar_causal_bridge for a concise object summary;
summary.gdpar_causal_bridge for a per-observation CATE table plus the
marginal ATE with credible bounds; predict.gdpar_causal_bridge for
re-evaluating the CATE on a fresh newdata without re-running the
compatibility checks.vignette("v08b_cate_ite_bridge_implementation")
canonizing the T-learner AMM-side: definition, identification assumptions
inherited from v02 plus residual no-confounding, estimator,
per-observation credible bounds by posterior quantiles, identifiability
per arm, a minimal reproducible example, limitations of the T-learner,
and open questions (O*-CATE) for Block 9. Continues the CATE/ITE
positioning of vignette("v08_cate_ite_positioning").Suggests: the comparator against external meta-learners
(grf, causalForest, EconML via reticulate) is queued for Sub-phase
8.5.B and adds no dependencies here.coef.gdpar_fit for K-individual fits (K > 1, distributional
regression on amm_distrib_K.stan). Decision E4.A of sub-phase
8.3.10: returns a named list of length K whose entries are
gdpar_coef objects (each with p = 1L), one per slot. The modulating
block W_raw is globally shared across slots in the K-individual
template (canonical decision "Scope of W: global", handoff 28); the
resulting W component is replicated identically across every slot's
gdpar_coef when any slot declared a non-NULL W. K-individual
fits with grouping (J_groups > 1) raise
gdpar_unsupported_feature_error and remain queued for Session 8.4.predict.gdpar_fit with newdata for K-individual fits. Mirrors
predict_from_newdata_multi per slot via the new internal helper
predict_from_newdata_K. Each slot's eta_k is rebuilt from
theta_ref_k[k], a_coef_k[k], c_b_k[k], and the globally-shared
W_raw + sigma_W (with the basis-difference vector evaluated at
the slot-specific theta_ref_k[k] and theta_anchor_K[k]). The
polynomial W branch is supported; B-spline W on new data remains
queued for Session 8.4 (the in-sample path through
apply_W_basis_diff() in Stan continues to support both). Grouping
(J_groups > 1) on new data is also queued for Session 8.4.vignette("vop04_amm_intermediate") — B-spline W bases and
heterogeneous families per slot.vignette("vop05_distributional_K_dharma") — distributional
regression K > 1 API, residual diagnostics G1 / G2 / G3, and
optional DHARMa integration via gdpar_dharma_object().GDPAR_GOLDEN_8_3_9_COMPARE_TIER1=1, ~9 min): 12 configs
across sub-phases 8.3.4 / 8.3.5a / 8.3.6 / 8.3.7 / 8.3.8 / 8.3.9.GDPAR_GOLDEN_8_3_9_COMPARE_TIER2=1, ~40 min): the
Tweedie K=3 config (~38 min by itself) and bspline_p2.inst/extdata/scripts/bootstrap_8_3_9_goldens.R
is now source()-able from the test file; its top-level driver
runs only when invoked as a top-level Rscript with
GDPAR_BOOTSTRAP_8_3_9=1.devtools::document()): coef.gdpar_fit.Rd
updated to describe the named-list return for K > 1; no NAMESPACE
changes.qnorm().pp_check.gdpar_fit() off the bayesplot::pp_check generic
(five PPC types: dens_overlay, hist, ecdf_overlay, stat,
intervals).residuals.gdpar_fit() exported with type ∈ {"quantile", "response", "pearson", "deviance"}; gdpar_posterior_predict()
exported for raw posterior-predictive draws; gdpar_dharma_object()
exported for optional DHARMa integration (DHARMa is a Suggests
dependency, decision E1.A; the package's minimalist Imports is
preserved).K=2, Student-t K=3, Tweedie K=3,
ZIP / ZINB / Hurdle-Poisson / Hurdle-NB, heterogeneous Gauss+Beta /
Gauss+Gamma / NB+Beta (K=2), Gaussian K=1 bspline p=1 / bspline
p=2 / polynomial p=2, and a heterogeneous Gauss+Beta K=2 bspline
config. Manifest CSV with 16 columns of reproducible metadata
(SHA256 fit_code_hash, seed, sub-phase, cmdstan / R / DHARMa
versions, etc.) at tests/testthat/data/golden_manifest.csv
(decision E3.A).W_basis(type = "bspline", knots, degree, boundary_knots) added.
Stan-side Cox-de Boor recursion via bspline_basis_eval() and the
helper apply_W_basis_diff(), used uniformly across the three Stan
templates (amm_main.stan, amm_distrib_multi.stan,
amm_distrib_K.stan).W paths are
bit-exact-preserved by constant propagation; the data field
W_is_polynomial is replaced by W_type_id + W_n_knots_full +
W_knots_full + W_degree.family argument for gdpar() with K = 2
(decision D5). Slot 1 is the canonical primary location family
(gaussian, poisson, neg_binomial_2, beta, gamma,
lognormal_loc_scale); slot 2 may declare any compatible auxiliary
family with coherent support.apply_inv_link_by_id() Stan helper added; refactors the K=2
branches (stan_id 1 / 3 / 5 / 6 / 7) to dispatch the per-slot
inverse link uniformly. Goldens K=2 Beta + Gamma were re-bootstrapped
to absorb the refactor (decision D7=(a)).K=2 (stan_id 10, mu log + pi logit),
ZINB K=3 (stan_id 11, mu log + phi log + pi logit),
Hurdle-Poisson K=2 (stan_id 12), and
Hurdle-NB K=3 (stan_id 13). Hurdle truncation enforced via
log1m_exp(neg_binomial_2_log_lpmf(0 | mu, phi)) for numerical
stability.K=3 (stan_id 8, mu identity + sigma log + nu log).K=3 (stan_id 9, mu log + phi log + p identity bounded (1.01, 1.99)). Custom tweedie_lpdf in the functions{} block via
a hybrid Dunn-Smyth series + Wood (2017) saddlepoint switchover at
tau = 0.4; tweedie_rng via compound Poisson-gamma.K=2 (stan_id 5, mu logit + phi log),
Gamma K=2 (stan_id 6, mu log + shape log),
log-normal K=2 (stan_id 7, custom-family wired via descriptor).gdpar_bf(y ~ a(x1), sigma ~ a(x2)) constructor + gdpar_formula_set S3 class. Three equivalent input forms canonicalise to the same internal representation.individual_params argument removed from gdpar_family() (decision 4C superseded by canonisation D). gdpar_family() now emits all eligible param_specs; the input to gdpar() is the single source of truth for K.amm_spec_designs(),
amm_spec_multi_designs()), and no longer by a sum-to-zero
reparametrization on the basis coefficients inside Stan.colMeans(Z_a) = 0 enforced empirically, the
expectation E_mu[a(X)] = colMeans(Z_a) * a_coef vanishes for any
a_coef, so (C2) is satisfied for the full J_a-dimensional space
of basis coefficients. The previous implementation enforced the
same centering twice (in R and again in Stan), and the second
enforcement, sum-to-zero on a_coef, was a strictly stronger
restriction not derived from (C2): it confined a_coef to the
J_a - 1-dimensional subspace {a : sum(a) = 0} for every fit.a_coef = (1.0, -0.8) on two
continuous covariates) are now recovered correctly. Previously they
were projected onto the sum-zero subspace at fit time, inflating
residual variance and degrading ELPD. The S7 (heterog_p3) scenario
of the Block 6 adversarial bench, where gdpar previously drifted
88 ELPD points outside 2*se versus the other four competitors, now
cleanly enters the indistinguishable cluster (gdpar = -1378.45 vs
best mgcv = -1375.50, delta 2.95 within 2*se = 70).inst/stan/amm_main.stan,
inst/stan/amm_main_multi.stan, R/stan_codegen.R
(generate_a_blocks_multi NCP/CP/mixed variants), R/preflight.R,
R/preflight_multi.R. Documentation revised in
vignette("v01_amm_identifiability") Section 10.1 and
vignette("v03_special_cases") Sections 6 and 8. Roxygen for
amm_spec() and gdpar_prior() reworded. The Block 6 multi golden
fits (smoke_p2_auto.rds, smoke_p2_cp.rds, smoke_p2_ncp.rds)
were regenerated against the new codegen.expect_false for the absent
sum-to-zero emission on the scalar template).gdpar() gains argument id_check_rigor = c("full", "fast") with
default "full" (backward-compatible). The flag is propagated to
.gdpar_multi() and to gdpar_check_identifiability(..., rigor = id_check_rigor). The "fast" value downgrades the C4-bis
cross-coordinate check to a single consolidated warning at the end
of the per-coordinate loop instead of aborting, enabling fits with
intentional overlap between the additive and modulating channels
combined with informative priors. Documented in ?gdpar and in
vignette("v01_amm_identifiability", package = "gdpar") Section
6.6.1.7.gdpar_loo() exported (flagged @keywords experimental): PSIS-LOO
approximate cross-validation via loo::loo() on the per-observation
log-likelihood emitted by the Stan template. Two aggregations:
"subject" (default; sums the coord-wise log-likelihoods to
per-subject scale, matching brms::brm() multivariate fits with
set_rescor(FALSE)) and "cell" (treats each (i, k) pair as an
independent observation, useful for per-coordinate Pareto-k
diagnostics).inst/benchmarks/:
inst/benchmarks/synthetic_hard_multi.R: 9 scenario generators
(S0..S8) covering sanity, adversarial, and stress tiers with
n_train = 800, n_test = 200 per scenario. Includes a
Bernoulli p = 2 adversarial scenario (S5) and a Gaussian
p = 2 residual-correlation scenario (S8) where brms is run a
second time with set_rescor(TRUE) under the label brms_rescor
to verify recovery of the true residual correlation.inst/benchmarks/scripts/bench_competitor_<method>.R wrappers
for gdpar, brms, mgcv, INLA, and rstanarm. PSIS-LOO is computed
on principled posterior draws for all five methods: gdpar /
brms / rstanarm use their native Stan posteriors; mgcv uses the
Bayesian interpretation of REML (Wood 2017, sec 6.10) via
mvnfast::rmvn(mu = coef(gam_obj), sigma = vcov(gam_obj)) to
sample from the gaussian Laplace approximation; INLA uses
inla.posterior.sample() with a CPO-sum cross-check that flags
inla_cpo_loo_divergent when the two values disagree by more
than 2 * se_elpd.inst/benchmarks/run_synthetic_hard_multi.R runner, env-gated
by GDPAR_BENCH_MULTI=1, writes the long-tidy CSV
synthetic_hard_multi_results.csv with columns including
max_pareto_k, n_pareto_k_above_07, convergence_flag, and
recovered_rescor (NA except for brms_rescor rows in S8).inst/benchmarks/scripts/report_synthetic_hard_multi.R
visualizer producing synthetic_hard_multi_plot.png,
synthetic_hard_multi_wall.png,
synthetic_hard_multi_pareto.png, and a markdown verdict table
synthetic_hard_multi_summary.md.vignettes/vop02_arbitrary_p.Rmd extended with Section 10
"PSIS-LOO via gdpar_loo()" documenting both aggregations,
Pareto-k caveats, and the experimental-status flag.Suggests: adds brms, INLA, rstanarm, scoringRules,
gratia, mvnfast for the Block 6 bench harness. All entries are
optional and gated by requireNamespace() at call sites.Additional_repositories: adds
https://inla.r-inla-download.org/R/stable so that INLA can be
installed from CRAN-compatible tooling.