Standardized-cell Z-ratio build and Cauchy-slab constants

NAMESPACE

@@ -90,6 +90,7 @@ importFrom(S7,new_property)
 importFrom(S7,prop)
 importFrom(methods,hasArg)
 importFrom(stats,approxfun)
+importFrom(stats,dcauchy)
 importFrom(stats,dnorm)
 importFrom(stats,dpois)
 importFrom(stats,integrate)

NEWS.md

@@ -30,8 +30,10 @@
 * Corrected inclusion-probability updates for `beta_bernoulli_prior()` on continuous (GGM) models: under the determinant-tilted precision prior, the conjugate Beta update omits a normalizing-constant factor and biases the sampled inclusion probability toward sparsity (at 5 variables with a uniform hyperprior its prior mean lands near 0.37 instead of 0.5). The update now draws from the corrected conditional using a table built from the prior distribution at the first fit of a model configuration and cached on disk (`tools::R_user_dir("bgms", "cache")`; a one-time cost of the order of minutes, announced when `verbose = TRUE`). The sampled inclusion probability is returned per chain in `fit$inclusion_parameter_samples`.
 * Corrected stochastic block model updates for `sbm_prior()` on continuous (GGM) models, using the same cached table: the block-probability draws, the block-allocation weights, and the new-cluster weight all carry the normalizing-constant correction. Without it the sampled partition collapses toward one block (at 20 variables the prior mean number of blocks lands near 1.0 instead of 1.87); with it a prior-only chain reproduces the model's partition prior (total variation 0.01-0.02 at 5-20 variables). `sample_ggm_prior()` accepts `edge_prior` objects (`bernoulli_prior()`, `beta_bernoulli_prior()`, `sbm_prior()`) and returns sampled allocations under the block model.
 * The normalizing-constant correction extends to mixed models with `beta_bernoulli_prior()` or `sbm_prior()`: the determinant tilt acts on the continuous precision block, so the correction table is built for the continuous variables and the block-structure corrections read continuous-continuous edges only. With fewer than two continuous variables no edge is tilted and the plain conjugate updates apply unchanged; with exactly two, the single tilted pair supports the inclusion-probability correction but not the block-model slope curve, so `sbm_prior()` warns and keeps the plain conjugate updates there. Prior-only mixed chains reproduce the Beta hyperprior on the inclusion probability and the partition prior on the number of blocks; without the correction the inclusion probability biases toward sparsity and the partition toward fewer blocks.
-* Hierarchical graph-prior specification for continuous (GGM) models: `bgm(precision_graph_prior = "hierarchical")` composes the edge prior and the precision prior as `p(Gamma) p(K | Gamma)` with `p(K | Gamma)` normalized per graph, so the graph marginal is exactly the edge prior (under the joint specification it is reweighted by the per-graph normalizer). Each between-edge move evaluates the normalizer ratio by a deterministic local Z-ratio approximation, calibrated online against a block-Gibbs oracle in an appended warm-up window (`calibration_window`) and frozen with a hull clamp before sampling. Requires `edge_selection = TRUE`, a `normal_prior()` interaction prior, and a shape-1 `gamma_prior()` (or `exponential_prior()`) precision scale prior. The same specification is available in `sample_ggm_prior(spec = "hierarchical")`. With-data simulation-based calibration at 50 variables, n = 25, Beta-Bernoulli(2, 4): posterior inclusion-probability bias +0.002 with uniform ranks (the additive approximation alone reads -0.03). Prior-only chains with a sampled inclusion probability at large dimension remain out of envelope (the alarm suite flags them). On mixed models the specification normalizes the continuous block `p(K_yy | Gamma_yy)`: the Z-ratio enters the continuous-continuous edge moves only (counts on the continuous subgraph), and at least two continuous variables are required.
+* Hierarchical graph-prior specification for continuous (GGM) models: `bgm(precision_graph_prior = "hierarchical")` composes the edge prior and the precision prior as `p(Gamma) p(K | Gamma)` with `p(K | Gamma)` normalized per graph, so the graph marginal is exactly the edge prior (under the joint specification it is reweighted by the per-graph normalizer). Each between-edge move evaluates the normalizer ratio by a deterministic local Z-ratio approximation, calibrated online against a block-Gibbs oracle in an appended warm-up window (`calibration_window`) and frozen with a hull clamp before sampling. Requires `edge_selection = TRUE`, a `normal_prior()` or `cauchy_prior()` interaction prior, and a shape-1 `gamma_prior()` (or `exponential_prior()`) precision scale prior. The same specification is available in `sample_ggm_prior(spec = "hierarchical")`. With-data simulation-based calibration at 50 variables, n = 25, Beta-Bernoulli(2, 4): posterior inclusion-probability bias +0.002 with uniform ranks (the additive approximation alone reads -0.03). Prior-only chains with a sampled inclusion probability at large dimension remain out of envelope (the alarm suite flags them). On mixed models the specification normalizes the continuous block `p(K_yy | Gamma_yy)`: the Z-ratio enters the continuous-continuous edge moves only (counts on the continuous subgraph), and at least two continuous variables are required.
 * Z-ratio alarm suite: `summarize_zratio_diagnostics()` audits the frozen Z-ratio kernel on graphs each chain visited (targeted, random, and additive-zone channels against a measurement-only block-Gibbs oracle), checks the calibration stream for end-of-warmup drift, and reports regime context. The verdict compares the maximum audit error to a per-regime threshold and prints like other sampler warnings; the summary is attached as `fit$zratio_diag` (and `$zratio_diagnostics` on `sample_ggm_prior()` output).
