The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
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).
sigma_a_k
compaction (Session B9.30, 2026-06-03)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
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 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
+ B base implementation + C multivariate partial paths + D full Path C +
E comparator and operational vignette). For the canonical opening +
canonization summary of Sub-fase 8.6 see
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).
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);laplace_control$ridge_init, default
1e-6);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;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 bases (2026-05-22)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 likelihoods (2026-05-21)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 likelihoods (2026-05-21)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).K > 1
(2026-05-20)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.These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.