+* The Z-ratio constants for the hierarchical specification are built in the standardized cell (slab scale 1, diagonal rate `eta = pairwise_scale * scale_rate`). The between-graph normalizer ratio is invariant to the slab scale at fixed `eta`, and the estimator's quadrature grids are sized for the unit-scale frame, so this makes the Z-ratio numerically independent of the user's scale choice. Previously the constants were built at the raw slab scale, where a `normal_prior(scale = 2.5)` fit read the saddle 0.10 low in `log Z(Gamma-)/Z(Gamma+)` and the `sample_ggm_prior()` hierarchical graph marginal at 6 variables sat at 0.292 for a 0.30 edge prior; unit-scale fits (`scale = 0.5`) are unchanged.
+* A `cauchy_prior()` interaction prior under the hierarchical specification now gets slab-family-specific Z-ratio constants: the per-graph normalizer is the marginal-Cauchy constant (the IG(1/2, 1/2) scale mixture of the Normal one), integrated per channel by an exact resolvent-form mixture rule, and the block-Gibbs calibration oracle runs omega-augmented with matching leg-dressed moments. The channel constants are verified against direct Monte Carlo of the same moments. Previously the Cauchy slab reused the Normal-slab constants at the base scale, which broke the hierarchical graph law: at 6 variables the graph marginal read 0.247 for a 0.30 edge prior; it now sits within Monte Carlo error of the target on both the adaptive-metropolis and gibbs update methods.
 * The NUTS diagnostics summary now prints the `warmup_incomplete` flag (energy not stationary) it already computed.
 * `extract_prior_inclusion_probabilities()`: prior edge-inclusion probabilities in the same matrix layout as `extract_posterior_inclusion_probabilities()`, for prior/posterior inclusion-odds computations. Under the joint spike-and-slab prior on a continuous block the graph marginal is reweighted by the per-graph normalizer (positive-definite-cone mass shaped by the determinant tilt), so continuous-continuous edges do not keep the edge-prior marginal at any `delta` — with a uniform Beta-Bernoulli hyperprior at 3 variables the prior edge probability is about 0.37 rather than 0.5, and a fixed `bernoulli_prior(0.5)` at `delta = 0` lands near 0.27. The values are read from the cached correction table (`bernoulli_prior()`, `beta_bernoulli_prior()`) or estimated by a prior-only chain with the fit's own correction settings (`sbm_prior()`, cached on the fit). Mixed models report per-edge-class values (discrete-discrete, continuous-continuous, cross); ordinal models use the analytic edge-prior marginals, including the exchangeable partition mixture for `sbm_prior()`.
 

R/RcppExports.R

@@ -161,8 +161,8 @@ zratio_scan_graph <- function(G, addc, tg, ihat, ghat, wt, psi0) {
     .Call(`_bgms_zratio_scan_graph`, G, addc, tg, ihat, ghat, wt, psi0)
 }
 
-zratio_audit_edges <- function(G, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, n_sweep, burn, seed) {
-    .Call(`_bgms_zratio_audit_edges`, G, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, n_sweep, burn, seed)
+zratio_audit_edges <- function(G, edges, addc, tg, ihat, ghat, wt, psi0, delta, eta, n_sweep, burn, seed, slab_cauchy = FALSE) {
+    .Call(`_bgms_zratio_audit_edges`, G, edges, addc, tg, ihat, ghat, wt, psi0, delta, eta, n_sweep, burn, seed, slab_cauchy)
 }
 
 zratio_test_eval <- function(G, edges, addc, tg, ihat, ghat, wt, psi0) {
@@ -173,8 +173,8 @@ zratio_test_saddle <- function(s1, s2, addc, tg, ihat, ghat, wt, psi0) {
     .Call(`_bgms_zratio_test_saddle`, s1, s2, addc, tg, ihat, ghat, wt, psi0)
 }
 
-zratio_test_calibrated_eval <- function(graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, seed, n_sweep, burn, freeze_after) {
-    .Call(`_bgms_zratio_test_calibrated_eval`, graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, seed, n_sweep, burn, freeze_after)
+zratio_test_calibrated_eval <- function(graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, eta, seed, n_sweep, burn, freeze_after, slab_cauchy = FALSE) {
+    .Call(`_bgms_zratio_test_calibrated_eval`, graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, eta, seed, n_sweep, burn, freeze_after, slab_cauchy)
 }
 
 zratio_test_precompute <- function(G, edges, addc, tg, ihat, ghat, wt, psi0, ncn_max, bre_max) {

R/bgms-package.R

@@ -64,7 +64,7 @@
 #' @docType package
 #' @keywords internal
 #' @useDynLib bgms, .registration=TRUE
-#' @importFrom stats approxfun dnorm integrate
+#' @importFrom stats approxfun dcauchy dnorm integrate
 #' @importFrom stats rbeta rexp rgamma rnorm rpois runif
 #' @references
 #' \insertAllCited{}

R/build_output_bgm.R

@@ -294,10 +294,9 @@ build_output_bgm = function(spec, raw) {
 
   # --- Z-ratio alarm suite (hierarchical graph-prior spec) ----------------------
   if(!is.null(zratio_chains)) {
-    zc = zratio_constants(
-      delta = p$delta,
-      sigma = 2 * p$pairwise_scale,
-      beta = p$scale_rate / 2
+    zc = zratio_cell_constants(
+      p$delta, p$pairwise_scale, p$scale_rate, p$scale_eta,
+      slab = p$interaction_prior_type
     )
     results$zratio_diag = summarize_zratio_diagnostics(
       zratio_chains, zc,

R/build_output_mixed_mrf.R

@@ -289,10 +289,9 @@ build_output_mixed_mrf = function(spec, raw) {
   # The indicator vector is [Gxx upper, Gyy upper, Gxy row-major], both
   # triangles without diagonals; the audit reads the Gyy segment.
   if(!is.null(zratio_chains)) {
-    zc = zratio_constants(
-      delta = pr$delta,
-      sigma = 2 * pr$pairwise_scale,
-      beta = pr$scale_rate / 2
+    zc = zratio_cell_constants(
+      pr$delta, pr$pairwise_scale, pr$scale_rate, pr$scale_eta,
+      slab = pr$interaction_prior_type
     )
     results$zratio_diag = summarize_zratio_diagnostics(
       zratio_chains, zc,

R/run_sampler.R

@@ -78,15 +78,14 @@ run_sampler_ggm = function(spec) {
   correction = NULL
   zratio = NULL
   if(identical(p$precision_graph_prior, "hierarchical")) {
-    zc = zratio_constants(
-      delta = p$delta,
-      sigma = 2 * p$pairwise_scale,
-      beta = p$scale_rate / 2
+    zc = zratio_cell_constants(
+      p$delta, p$pairwise_scale, p$scale_rate, p$scale_eta,
+      slab = p$interaction_prior_type
     )
     zratio = list(
       addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat,
       wt = zc$wt, psi0 = zc$psi0,
-      delta = zc$delta, sigma = zc$sigma, beta = zc$beta,
+      delta = zc$delta, eta = zc$eta, slab = zc$slab,
       calibration_window = resolve_zratio_calibration_window(
         p$calibration_window, d$num_variables, s$warmup
       )
@@ -222,15 +221,14 @@ run_sampler_mixed_mrf = function(spec) {
   correction = NULL
   zratio = NULL
   if(identical(p$precision_graph_prior, "hierarchical")) {
-    zc = zratio_constants(
-      delta = p$delta,
-      sigma = 2 * p$pairwise_scale,
-      beta = p$scale_rate / 2
+    zc = zratio_cell_constants(
+      p$delta, p$pairwise_scale, p$scale_rate, p$scale_eta,
+      slab = p$interaction_prior_type
     )
     zratio = list(
       addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat,
       wt = zc$wt, psi0 = zc$psi0,
-      delta = zc$delta, sigma = zc$sigma, beta = zc$beta,
+      delta = zc$delta, eta = zc$eta, slab = zc$slab,
       calibration_window = resolve_zratio_calibration_window(
         p$calibration_window, d$num_continuous, s$warmup
       )

R/sample_ggm_prior.R

@@ -336,15 +336,14 @@ sample_ggm_prior = function(
         format(sp$scale_shape)
       ))
     }
-    zc = zratio_constants(
-      delta = delta,
-      sigma = 2 * ip$pairwise_scale,
-      beta = sp$scale_rate / 2
+    zc = zratio_cell_constants(
+      delta, ip$pairwise_scale, sp$scale_rate, sp$scale_eta,
+      slab = ip$interaction_prior_type
     )
     zratio = list(
       addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat,
       wt = zc$wt, psi0 = zc$psi0,
-      delta = zc$delta, sigma = zc$sigma, beta = zc$beta,
+      delta = zc$delta, eta = zc$eta, slab = zc$slab,
       calibration_window = resolve_zratio_calibration_window(
         calibration_window, p, n_warmup
       )

R/zratio_diagnostics.R

@@ -13,15 +13,15 @@
 # ------------------------------------------------------------------------------
 # zratio_eta
 # ------------------------------------------------------------------------------
-# Standardized prior-scale regime index: eta = sigma * beta in the bare frame
-# (sigma = 1, beta = eta after standardization). Drives the verdict threshold.
+# Standardized prior-scale regime index: the diagonal rate eta in the
+# standardized frame (unit slab scale). Drives the verdict threshold.
 #
-# @param zratio_spec  Z-ratio spec list with bare-scale sigma and beta.
+# @param zratio_spec  Z-ratio spec list carrying the standardized eta.
 #
 # Returns: Numeric scalar eta.
 # ------------------------------------------------------------------------------
 zratio_eta = function(zratio_spec) {
-  zratio_spec$sigma * zratio_spec$beta
+  zratio_spec$eta
 }
 
 # ------------------------------------------------------------------------------
@@ -277,8 +277,9 @@ zratio_audit_chain = function(chain, zratio_spec, num_nodes, n_graphs, top_k,
       graphs[[g]], edges,
       addc, zratio_spec$tg, zratio_spec$ihat, zratio_spec$ghat,
       zratio_spec$wt, zratio_spec$psi0,
-      zratio_spec$delta, zratio_spec$sigma, zratio_spec$beta,
-      as.integer(audit_sweep), 30L, as.integer(seed + g)
+      zratio_spec$delta, zratio_spec$eta,
+      as.integer(audit_sweep), 30L, as.integer(seed + g),
+      identical(zratio_spec$slab, "cauchy")
     )
     err[sel] = audit$err
   }
@@ -318,8 +319,9 @@ zratio_audit_chain = function(chain, zratio_spec, num_nodes, n_graphs, top_k,
 #' visited-density band). The verdict flags a chain when the maximum of its
 #' targeted (or, without a fit, random) and additive-zone audit errors
 #' exceeds the regime threshold: 0.01 at \eqn{\eta = 1}, 0.02 at
-#' \eqn{\eta = 2}, 0.04 at \eqn{\eta \ge 3}, with
-#' \eqn{\eta = \sigma \beta} in the bare prior scale. A quiet verdict bounds
+#' \eqn{\eta = 2}, 0.04 at \eqn{\eta \ge 3}, where \eqn{\eta} is the
+#' diagonal rate in the standardized (unit slab scale) frame. A quiet
+#' verdict bounds
 #' the pointwise approximation error on visited graphs; coherent sub-margin
 #' bias that accumulates through inclusion-probability feedback (prior-only
 #' chains at large \eqn{p}) is outside its reach.
@@ -331,8 +333,8 @@ zratio_audit_chain = function(chain, zratio_spec, num_nodes, n_graphs, top_k,
 #'   active.
 #' @param zratio_spec The Z-ratio specification list used for the run:
 #'   quadrature tables \code{tg}, \code{ihat}, \code{ghat}, \code{wt},
-#'   isolated-edge ratio \code{psi0}, and the bare-scale constants
-#'   \code{delta}, \code{sigma}, \code{beta}.
+#'   isolated-edge ratio \code{psi0}, and the standardized-cell constants
+#'   \code{delta} and \code{eta}.
 #' @param num_nodes Integer: number of nodes \eqn{p}.
 #' @param n_graphs Integer: visited graphs scanned per chain (default 12).
 #' @param top_k Integer: picks per targeted rule and additive-zone top-m
@@ -361,7 +363,7 @@ zratio_audit_chain = function(chain, zratio_spec, num_nodes, n_graphs, top_k,
 #'       flags (\code{drift_density}, \code{drift_theta}), and regime
 #'       context (\code{a1a_out_frac}, \code{m32_frac},
 #'       \code{dens_out_frac}).}
-#'     \item{\code{eta}}{Regime index \eqn{\sigma \beta}.}
+#'     \item{\code{eta}}{Regime index: the standardized diagonal rate.}
 #'     \item{\code{tau}}{Verdict threshold at this \code{eta}.}
 #'     \item{\code{verdict_flagged}}{Logical: any chain flagged.}
 #'     \item{\code{calibration_incomplete}}{Logical: any chain's