From 03b5fe56412fd99f79a2fa4d88a958b808338248 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Mon, 6 Jul 2026 23:29:11 +0200 Subject: [PATCH 1/4] Build the Z-ratio constants in the standardized cell The between-graph normalizer ratio is invariant under the diagonal congruence Theta = A K A, so it depends on (delta, eta) with eta = pairwise_scale * scale_rate, and the estimator's fixed quadrature grids are sized for the sigma = 1 frame. All call sites now resolve the constants at (delta, sigma = 1, beta = eta) via zratio_cell_constants(), and the calibration oracle receives the same cell. At the previous bare-scale build a scale-2.5 fit read the saddle 0.10 low in log Z(Gamma-)/Z(Gamma+) and the q = 6 hierarchical graph marginal sat at 0.292 for a 0.30 edge prior; unit-scale cells are bit-identical before and after. --- NEWS.md | 1 + R/build_output_bgm.R | 6 ++-- R/build_output_mixed_mrf.R | 6 ++-- R/run_sampler.R | 12 +++----- R/sample_ggm_prior.R | 6 ++-- R/zratio_tables.R | 26 +++++++++++++--- src/models/ggm/zratio_engine.h | 23 ++++++++------ tests/testthat/test-hier-zratio-identity.R | 35 ++++++++++++++++++++++ 8 files changed, 82 insertions(+), 33 deletions(-) diff --git a/NEWS.md b/NEWS.md index 3087e7db..60d0b12e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -32,6 +32,7 @@ * 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. * 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. * 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()`. diff --git a/R/build_output_bgm.R b/R/build_output_bgm.R index b632f4d7..f2b17eb2 100644 --- a/R/build_output_bgm.R +++ b/R/build_output_bgm.R @@ -294,10 +294,8 @@ 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 ) results$zratio_diag = summarize_zratio_diagnostics( zratio_chains, zc, diff --git a/R/build_output_mixed_mrf.R b/R/build_output_mixed_mrf.R index 0b42c0b3..21030093 100644 --- a/R/build_output_mixed_mrf.R +++ b/R/build_output_mixed_mrf.R @@ -289,10 +289,8 @@ 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 ) results$zratio_diag = summarize_zratio_diagnostics( zratio_chains, zc, diff --git a/R/run_sampler.R b/R/run_sampler.R index f48cd716..c42388cc 100644 --- a/R/run_sampler.R +++ b/R/run_sampler.R @@ -78,10 +78,8 @@ 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 ) zratio = list( addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, @@ -222,10 +220,8 @@ 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 ) zratio = list( addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R index 7311abd6..64c7e942 100644 --- a/R/sample_ggm_prior.R +++ b/R/sample_ggm_prior.R @@ -336,10 +336,8 @@ 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 ) zratio = list( addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, diff --git a/R/zratio_tables.R b/R/zratio_tables.R index 69b6f325..06461975 100644 --- a/R/zratio_tables.R +++ b/R/zratio_tables.R @@ -10,10 +10,16 @@ # - the per-channel moment constants addc[0..5] for the additive counts # (common-neighbour node, CN-CN edge, bridge edge), # - psi0, the isolated-edge ratio I_spike(0)/G(0). -# Conventions (bare scale): K_ii ~ Exp(beta), slab K_ij ~ N(0, sigma^2), -# tilt |K|^delta. In bgms parameter units: sigma = 2 * pairwise_scale and -# beta = scale_rate / 2 (priors act on K/2). Reference implementation: -# SV/Z don-validation (sd_marginal_helpers.R, ks_validation_grid.R). +# Conventions: K_ii ~ Exp(beta), slab K_ij ~ N(0, sigma^2), tilt +# |K|^delta. The between-graph ratio Z(Gamma-)/Z(Gamma+) is invariant under +# the diagonal congruence Theta = A K A (scale standardization of the +# normalizer), so it depends on (delta, eta) alone, where in bgms parameter +# units eta = (2 * pairwise_scale) * (scale_rate / 2) = pairwise_scale * +# scale_rate (priors act on K/2). The fixed quadrature grids below (cmax, +# Cmax, Tmax, Laguerre ranges) are sized for the sigma = 1 frame, so every +# consumer builds in the standardized cell (delta, sigma = 1, beta = eta) +# via zratio_cell_constants(). Reference implementation: SV/Z don-validation +# (sd_marginal_helpers.R, ks_validation_grid.R). # Golub-Welsch Gauss quadrature nodes/weights. kind: "laguerre" (weight # e^{-x} on (0, Inf)), "hermite" (weight e^{-x^2} on (-Inf, Inf)), @@ -213,6 +219,18 @@ zratio_bridge_channel = function(delta, sigma, beta, nlag = 64, nher = 80) { c(n1 / den, n2 / den) } +# Standardized cell for one fit's Z-ratio constants. The between-graph +# ratio depends on (delta, eta) only, so the constants are built at +# sigma = 1, beta = eta, the frame the quadrature grids are sized for; user +# scale choices with the same eta share one constant set. eta is the +# user-specified standardized rate when the scale prior carries one, else +# pairwise_scale * scale_rate (the same number up to rounding). +zratio_cell_constants = function(delta, pairwise_scale, scale_rate, + scale_eta = NA_real_) { + eta = if(is.finite(scale_eta)) scale_eta else pairwise_scale * scale_rate + zratio_constants(delta, sigma = 1, beta = eta) +} + # Session cache for zratio_constants: the constant set is deterministic per # (delta, sigma, beta) cell, and one fit resolves the same cell more than # once (sampler dispatch and diagnostics assembly). diff --git a/src/models/ggm/zratio_engine.h b/src/models/ggm/zratio_engine.h index 8bad3219..92c27d7c 100644 --- a/src/models/ggm/zratio_engine.h +++ b/src/models/ggm/zratio_engine.h @@ -53,10 +53,14 @@ struct ZRatioBlock { * design, order bre, m, cne, maxbd, dens); outside the box the * correction is zeroed so the frozen kernel never extrapolates. * - * Conventions (bare scale): K_ii ~ Exp(beta), slab K_ij ~ N(0, sigma^2), - * tilt |K|^delta; sigma = 2 * pairwise_scale, beta = scale_rate / 2 in bgms - * parameter units. Reference: SV/Z sbc_prior_chain_exact.cpp and the - * z_graph_prior deployed kernel (hier_chain_data.cpp). + * Conventions (standardized cell): K_ii ~ Exp(beta), slab K_ij ~ N(0, + * sigma^2), tilt |K|^delta. The between-graph ratio is invariant under the + * diagonal congruence Theta = A K A, so the constants are built at + * sigma = 1, beta = eta = pairwise_scale * scale_rate in bgms parameter + * units (R/zratio_tables.R, zratio_cell_constants), and the same cell + * serves every user scale choice. Reference: SV/Z + * sbc_prior_chain_exact.cpp and the z_graph_prior deployed kernel + * (hier_chain_data.cpp). */ class ZRatioEngine { public: @@ -113,9 +117,9 @@ class ZRatioEngine { double& oracle_out, ZRatioBlock& bl_out); /** - * Set the bare-scale prior constants and RNG the block-Gibbs oracle - * samples under, without entering calibration mode. rng must outlive - * the engine. + * Set the standardized-cell prior constants and RNG the block-Gibbs + * oracle samples under, without entering calibration mode. rng must + * outlive the engine. */ void set_oracle_params(double delta, double sigma, double beta, SafeRNG* rng, int n_sweep = 300, int burn = 30) { @@ -138,8 +142,9 @@ class ZRatioEngine { * box into the addc layout, after which the engine behaves exactly * like one constructed with a full 23-slot constant block. * - * (delta, sigma, beta) are the bare-scale prior constants the oracle - * samples under; rng must outlive the engine (the model's chain RNG). + * (delta, sigma, beta) are the standardized-cell prior constants the + * oracle samples under; rng must outlive the engine (the model's + * chain RNG). */ void enable_calibration(double delta, double sigma, double beta, SafeRNG* rng, int n_sweep = 300, int burn = 30, diff --git a/tests/testthat/test-hier-zratio-identity.R b/tests/testthat/test-hier-zratio-identity.R index abc9c877..ce6bacb3 100644 --- a/tests/testthat/test-hier-zratio-identity.R +++ b/tests/testthat/test-hier-zratio-identity.R @@ -53,6 +53,41 @@ test_that("hierarchical graph marginal matches Bernoulli(p); joint does not", { } }) +test_that("Z-ratio constants build in the standardized cell", { + skip_on_cran() + # The between-graph ratio depends on (delta, eta) only, so every user + # frame with the same eta = pairwise_scale * scale_rate must resolve to + # one constant set, built at sigma = 1 (the frame the quadrature grids + # are sized for). At the old bare-scale build the scale-2.5 frame drifted + # the saddle by -0.10 in log J and the q = 6 identity by -0.008. + d = 0.5 * log(6) + a = bgms:::zratio_cell_constants(d, pairwise_scale = 0.5, scale_rate = 2) + b = bgms:::zratio_cell_constants(d, pairwise_scale = 0.25, scale_rate = 4) + e = bgms:::zratio_cell_constants( + d, pairwise_scale = 2.5, scale_rate = 0.4, scale_eta = 1 + ) + expect_identical(a, b) + expect_identical(a, e) + expect_identical(a$sigma, 1) + expect_identical(a$beta, 1) +}) + +test_that("hierarchical graph marginal holds at a non-unit slab scale", { + skip_on_cran() + # pairwise scale 2.5 (the sample_ggm_prior default) with eta = 1 is the + # same physical cell as scale 0.5 / rate 2; the graph law must hold there + # identically. + d = sample_ggm_prior( + p = 6L, n_samples = 6000L, n_warmup = 1500L, + interaction_prior = normal_prior(scale = 2.5), + precision_scale_prior = gamma_prior(shape = 1, eta = 1), + spec = "hierarchical", edge_inclusion_prob = 0.3, + update_method = "adaptive-metropolis", delta = 0.5 * log(6), + seed = 11L, verbose = FALSE + ) + expect_lt(abs(mean(d$edge_indicators) - 0.3), 0.02) +}) + test_that("hierarchical BB identity: theta ~ Beta(a, b), PIP = a/(a+b)", { skip_on_cran() d = hier_prior_run( From 385ebee06117d46ebac3d181ff7f992bf2433a50 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Mon, 6 Jul 2026 23:53:08 +0200 Subject: [PATCH 2/4] Slab-family-specific Z-ratio constants for the Cauchy prior Under the hierarchical specification a Cauchy interaction prior reused the Normal-slab Z-ratio constants at the base scale, which broke the graph-law identity (q = 6 marginal 0.247 for a 0.30 edge prior). The per-graph normalizer for the Cauchy slab is the marginal-Cauchy constant Z_C(Gamma) = E_omega[Z_N(Gamma; sigma sqrt(omega))], which keeps the omega updates conjugate and the between-move structure unchanged; only the constants and the calibration oracle are family-specific. zratio_constants() gains a slab argument. The Cauchy builders swap the slab density in the pair integral and mix each channel over the IG(1/2, 1/2) leg weights through the exact resolvent-form rule E[omega^k (A + B omega)^-(k+1/2)] = sqrt(2/pi) sum w_i (2 t_i A + B)^-(k+1/2) on the plain Gauss-Laguerre grid; the clique-2 chain runs omega-augmented with leg-dressed block moments. The C++ block-Gibbs oracle gets the same augmentation (conjugate block-edge weights, fresh per-sweep leg weights) behind a slab flag carried in zratio_spec and the audit interfaces. Verification: node, bridge, and psi0 constants within 0.14% of direct Monte Carlo of the same F-measure moments (standing tests in test-zratio-cauchy.R); the omega-augmented oracle sits at the additive prediction on the exact no-edge baseline block; the q = 6 Cauchy hierarchical graph marginal reads +0.0005 (adaptive-metropolis) and -0.0038 (gibbs) against the 0.30 target; the with-data Cauchy battery (three update methods, calibration window, alarm suite) passes. --- NAMESPACE | 1 + NEWS.md | 3 +- R/RcppExports.R | 8 +- R/bgms-package.R | 2 +- R/build_output_bgm.R | 3 +- R/build_output_mixed_mrf.R | 3 +- R/run_sampler.R | 10 +- R/sample_ggm_prior.R | 5 +- R/zratio_diagnostics.R | 3 +- R/zratio_tables.R | 161 ++++++++++++++++++--- src/RcppExports.cpp | 18 ++- src/models/ggm/zratio_engine.cpp | 71 ++++++++- src/models/ggm/zratio_engine.h | 21 ++- src/sample_ggm.cpp | 5 +- src/sample_mixed.cpp | 5 +- src/zratio_audit_interface.cpp | 9 +- src/zratio_test_interface.cpp | 6 +- tests/testthat/test-hier-zratio-identity.R | 19 ++- tests/testthat/test-zratio-cauchy.R | 116 +++++++++++++++ 19 files changed, 400 insertions(+), 69 deletions(-) create mode 100644 tests/testthat/test-zratio-cauchy.R diff --git a/NAMESPACE b/NAMESPACE index c7e82c2f..dc89b5bb 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index 60d0b12e..24fc598d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -30,9 +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()`. diff --git a/R/RcppExports.R b/R/RcppExports.R index 0904e393..30477814 100644 --- a/R/RcppExports.R +++ b/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, sigma, beta, n_sweep, burn, seed, slab_cauchy = FALSE) { + .Call(`_bgms_zratio_audit_edges`, G, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, 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, sigma, beta, seed, n_sweep, burn, freeze_after, slab_cauchy = FALSE) { + .Call(`_bgms_zratio_test_calibrated_eval`, graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, seed, n_sweep, burn, freeze_after, slab_cauchy) } zratio_test_precompute <- function(G, edges, addc, tg, ihat, ghat, wt, psi0, ncn_max, bre_max) { diff --git a/R/bgms-package.R b/R/bgms-package.R index 7f995cbc..0c6d7ecf 100644 --- a/R/bgms-package.R +++ b/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{} diff --git a/R/build_output_bgm.R b/R/build_output_bgm.R index f2b17eb2..e3a8a4df 100644 --- a/R/build_output_bgm.R +++ b/R/build_output_bgm.R @@ -295,7 +295,8 @@ build_output_bgm = function(spec, raw) { # --- Z-ratio alarm suite (hierarchical graph-prior spec) ---------------------- if(!is.null(zratio_chains)) { zc = zratio_cell_constants( - p$delta, p$pairwise_scale, p$scale_rate, p$scale_eta + 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, diff --git a/R/build_output_mixed_mrf.R b/R/build_output_mixed_mrf.R index 21030093..c6171c74 100644 --- a/R/build_output_mixed_mrf.R +++ b/R/build_output_mixed_mrf.R @@ -290,7 +290,8 @@ build_output_mixed_mrf = function(spec, raw) { # triangles without diagonals; the audit reads the Gyy segment. if(!is.null(zratio_chains)) { zc = zratio_cell_constants( - pr$delta, pr$pairwise_scale, pr$scale_rate, pr$scale_eta + 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, diff --git a/R/run_sampler.R b/R/run_sampler.R index c42388cc..5715591a 100644 --- a/R/run_sampler.R +++ b/R/run_sampler.R @@ -79,12 +79,13 @@ run_sampler_ggm = function(spec) { zratio = NULL if(identical(p$precision_graph_prior, "hierarchical")) { zc = zratio_cell_constants( - p$delta, p$pairwise_scale, p$scale_rate, p$scale_eta + 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, sigma = zc$sigma, beta = zc$beta, slab = zc$slab, calibration_window = resolve_zratio_calibration_window( p$calibration_window, d$num_variables, s$warmup ) @@ -221,12 +222,13 @@ run_sampler_mixed_mrf = function(spec) { zratio = NULL if(identical(p$precision_graph_prior, "hierarchical")) { zc = zratio_cell_constants( - p$delta, p$pairwise_scale, p$scale_rate, p$scale_eta + 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, sigma = zc$sigma, beta = zc$beta, slab = zc$slab, calibration_window = resolve_zratio_calibration_window( p$calibration_window, d$num_continuous, s$warmup ) diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R index 64c7e942..ea960c1a 100644 --- a/R/sample_ggm_prior.R +++ b/R/sample_ggm_prior.R @@ -337,12 +337,13 @@ sample_ggm_prior = function( )) } zc = zratio_cell_constants( - delta, ip$pairwise_scale, sp$scale_rate, sp$scale_eta + 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, sigma = zc$sigma, beta = zc$beta, slab = zc$slab, calibration_window = resolve_zratio_calibration_window( calibration_window, p, n_warmup ) diff --git a/R/zratio_diagnostics.R b/R/zratio_diagnostics.R index 8c2b565a..2c8967d6 100644 --- a/R/zratio_diagnostics.R +++ b/R/zratio_diagnostics.R @@ -278,7 +278,8 @@ zratio_audit_chain = function(chain, zratio_spec, num_nodes, n_graphs, top_k, 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) + as.integer(audit_sweep), 30L, as.integer(seed + g), + identical(zratio_spec$slab, "cauchy") ) err[sel] = audit$err } diff --git a/R/zratio_tables.R b/R/zratio_tables.R index 06461975..092bbc37 100644 --- a/R/zratio_tables.R +++ b/R/zratio_tables.R @@ -22,7 +22,8 @@ # (sd_marginal_helpers.R, ks_validation_grid.R). # Golub-Welsch Gauss quadrature nodes/weights. kind: "laguerre" (weight -# e^{-x} on (0, Inf)), "hermite" (weight e^{-x^2} on (-Inf, Inf)), +# e^{-x} on (0, Inf)), "laguerre_half" (generalized Laguerre, weight +# x^{-1/2} e^{-x} on (0, Inf)), "hermite" (weight e^{-x^2} on (-Inf, Inf)), # "legendre" (weight 1 on (-1, 1)). zratio_gauss_quad = function(n, kind) { k = seq_len(n - 1) @@ -30,6 +31,10 @@ zratio_gauss_quad = function(n, kind) { a = 2 * seq_len(n) - 1 b = k mu0 = 1 + } else if(kind == "laguerre_half") { + a = 2 * seq_len(n) - 1.5 + b = sqrt(k * (k - 0.5)) + mu0 = gamma(0.5) } else if(kind == "hermite") { a = rep(0, n) b = sqrt(k / 2) @@ -54,6 +59,19 @@ zratio_gauss_quad = function(n, kind) { ) } +# Cauchy-slab leg mixture in resolvent form: with omega ~ IG(1/2, 1/2) +# (equivalently omega = 1/(2t), t ~ Gamma(1/2, 1)), +# E[omega^k (A + B omega)^-(k + 1/2)] +# = sqrt(2/pi) sum_i w_i (2 t_i A + B)^-(k + 1/2) +# exactly, on the plain Gauss-Laguerre grid (t_i, w_i). The transformed +# integrand is analytic in t (the omega-side form carries a sqrt(t) factor +# that defeats polynomial rules), and every leg term in the channel +# moments has this (k, k + 1/2) power pairing. +zratio_omega_mixture = function(n) { + q = zratio_gauss_quad(n, "laguerre") + list(t = q$nodes, w = sqrt(2 / pi) * q$weights) +} + # Closed-form spike pair integral # I_spike(c) = 2 Gamma(nu) beta^{-nu} |c|^nu K_nu(2 beta |c|), nu = delta + 1, # with limit Gamma(nu)^2 beta^{-2 nu} at c = 0. @@ -73,9 +91,11 @@ zratio_ispike = function(c_val, delta, beta) { # Slab pair integral G(c) tabulated on [0, cmax] with linear interpolation. # Inner S12 integral by the sin substitution (removes the boundary # singularity of (S11 S22 - S12^2)^delta); outer (S11, S22) by -# Gauss-Laguerre with weight e^{-beta S}. +# Gauss-Laguerre with weight e^{-beta S}. The slab density multiplies at +# the raw entry (the tilt sees the Schur-shifted entry), so the Cauchy +# variant only swaps the density factor. zratio_pair_integrals = function( - delta, sigma, beta, + delta, sigma, beta, slab = "normal", cmax = 18, ngrid = 121, nlag = 48, nleg = 64 ) { gl = zratio_gauss_quad(nlag, "laguerre") @@ -90,9 +110,14 @@ zratio_pair_integrals = function( s2 = rep(xq, times = nlag) w_pair = rep(wq, each = nlag) * rep(wq, times = nlag) b = sqrt(s1 * s2) + slab_dens = if(identical(slab, "cauchy")) { + function(x) dcauchy(x, 0, sigma) + } else { + function(x) dnorm(x, 0, sigma) + } g_one = function(c_val) { cm = outer(b, sth) - dn = dnorm(cm + c_val, 0, sigma) + dn = slab_dens(cm + c_val) inner = (b^(2 * delta + 1)) * as.numeric(dn %*% (wt * cth^(2 * delta + 1))) sum(w_pair * inner) } @@ -133,8 +158,26 @@ zratio_saddle_grid = function( # Two-moment constants for the common-neighbour node channel: weighted # moments of the single-node resolvent under the tilted diagonal prior, # w_k = sigma^{4k} int x^{delta+1} (x + t2)^{-(2k+1)} e^{-beta x} dx / I(1). -zratio_node_channel = function(delta, sigma, beta) { +# Cauchy slab: the two legs from the node to the toggled endpoints carry +# independent mixture weights, K_leg | omega ~ N(0, sigma^2 omega), so the +# resolvent splits per leg, (x + t2)^{-(2k+1)} -> +# (x + t2 w_a)^{-(2k+1)/2} (x + t2 w_b)^{-(2k+1)/2} with prefactor +# (w_a w_b)^k, and each leg mixes on the omega grid. +zratio_node_channel = function(delta, sigma, beta, slab = "normal") { t2 = 2 * beta * sigma^2 + if(identical(slab, "cauchy")) { + mx = zratio_omega_mixture(48) + gl = zratio_gauss_quad(96, "laguerre") + x = gl$nodes / beta + cx = (gl$weights / beta) * x^(delta + 1) + tx = 2 * outer(x, mx$t) + mix = function(p) as.numeric((tx + t2)^(-p) %*% mx$w) + den = sum(cx * mix(0.5)^2) + return(c( + sigma^4 * sum(cx * mix(1.5)^2) / den, + sigma^8 * sum(cx * mix(2.5)^2) / den + )) + } ip = function(p) { integrate( function(x) x^(delta + 1) * (x + t2)^(-p) * exp(-beta * x), @@ -149,35 +192,59 @@ zratio_node_channel = function(delta, sigma, beta) { # Excess two-moment constants for the CN-CN edge channel: moments of the # connected 2-clique block minus twice the single-node constants. Estimated # by a seeded within-block Gibbs run (matches the reference implementation -# draw for draw at the same seed). +# draw for draw at the same seed). Cauchy slab: the block coupling runs +# omega-augmented (conjugate IG(1, 1/2 + k^2/(2 sigma^2)) refresh after +# each draw), the four legs draw fresh prior weights sqrt(omega) = 1/|z| +# per kept sweep, and the moments use the leg-dressed block recipe +# Mi = (I + t2 Wi R Wi)^{-1}, Wt = sqrt(det Mi det Mj), +# P = (Wi Mi Wi) R (Wj Mj Wj) R, p1 = sigma^4 tr P, p2 = sigma^8 tr P^2, +# which reduces to the a2-resolvent form at omega = 1. zratio_clique2_moments = function( - delta, sigma, beta, n_mc = 20000, burn = 60, seed = 7 + delta, sigma, beta, slab = "normal", n_mc = 20000, burn = 60, seed = 7 ) { t2 = 2 * beta * sigma^2 + cauchy = identical(slab, "cauchy") set.seed(seed) k_mat = diag(rexp(2, beta) + 2, 2) s2i = 1 / sigma^2 + om_e = 1 p1 = p2 = w = numeric(n_mc) for(s in 1:(burn + n_mc)) { for(i in 1:2) { rest = setdiff(1:2, i) c_inv = solve(k_mat[rest, rest, drop = FALSE]) m = 2 * beta * c_inv - diag(m) = diag(m) + s2i + diag(m) = diag(m) + if(cauchy) 1 / (sigma^2 * om_e) else s2i r = chol(m) bvec = backsolve(r, rnorm(1)) k_mat[i, rest] = bvec k_mat[rest, i] = bvec k_mat[i, i] = rgamma(1, delta + 1, beta) + as.numeric(t(bvec) %*% c_inv %*% bvec) + if(cauchy) { + om_e = (0.5 + bvec^2 / (2 * sigma^2)) / rexp(1) + } } if(s > burn) { k = s - burn - a2 = k_mat + t2 * diag(2) - ri = solve(a2) - w[k] = det(k_mat) / det(a2) - p1[k] = sigma^4 * sum(diag(ri %*% ri)) - p2[k] = sigma^8 * sum(diag(ri %*% ri %*% ri %*% ri)) + if(cauchy) { + wsi = 1 / abs(rnorm(2)) + wsj = 1 / abs(rnorm(2)) + r_blk = solve(k_mat) + mi = solve(diag(2) + t2 * (r_blk * outer(wsi, wsi))) + mj = solve(diag(2) + t2 * (r_blk * outer(wsj, wsj))) + p_mat = (mi * outer(wsi, wsi)) %*% r_blk %*% + (mj * outer(wsj, wsj)) %*% r_blk + w[k] = sqrt(det(mi) * det(mj)) + p1[k] = sigma^4 * sum(diag(p_mat)) + p2[k] = sigma^8 * sum(p_mat * t(p_mat)) + } else { + a2 = k_mat + t2 * diag(2) + ri = solve(a2) + w[k] = det(k_mat) / det(a2) + p1[k] = sigma^4 * sum(diag(ri %*% ri)) + p2[k] = sigma^8 * sum(diag(ri %*% ri %*% ri %*% ri)) + } } } ok = is.finite(w) & is.finite(p1) & is.finite(p2) @@ -188,12 +255,50 @@ zratio_clique2_moments = function( # Two-moment constants for the bridge channel (edge from Si\Sj to Sj\Si): # 2-node quadrature, Gauss-Laguerre on the diagonals x Gauss-Hermite on the -# coupling. -zratio_bridge_channel = function(delta, sigma, beta, nlag = 64, nher = 80) { +# coupling. Cauchy slab: the coupling integrates on its exact PD support by +# the sin substitution u = sqrt(kaa kbb) sin(theta) against the Cauchy +# density, and the two legs mix independently on the omega grid; at fixed +# theta the leg resolvents factor, Da = d + t2 w_a kbb, Db = d + t2 w_b kaa, +# with per-moment prefactors (w_a w_b)^k folded into the leg mixtures. +zratio_bridge_channel = function(delta, sigma, beta, slab = "normal", + nlag = 64, nher = 80) { t2 = 2 * beta * sigma^2 gl = zratio_gauss_quad(nlag, "laguerre") xa = gl$nodes / beta wa = gl$weights + if(identical(slab, "cauchy")) { + lgq = zratio_gauss_quad(64, "legendre") + th = lgq$nodes * (pi / 2) + wth = lgq$weights * (pi / 2) + sth = sin(th) + cth = cos(th) + mx = zratio_omega_mixture(32) + den = n1 = n2 = 0 + for(ia in 1:nlag) { + kaa = xa[ia] + for(ib in 1:nlag) { + kbb = xa[ib] + pw = wa[ia] * wa[ib] + bmax = sqrt(kaa * kbb) + u = bmax * sth + d = bmax^2 * cth^2 + wu = pw * wth * dcauchy(u, 0, sigma) * bmax * cth + da = 2 * outer(d, mx$t) + t2 * kbb + db = 2 * outer(d, mx$t) + t2 * kaa + fa0 = as.numeric(da^(-0.5) %*% mx$w) + fb0 = as.numeric(db^(-0.5) %*% mx$w) + fa1 = as.numeric(da^(-1.5) %*% mx$w) + fb1 = as.numeric(db^(-1.5) %*% mx$w) + fa2 = as.numeric(da^(-2.5) %*% mx$w) + fb2 = as.numeric(db^(-2.5) %*% mx$w) + base = wu * d^delta * d + den = den + sum(base * fa0 * fb0) + n1 = n1 + sum(base * sigma^4 * u^2 * fa1 * fb1) + n2 = n2 + sum(base * sigma^8 * u^4 * fa2 * fb2) + } + } + return(c(n1 / den, n2 / den)) + } gh = zratio_gauss_quad(nher, "hermite") zk = gh$nodes wh = gh$weights / sqrt(pi) @@ -226,9 +331,9 @@ zratio_bridge_channel = function(delta, sigma, beta, nlag = 64, nher = 80) { # user-specified standardized rate when the scale prior carries one, else # pairwise_scale * scale_rate (the same number up to rounding). zratio_cell_constants = function(delta, pairwise_scale, scale_rate, - scale_eta = NA_real_) { + scale_eta = NA_real_, slab = "normal") { eta = if(is.finite(scale_eta)) scale_eta else pairwise_scale * scale_rate - zratio_constants(delta, sigma = 1, beta = eta) + zratio_constants(delta, sigma = 1, beta = eta, slab = slab) } # Session cache for zratio_constants: the constant set is deterministic per @@ -236,17 +341,18 @@ zratio_cell_constants = function(delta, pairwise_scale, scale_rate, # once (sampler dispatch and diagnostics assembly). zratio_constants_cache = new.env(parent = emptyenv()) -# Full fit-time constant set for one (delta, sigma, beta) cell: +# Full fit-time constant set for one (delta, sigma, beta, slab) cell: # addc[1..6] (R indexing) = (w1, w2, ce1, ce2, cb1, cb2), the saddle grid, # and psi0 = I_spike(0)/G(0). The OLS-correction slots (7..13) and the hull # box (14..23) are absent here; the warm-up calibrator appends them. # The clique-2 channel draws seeded Monte Carlo samples, so the caller's RNG # state is saved and restored. Results are served from a session cache keyed # on the cell. -zratio_constants = function(delta, sigma, beta) { +zratio_constants = function(delta, sigma, beta, slab = "normal") { + slab = match.arg(slab, c("normal", "cauchy")) key = paste( format(delta, digits = 17), format(sigma, digits = 17), - format(beta, digits = 17), + format(beta, digits = 17), slab, sep = "_" ) cached = zratio_constants_cache[[key]] @@ -265,16 +371,23 @@ zratio_constants = function(delta, sigma, beta) { add = TRUE ) } - pair = zratio_pair_integrals(delta, sigma, beta) + pair = zratio_pair_integrals(delta, sigma, beta, slab) grid = zratio_saddle_grid(pair) - w12 = zratio_node_channel(delta, sigma, beta) - e2 = zratio_clique2_moments(delta, sigma, beta) - cb = zratio_bridge_channel(delta, sigma, beta) + w12 = zratio_node_channel(delta, sigma, beta, slab) + # The omega-augmented chain mixes slower than the Normal one, so the + # Cauchy cell runs longer; the Normal cell keeps its draw-for-draw + # reference length. + e2 = zratio_clique2_moments( + delta, sigma, beta, slab, + n_mc = if(identical(slab, "cauchy")) 60000 else 20000 + ) + cb = zratio_bridge_channel(delta, sigma, beta, slab) addc = c(w12[1], w12[2], e2[1] - 2 * w12[1], e2[2] - 2 * w12[2], cb[1], cb[2]) out = list( delta = delta, sigma = sigma, beta = beta, + slab = slab, addc = addc, tg = grid$tg, ihat = grid$ihat, diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 51f00e21..bf239ecc 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -775,8 +775,8 @@ BEGIN_RCPP END_RCPP } // zratio_audit_edges -Rcpp::List zratio_audit_edges(arma::imat G, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double sigma, double beta, int n_sweep, int burn, int seed); -RcppExport SEXP _bgms_zratio_audit_edges(SEXP GSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP sigmaSEXP, SEXP betaSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP seedSEXP) { +Rcpp::List zratio_audit_edges(arma::imat G, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double sigma, double beta, int n_sweep, int burn, int seed, bool slab_cauchy); +RcppExport SEXP _bgms_zratio_audit_edges(SEXP GSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP sigmaSEXP, SEXP betaSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP seedSEXP, SEXP slab_cauchySEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -794,7 +794,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type n_sweep(n_sweepSEXP); Rcpp::traits::input_parameter< int >::type burn(burnSEXP); Rcpp::traits::input_parameter< int >::type seed(seedSEXP); - rcpp_result_gen = Rcpp::wrap(zratio_audit_edges(G, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, n_sweep, burn, seed)); + Rcpp::traits::input_parameter< bool >::type slab_cauchy(slab_cauchySEXP); + rcpp_result_gen = Rcpp::wrap(zratio_audit_edges(G, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, n_sweep, burn, seed, slab_cauchy)); return rcpp_result_gen; END_RCPP } @@ -835,8 +836,8 @@ BEGIN_RCPP END_RCPP } // zratio_test_calibrated_eval -Rcpp::List zratio_test_calibrated_eval(Rcpp::List graphs, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double sigma, double beta, int seed, int n_sweep, int burn, int freeze_after); -RcppExport SEXP _bgms_zratio_test_calibrated_eval(SEXP graphsSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP sigmaSEXP, SEXP betaSEXP, SEXP seedSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP freeze_afterSEXP) { +Rcpp::List zratio_test_calibrated_eval(Rcpp::List graphs, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double sigma, double beta, int seed, int n_sweep, int burn, int freeze_after, bool slab_cauchy); +RcppExport SEXP _bgms_zratio_test_calibrated_eval(SEXP graphsSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP sigmaSEXP, SEXP betaSEXP, SEXP seedSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP freeze_afterSEXP, SEXP slab_cauchySEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -855,7 +856,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type n_sweep(n_sweepSEXP); Rcpp::traits::input_parameter< int >::type burn(burnSEXP); Rcpp::traits::input_parameter< int >::type freeze_after(freeze_afterSEXP); - rcpp_result_gen = Rcpp::wrap(zratio_test_calibrated_eval(graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, seed, n_sweep, burn, freeze_after)); + Rcpp::traits::input_parameter< bool >::type slab_cauchy(slab_cauchySEXP); + rcpp_result_gen = Rcpp::wrap(zratio_test_calibrated_eval(graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, seed, n_sweep, burn, freeze_after, slab_cauchy)); return rcpp_result_gen; END_RCPP } @@ -921,10 +923,10 @@ static const R_CallMethodDef CallEntries[] = { {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4}, {"_bgms_test_warmup_schedule", (DL_FUNC) &_bgms_test_warmup_schedule, 5}, {"_bgms_zratio_scan_graph", (DL_FUNC) &_bgms_zratio_scan_graph, 7}, - {"_bgms_zratio_audit_edges", (DL_FUNC) &_bgms_zratio_audit_edges, 14}, + {"_bgms_zratio_audit_edges", (DL_FUNC) &_bgms_zratio_audit_edges, 15}, {"_bgms_zratio_test_eval", (DL_FUNC) &_bgms_zratio_test_eval, 8}, {"_bgms_zratio_test_saddle", (DL_FUNC) &_bgms_zratio_test_saddle, 8}, - {"_bgms_zratio_test_calibrated_eval", (DL_FUNC) &_bgms_zratio_test_calibrated_eval, 15}, + {"_bgms_zratio_test_calibrated_eval", (DL_FUNC) &_bgms_zratio_test_calibrated_eval, 16}, {"_bgms_zratio_test_precompute", (DL_FUNC) &_bgms_zratio_test_precompute, 10}, {NULL, NULL, 0} }; diff --git a/src/models/ggm/zratio_engine.cpp b/src/models/ggm/zratio_engine.cpp index 5bc51ae2..c270d02b 100644 --- a/src/models/ggm/zratio_engine.cpp +++ b/src/models/ggm/zratio_engine.cpp @@ -247,7 +247,8 @@ bool ZRatioEngine::audit_edge(const arma::imat& G, int i, int j, void ZRatioEngine::enable_calibration(double delta, double sigma, double beta, SafeRNG* rng, int n_sweep, int burn, - double maha_thresh, int min_anchors) { + double maha_thresh, int min_anchors, + bool slab_cauchy) { calibration_enabled_ = true; frozen_ = false; delta_ = delta; @@ -258,6 +259,7 @@ void ZRatioEngine::enable_calibration(double delta, double sigma, double beta, burn_ = burn; maha_thresh_ = maha_thresh; min_anchors_ = min_anchors; + oracle_slab_cauchy_ = slab_cauchy; } void ZRatioEngine::refit_() { @@ -290,7 +292,7 @@ void ZRatioEngine::freeze_calibration() { addc_ = packed; } -void ZRatioEngine::gibbs_sweep_(arma::mat& k_blk, +void ZRatioEngine::gibbs_sweep_(arma::mat& k_blk, arma::mat& omega_blk, const std::vector& nbr) const { const int m = static_cast(k_blk.n_rows); const double s2i = 1.0 / (sigma_ * sigma_); @@ -312,7 +314,13 @@ void ZRatioEngine::gibbs_sweep_(arma::mat& k_blk, } arma::mat c_mat = a_inv.submat(idx_a, idx_a); arma::mat m_mat = 2.0 * beta_ * c_mat; - m_mat.diag() += s2i; + if (oracle_slab_cauchy_) { + for (int j = 0; j < nq; ++j) { + m_mat(j, j) += s2i / omega_blk(i, ni[j]); + } + } else { + m_mat.diag() += s2i; + } arma::mat r_chol; if (!arma::chol(r_chol, m_mat)) continue; arma::vec z(nq); @@ -325,6 +333,16 @@ void ZRatioEngine::gibbs_sweep_(arma::mat& k_blk, k_blk(i, ni[j]) = bvec[j]; } k_blk(i, i) = xi + quad; + if (oracle_slab_cauchy_) { + // Conjugate omega | k ~ IG(1, 1/2 + k^2 / (2 sigma^2)). + for (int j = 0; j < nq; ++j) { + const double b = bvec[j]; + const double ig_rate = 0.5 + 0.5 * b * b * s2i; + const double wo = ig_rate / rexp(*rng_, 1.0); + omega_blk(i, ni[j]) = wo; + omega_blk(ni[j], i) = wo; + } + } } else { k_blk(i, i) = rgamma(*rng_, delta_ + 1.0, beta_); } @@ -332,13 +350,36 @@ void ZRatioEngine::gibbs_sweep_(arma::mat& k_blk, } bool ZRatioEngine::inner_moments_(const arma::mat& k_blk, const arma::uvec& si, - const arma::uvec& sj, double& w, double& p1, + const arma::uvec& sj, const arma::vec& wsi, + const arma::vec& wsj, double& w, double& p1, double& p2) const { const double t2 = 2.0 * beta_ * sigma_ * sigma_; arma::mat r_inv; if (!arma::inv_sympd(r_inv, k_blk)) return false; arma::mat rii = r_inv.submat(si, si), rjj = r_inv.submat(sj, sj), rij = r_inv.submat(si, sj); + const double s4 = sigma_ * sigma_ * sigma_ * sigma_, s8 = s4 * s4; + if (oracle_slab_cauchy_) { + // Leg-dressed recipe under the scale-mixture slab: with + // Wi = diag(sqrt(omega_leg)), Mi = (I + t2 Wi Rii Wi)^{-1} and + // P = (Wi Mi Wi) Rij (Wj Mj Wj) Rij^T; reduces to the plain + // resolvent form at omega = 1. + arma::mat di = arma::diagmat(wsi), dj = arma::diagmat(wsj); + arma::mat mi, mj; + if (!arma::inv_sympd(mi, arma::eye(si.n_elem, si.n_elem) + + t2 * di * rii * di)) { + return false; + } + if (!arma::inv_sympd(mj, arma::eye(sj.n_elem, sj.n_elem) + + t2 * dj * rjj * dj)) { + return false; + } + arma::mat p_mat = (di * mi * di) * rij * (dj * mj * dj) * rij.t(); + w = std::sqrt(arma::det(mi) * arma::det(mj)); + p1 = s4 * arma::trace(p_mat); + p2 = s8 * arma::accu(p_mat % p_mat.t()); + return true; + } arma::mat mi, mj; if (!arma::inv_sympd(mi, arma::eye(si.n_elem, si.n_elem) + t2 * rii)) { return false; @@ -350,7 +391,6 @@ bool ZRatioEngine::inner_moments_(const arma::mat& k_blk, const arma::uvec& si, // values of Mi^.5 Rij Mj^.5, sum u_k = sigma^4 tr(P) and sum u_k^2 = // sigma^8 tr(P^2) for P = Mi Rij Mj Rij^T. arma::mat p_mat = mi * rij * mj * rij.t(); - const double s4 = sigma_ * sigma_ * sigma_ * sigma_, s8 = s4 * s4; w = std::sqrt(arma::det(mi) * arma::det(mj)); p1 = s4 * arma::trace(p_mat); p2 = s8 * arma::accu(p_mat % p_mat.t()); @@ -371,16 +411,31 @@ bool ZRatioEngine::block_oracle_moments(const arma::imat& a_blk, nbr[i] = arma::uvec(v); } arma::mat k_blk(m, m, arma::fill::zeros); + arma::mat omega_blk; + if (oracle_slab_cauchy_) omega_blk.ones(m, m); for (int l = 0; l < m; ++l) { k_blk(l, l) = rexp(*rng_, beta_) + m; } - for (int s = 0; s < burn_; ++s) gibbs_sweep_(k_blk, nbr); + for (int s = 0; s < burn_; ++s) gibbs_sweep_(k_blk, omega_blk, nbr); double sw = 0, sw1 = 0, sw2 = 0; long kept = 0; + arma::vec wsi, wsj; for (int s = 0; s < n_sweep_; ++s) { - gibbs_sweep_(k_blk, nbr); + gibbs_sweep_(k_blk, omega_blk, nbr); + if (oracle_slab_cauchy_) { + // Fresh leg weights per kept sweep: sqrt(omega) = 1/|z| with z + // standard normal; the sweep average carries the mixture. + wsi.set_size(si.n_elem); + wsj.set_size(sj.n_elem); + for (arma::uword k = 0; k < wsi.n_elem; ++k) { + wsi[k] = 1.0 / std::abs(rnorm(*rng_, 0.0, 1.0)); + } + for (arma::uword k = 0; k < wsj.n_elem; ++k) { + wsj[k] = 1.0 / std::abs(rnorm(*rng_, 0.0, 1.0)); + } + } double w, p1, p2; - if (inner_moments_(k_blk, si, sj, w, p1, p2)) { + if (inner_moments_(k_blk, si, sj, wsi, wsj, w, p1, p2)) { sw += w; sw1 += w * p1; sw2 += w * p2; diff --git a/src/models/ggm/zratio_engine.h b/src/models/ggm/zratio_engine.h index 92c27d7c..02f8cf8d 100644 --- a/src/models/ggm/zratio_engine.h +++ b/src/models/ggm/zratio_engine.h @@ -119,16 +119,21 @@ class ZRatioEngine { /** * Set the standardized-cell prior constants and RNG the block-Gibbs * oracle samples under, without entering calibration mode. rng must - * outlive the engine. + * outlive the engine. slab_cauchy selects the Cauchy slab family: the + * block couplings then run omega-augmented (scale-mixture of normals) + * and the endpoint legs mix per sweep, matching the marginal-Cauchy + * normalizer the tables integrate. */ void set_oracle_params(double delta, double sigma, double beta, - SafeRNG* rng, int n_sweep = 300, int burn = 30) { + SafeRNG* rng, int n_sweep = 300, int burn = 30, + bool slab_cauchy = false) { delta_ = delta; sigma_ = sigma; beta_ = beta; rng_ = rng; n_sweep_ = n_sweep; burn_ = burn; + oracle_slab_cauchy_ = slab_cauchy; } /** @@ -144,11 +149,13 @@ class ZRatioEngine { * * (delta, sigma, beta) are the standardized-cell prior constants the * oracle samples under; rng must outlive the engine (the model's - * chain RNG). + * chain RNG). slab_cauchy selects the Cauchy slab family for the + * oracle (see set_oracle_params). */ void enable_calibration(double delta, double sigma, double beta, SafeRNG* rng, int n_sweep = 300, int burn = 30, - double maha_thresh = 9.0, int min_anchors = 6); + double maha_thresh = 9.0, int min_anchors = 6, + bool slab_cauchy = false); /** Refit and freeze: pack coefficients + hull box into addc[6..22]. */ void freeze_calibration(); @@ -183,10 +190,11 @@ class ZRatioEngine { const arma::vec& anchors_y() const { return ay_; } private: - void gibbs_sweep_(arma::mat& k_blk, + void gibbs_sweep_(arma::mat& k_blk, arma::mat& omega_blk, const std::vector& nbr) const; bool inner_moments_(const arma::mat& k_blk, const arma::uvec& si, - const arma::uvec& sj, double& w, double& p1, + const arma::uvec& sj, const arma::vec& wsi, + const arma::vec& wsj, double& w, double& p1, double& p2) const; void refit_(); @@ -207,6 +215,7 @@ class ZRatioEngine { // Online-calibration state (inert unless enable_calibration ran). bool calibration_enabled_ = false; bool frozen_ = true; + bool oracle_slab_cauchy_ = false; double delta_ = 0.0, sigma_ = 1.0, beta_ = 0.5; SafeRNG* rng_ = nullptr; int n_sweep_ = 300, burn_ = 30; diff --git a/src/sample_ggm.cpp b/src/sample_ggm.cpp index d32cfca3..7b465c86 100644 --- a/src/sample_ggm.cpp +++ b/src/sample_ggm.cpp @@ -124,11 +124,14 @@ Rcpp::List sample_ggm( zratio_window = Rcpp::as(zs["calibration_window"]); } if (zratio_window > 0) { + bool zr_cauchy = zs.containsElementNamed("slab") && + Rcpp::as(zs["slab"]) == "cauchy"; // The rng pointer is rebound per chain clone by GGMModel. engine->enable_calibration( Rcpp::as(zs["delta"]), Rcpp::as(zs["sigma"]), - Rcpp::as(zs["beta"]), nullptr); + Rcpp::as(zs["beta"]), nullptr, 300, 30, 9.0, 6, + zr_cauchy); } model.set_zratio_engine(std::move(engine)); } diff --git a/src/sample_mixed.cpp b/src/sample_mixed.cpp index 3b1d1336..10c51d6a 100644 --- a/src/sample_mixed.cpp +++ b/src/sample_mixed.cpp @@ -182,11 +182,14 @@ Rcpp::List sample_mixed_mrf( zratio_window = Rcpp::as(zs["calibration_window"]); } if (zratio_window > 0) { + bool zr_cauchy = zs.containsElementNamed("slab") && + Rcpp::as(zs["slab"]) == "cauchy"; // The rng pointer is rebound per chain clone by MixedMRFModel. engine->enable_calibration( Rcpp::as(zs["delta"]), Rcpp::as(zs["sigma"]), - Rcpp::as(zs["beta"]), nullptr); + Rcpp::as(zs["beta"]), nullptr, 300, 30, 9.0, 6, + zr_cauchy); } model.set_zratio_engine(std::move(engine)); } diff --git a/src/zratio_audit_interface.cpp b/src/zratio_audit_interface.cpp index 536ccb00..8402b44c 100644 --- a/src/zratio_audit_interface.cpp +++ b/src/zratio_audit_interface.cpp @@ -73,7 +73,8 @@ arma::mat zratio_scan_graph( // the deployed correction under `addc` versus the block-Gibbs local // oracle discrepancy log(saddle on oracle moments) - log(additive), and // the audit score |pred - oracle|. (delta, sigma, beta) are the -// bare-scale prior constants. Rows with ok = 0 had a one-sided block, +// standardized-cell prior constants; slab_cauchy selects the Cauchy +// slab family for the oracle. Rows with ok = 0 had a one-sided block, // non-positive additive moments, or an oracle with no finite sweep; // their scores are NA. // ----------------------------------------------------------------------------- @@ -93,11 +94,13 @@ Rcpp::List zratio_audit_edges( double beta, int n_sweep, int burn, - int seed + int seed, + bool slab_cauchy = false ) { ZRatioEngine engine(addc, tg, ihat, ghat, wt, psi0); SafeRNG rng(seed); - engine.set_oracle_params(delta, sigma, beta, &rng, n_sweep, burn); + engine.set_oracle_params(delta, sigma, beta, &rng, n_sweep, burn, + slab_cauchy); const arma::uword n = edges.n_rows; arma::vec err(n), pred(n), oracle(n); arma::ivec ok(n); diff --git a/src/zratio_test_interface.cpp b/src/zratio_test_interface.cpp index 8e2a39d4..64208cfb 100644 --- a/src/zratio_test_interface.cpp +++ b/src/zratio_test_interface.cpp @@ -95,11 +95,13 @@ Rcpp::List zratio_test_calibrated_eval( int seed, int n_sweep, int burn, - int freeze_after + int freeze_after, + bool slab_cauchy = false ) { ZRatioEngine engine(addc, tg, ihat, ghat, wt, psi0); SafeRNG rng(seed); - engine.enable_calibration(delta, sigma, beta, &rng, n_sweep, burn); + engine.enable_calibration(delta, sigma, beta, &rng, n_sweep, burn, 9.0, 6, + slab_cauchy); arma::vec out(edges.n_rows); for (arma::uword e = 0; e < edges.n_rows; ++e) { if (freeze_after > 0 && diff --git a/tests/testthat/test-hier-zratio-identity.R b/tests/testthat/test-hier-zratio-identity.R index ce6bacb3..78008f3f 100644 --- a/tests/testthat/test-hier-zratio-identity.R +++ b/tests/testthat/test-hier-zratio-identity.R @@ -64,7 +64,8 @@ test_that("Z-ratio constants build in the standardized cell", { a = bgms:::zratio_cell_constants(d, pairwise_scale = 0.5, scale_rate = 2) b = bgms:::zratio_cell_constants(d, pairwise_scale = 0.25, scale_rate = 4) e = bgms:::zratio_cell_constants( - d, pairwise_scale = 2.5, scale_rate = 0.4, scale_eta = 1 + d, + pairwise_scale = 2.5, scale_rate = 0.4, scale_eta = 1 ) expect_identical(a, b) expect_identical(a, e) @@ -88,6 +89,22 @@ test_that("hierarchical graph marginal holds at a non-unit slab scale", { expect_lt(abs(mean(d$edge_indicators) - 0.3), 0.02) }) +test_that("hierarchical graph marginal holds for the Cauchy slab", { + skip_on_cran() + # The Cauchy slab needs its own (marginal-Cauchy) Z-ratio constants; with + # the Normal tables this cell read 0.247 for a 0.30 edge prior. + for(um in c("adaptive-metropolis", "gibbs")) { + d = sample_ggm_prior( + p = 6L, n_samples = 6000L, n_warmup = 1500L, + interaction_prior = cauchy_prior(scale = 0.5), + precision_scale_prior = gamma_prior(shape = 1, rate = 2), + spec = "hierarchical", edge_inclusion_prob = 0.3, + update_method = um, delta = 0.5 * log(6), seed = 7L, verbose = FALSE + ) + expect_lt(abs(mean(d$edge_indicators) - 0.3), 0.02, label = um) + } +}) + test_that("hierarchical BB identity: theta ~ Beta(a, b), PIP = a/(a+b)", { skip_on_cran() d = hier_prior_run( diff --git a/tests/testthat/test-zratio-cauchy.R b/tests/testthat/test-zratio-cauchy.R new file mode 100644 index 00000000..d28e2911 --- /dev/null +++ b/tests/testthat/test-zratio-cauchy.R @@ -0,0 +1,116 @@ +# Cauchy-slab Z-ratio constants. The hierarchical spec with a Cauchy +# interaction prior normalizes by the marginal-Cauchy per-graph constant +# Z_C(Gamma) = E_omega[Z_N(Gamma; sigma sqrt(omega))]: the omega updates +# stay conjugate, the between-move slab factor stays omega-conditional, and +# the Z-ratio becomes a fixed function of Gamma with slab-family-specific +# channel constants. Each quadrature channel is checked here against a +# direct Monte Carlo of the same F-measure moment ratio; the graph-law +# identity itself is gated in test-hier-zratio-identity.R. + +test_that("Cauchy node channel matches direct Monte Carlo", { + skip_on_cran() + delta = 0.5 * log(6) + sigma = 1 + beta = 1 + t2 = 2 * beta * sigma^2 + w_quad = bgms:::zratio_node_channel(delta, sigma, beta, slab = "cauchy") + set.seed(101) + n = 1e6 + x = rgamma(n, delta + 2, rate = beta) + oa = 1 / rnorm(n)^2 + ob = 1 / rnorm(n)^2 + da = x + t2 * oa + db = x + t2 * ob + den = mean((da * db)^-0.5) + w1 = sigma^4 * mean(oa * ob * (da * db)^-1.5) / den + w2 = sigma^8 * mean((oa * ob)^2 * (da * db)^-2.5) / den + expect_lt(abs(w1 / w_quad[1] - 1), 0.01) + expect_lt(abs(w2 / w_quad[2] - 1), 0.02) +}) + +test_that("Cauchy bridge channel matches direct Monte Carlo", { + skip_on_cran() + delta = 0.5 * log(6) + sigma = 1 + beta = 1 + t2 = 2 * beta * sigma^2 + cb_quad = bgms:::zratio_bridge_channel(delta, sigma, beta, slab = "cauchy") + set.seed(202) + n = 2e6 + kaa = rexp(n, beta) + kbb = rexp(n, beta) + u = rcauchy(n, 0, sigma) + oa = 1 / rnorm(n)^2 + ob = 1 / rnorm(n)^2 + d = kaa * kbb - u^2 + pd = d > 0 + d = d[pd] + da = d + t2 * oa[pd] * kbb[pd] + db = d + t2 * ob[pd] * kaa[pd] + base = d^delta * d + den = mean(base * (da * db)^-0.5) + cb1 = mean(base * sigma^4 * oa[pd] * ob[pd] * u[pd]^2 * (da * db)^-1.5) / den + cb2 = mean( + base * sigma^8 * (oa[pd] * ob[pd])^2 * u[pd]^4 * (da * db)^-2.5 + ) / den + expect_lt(abs(cb1 / cb_quad[1] - 1), 0.02) + expect_lt(abs(cb2 / cb_quad[2] - 1), 0.04) +}) + +test_that("Cauchy isolated-edge ratio psi0 matches direct Monte Carlo", { + skip_on_cran() + delta = 0.5 * log(6) + pair = bgms:::zratio_pair_integrals(delta, + sigma = 1, beta = 1, + slab = "cauchy" + ) + psi0 = pair$ispike(0) / pair$g(0) + set.seed(303) + n = 2e6 + x1 = rexp(n, 1) + x2 = rexp(n, 1) + k = rcauchy(n, 0, 1) + i0 = mean((x1 * x2)^delta) + g0 = mean(ifelse(k^2 < x1 * x2, (x1 * x2 - k^2)^delta, 0)) + expect_lt(abs((i0 / g0) / psi0 - 1), 0.01) +}) + +test_that("Cauchy constants live in their own cache cell", { + skip_on_cran() + d = 0.5 * log(6) + zn = bgms:::zratio_cell_constants(d, 0.5, 2) + zc = bgms:::zratio_cell_constants(d, 0.5, 2, slab = "cauchy") + expect_identical(zn$slab, "normal") + expect_identical(zc$slab, "cauchy") + expect_false(isTRUE(all.equal(zn$addc, zc$addc))) + expect_false(isTRUE(all.equal(zn$psi0, zc$psi0))) + # Scale invariance holds per family: same eta, different frames. + zc2 = bgms:::zratio_cell_constants(d, 0.25, 4, slab = "cauchy") + expect_identical(zc, zc2) +}) + +test_that("Cauchy block-Gibbs oracle tracks the additive prediction", { + skip_on_cran() + # Two common neighbours, no CN-CN edge: the additive form is the exact + # no-edge baseline (kappa_2 = m * w1), so the oracle correction + # log(saddle on oracle moments) - log(additive) must sit at MC noise. + # Exercises the omega-augmented sweep and the leg-dressed moments. + q = 4L + G = matrix(0L, q, q) + diag(G) = 1L + for(e in list(c(1, 3), c(2, 3), c(1, 4), c(2, 4))) { + G[e[1], e[2]] = 1L + G[e[2], e[1]] = 1L + } + for(slab in c("normal", "cauchy")) { + zc = bgms:::zratio_cell_constants(0.5 * log(6), 0.5, 2, slab = slab) + audit = zratio_audit_edges( + G, matrix(c(1L, 2L), 1, 2), + zc$addc, zc$tg, zc$ihat, zc$ghat, zc$wt, zc$psi0, + zc$delta, zc$sigma, zc$beta, 2000L, 50L, 42L, + identical(slab, "cauchy") + ) + expect_identical(as.integer(audit$ok), 1L, label = slab) + expect_lt(abs(audit$oracle), 0.06, label = slab) + } +}) From 286072bf2113e4b53413a8c2cc80777bc44294c3 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Tue, 7 Jul 2026 08:49:41 +0200 Subject: [PATCH 3/4] Collapse the Z-ratio signature to (delta, eta); drop dead integrator The hierarchical Z-ratio depends only on (delta, eta): the diagonal congruence factors the slab scale out of every between-graph ratio, and the quadrature grids are sized for the standardized (unit slab scale) frame. zratio_constants() now takes (delta, eta) and fixes sigma = 1 internally, so an off-frame build is unrepresentable rather than merely avoided by convention (Fix A already pinned every call site). The spec list, the C++ oracle params (set_oracle_params/enable_calibration), and the audit/test interfaces carry eta in place of (sigma, beta). The reference-builder parity test now matches the two sigma = 1 fixture cells (the sigma = 2 cells are the same eta = 1 physical cell up to grid error and are still covered by the engine-level reproduction test); the online-calibrator test runs at eta = 1 and still beats additive-only. Also removes the orphaned laguerre_half quadrature kind, left unused when the omega mixture was rewritten to the exact resolvent-form rule on the plain Gauss-Laguerre grid. --- R/RcppExports.R | 8 ++-- R/run_sampler.R | 4 +- R/sample_ggm_prior.R | 2 +- R/zratio_diagnostics.R | 21 +++++----- R/zratio_tables.R | 45 ++++++++++------------ src/RcppExports.cpp | 22 +++++------ src/models/ggm/zratio_engine.cpp | 8 ++-- src/models/ggm/zratio_engine.h | 31 +++++++-------- src/sample_ggm.cpp | 3 +- src/sample_mixed.cpp | 3 +- src/zratio_audit_interface.cpp | 12 +++--- src/zratio_test_interface.cpp | 5 +-- tests/testthat/test-hier-zratio-identity.R | 3 +- tests/testthat/test-zratio-alarms.R | 4 +- tests/testthat/test-zratio-cauchy.R | 2 +- tests/testthat/test-zratio-engine.R | 15 ++++---- 16 files changed, 90 insertions(+), 98 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 30477814..40e76921 100644 --- a/R/RcppExports.R +++ b/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, slab_cauchy = FALSE) { - .Call(`_bgms_zratio_audit_edges`, G, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, n_sweep, burn, seed, slab_cauchy) +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, slab_cauchy = FALSE) { - .Call(`_bgms_zratio_test_calibrated_eval`, graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, seed, n_sweep, burn, freeze_after, slab_cauchy) +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) { diff --git a/R/run_sampler.R b/R/run_sampler.R index 5715591a..4b161bd4 100644 --- a/R/run_sampler.R +++ b/R/run_sampler.R @@ -85,7 +85,7 @@ run_sampler_ggm = function(spec) { 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, slab = zc$slab, + delta = zc$delta, eta = zc$eta, slab = zc$slab, calibration_window = resolve_zratio_calibration_window( p$calibration_window, d$num_variables, s$warmup ) @@ -228,7 +228,7 @@ run_sampler_mixed_mrf = function(spec) { 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, slab = zc$slab, + delta = zc$delta, eta = zc$eta, slab = zc$slab, calibration_window = resolve_zratio_calibration_window( p$calibration_window, d$num_continuous, s$warmup ) diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R index ea960c1a..7fa86f0f 100644 --- a/R/sample_ggm_prior.R +++ b/R/sample_ggm_prior.R @@ -343,7 +343,7 @@ sample_ggm_prior = function( 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, slab = zc$slab, + delta = zc$delta, eta = zc$eta, slab = zc$slab, calibration_window = resolve_zratio_calibration_window( calibration_window, p, n_warmup ) diff --git a/R/zratio_diagnostics.R b/R/zratio_diagnostics.R index 2c8967d6..ae03a2b3 100644 --- a/R/zratio_diagnostics.R +++ b/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,7 +277,7 @@ 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, + zratio_spec$delta, zratio_spec$eta, as.integer(audit_sweep), 30L, as.integer(seed + g), identical(zratio_spec$slab, "cauchy") ) @@ -319,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. @@ -332,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 @@ -362,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 diff --git a/R/zratio_tables.R b/R/zratio_tables.R index 092bbc37..a87da148 100644 --- a/R/zratio_tables.R +++ b/R/zratio_tables.R @@ -22,8 +22,7 @@ # (sd_marginal_helpers.R, ks_validation_grid.R). # Golub-Welsch Gauss quadrature nodes/weights. kind: "laguerre" (weight -# e^{-x} on (0, Inf)), "laguerre_half" (generalized Laguerre, weight -# x^{-1/2} e^{-x} on (0, Inf)), "hermite" (weight e^{-x^2} on (-Inf, Inf)), +# e^{-x} on (0, Inf)), "hermite" (weight e^{-x^2} on (-Inf, Inf)), # "legendre" (weight 1 on (-1, 1)). zratio_gauss_quad = function(n, kind) { k = seq_len(n - 1) @@ -31,10 +30,6 @@ zratio_gauss_quad = function(n, kind) { a = 2 * seq_len(n) - 1 b = k mu0 = 1 - } else if(kind == "laguerre_half") { - a = 2 * seq_len(n) - 1.5 - b = sqrt(k * (k - 0.5)) - mu0 = gamma(0.5) } else if(kind == "hermite") { a = rep(0, n) b = sqrt(k / 2) @@ -325,34 +320,37 @@ zratio_bridge_channel = function(delta, sigma, beta, slab = "normal", } # Standardized cell for one fit's Z-ratio constants. The between-graph -# ratio depends on (delta, eta) only, so the constants are built at -# sigma = 1, beta = eta, the frame the quadrature grids are sized for; user -# scale choices with the same eta share one constant set. eta is the -# user-specified standardized rate when the scale prior carries one, else -# pairwise_scale * scale_rate (the same number up to rounding). +# ratio depends on (delta, eta) only, so the constants take just those two: +# eta is the user-specified standardized rate when the scale prior carries +# one, else pairwise_scale * scale_rate (the same number up to rounding). zratio_cell_constants = function(delta, pairwise_scale, scale_rate, scale_eta = NA_real_, slab = "normal") { eta = if(is.finite(scale_eta)) scale_eta else pairwise_scale * scale_rate - zratio_constants(delta, sigma = 1, beta = eta, slab = slab) + zratio_constants(delta, eta, slab = slab) } # Session cache for zratio_constants: the constant set is deterministic per -# (delta, sigma, beta) cell, and one fit resolves the same cell more than -# once (sampler dispatch and diagnostics assembly). +# (delta, eta, slab) cell, and one fit resolves the same cell more than once +# (sampler dispatch and diagnostics assembly). zratio_constants_cache = new.env(parent = emptyenv()) -# Full fit-time constant set for one (delta, sigma, beta, slab) cell: -# addc[1..6] (R indexing) = (w1, w2, ce1, ce2, cb1, cb2), the saddle grid, -# and psi0 = I_spike(0)/G(0). The OLS-correction slots (7..13) and the hull -# box (14..23) are absent here; the warm-up calibrator appends them. -# The clique-2 channel draws seeded Monte Carlo samples, so the caller's RNG +# Full fit-time constant set for one (delta, eta, slab) cell. The estimator +# is defined in the standardized frame, so the slab has unit scale and the +# diagonal rate is eta: the internal channel builders are evaluated at +# sigma = 1, beta = eta (the frame the quadrature grids are sized for), and +# no other scale is representable. addc[1..6] (R indexing) = +# (w1, w2, ce1, ce2, cb1, cb2), plus the saddle grid and +# psi0 = I_spike(0)/G(0). The OLS-correction slots (7..13) and the hull box +# (14..23) are absent here; the warm-up calibrator appends them. The +# clique-2 channel draws seeded Monte Carlo samples, so the caller's RNG # state is saved and restored. Results are served from a session cache keyed # on the cell. -zratio_constants = function(delta, sigma, beta, slab = "normal") { +zratio_constants = function(delta, eta, slab = "normal") { slab = match.arg(slab, c("normal", "cauchy")) + sigma = 1 + beta = eta key = paste( - format(delta, digits = 17), format(sigma, digits = 17), - format(beta, digits = 17), slab, + format(delta, digits = 17), format(eta, digits = 17), slab, sep = "_" ) cached = zratio_constants_cache[[key]] @@ -385,8 +383,7 @@ zratio_constants = function(delta, sigma, beta, slab = "normal") { addc = c(w12[1], w12[2], e2[1] - 2 * w12[1], e2[2] - 2 * w12[2], cb[1], cb[2]) out = list( delta = delta, - sigma = sigma, - beta = beta, + eta = eta, slab = slab, addc = addc, tg = grid$tg, diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index bf239ecc..de1e2586 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -775,8 +775,8 @@ BEGIN_RCPP END_RCPP } // zratio_audit_edges -Rcpp::List zratio_audit_edges(arma::imat G, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double sigma, double beta, int n_sweep, int burn, int seed, bool slab_cauchy); -RcppExport SEXP _bgms_zratio_audit_edges(SEXP GSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP sigmaSEXP, SEXP betaSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP seedSEXP, SEXP slab_cauchySEXP) { +Rcpp::List zratio_audit_edges(arma::imat G, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double eta, int n_sweep, int burn, int seed, bool slab_cauchy); +RcppExport SEXP _bgms_zratio_audit_edges(SEXP GSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP etaSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP seedSEXP, SEXP slab_cauchySEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -789,13 +789,12 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::vec >::type wt(wtSEXP); Rcpp::traits::input_parameter< double >::type psi0(psi0SEXP); Rcpp::traits::input_parameter< double >::type delta(deltaSEXP); - Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP); - Rcpp::traits::input_parameter< double >::type beta(betaSEXP); + Rcpp::traits::input_parameter< double >::type eta(etaSEXP); Rcpp::traits::input_parameter< int >::type n_sweep(n_sweepSEXP); Rcpp::traits::input_parameter< int >::type burn(burnSEXP); Rcpp::traits::input_parameter< int >::type seed(seedSEXP); Rcpp::traits::input_parameter< bool >::type slab_cauchy(slab_cauchySEXP); - rcpp_result_gen = Rcpp::wrap(zratio_audit_edges(G, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, n_sweep, burn, seed, slab_cauchy)); + rcpp_result_gen = Rcpp::wrap(zratio_audit_edges(G, edges, addc, tg, ihat, ghat, wt, psi0, delta, eta, n_sweep, burn, seed, slab_cauchy)); return rcpp_result_gen; END_RCPP } @@ -836,8 +835,8 @@ BEGIN_RCPP END_RCPP } // zratio_test_calibrated_eval -Rcpp::List zratio_test_calibrated_eval(Rcpp::List graphs, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double sigma, double beta, int seed, int n_sweep, int burn, int freeze_after, bool slab_cauchy); -RcppExport SEXP _bgms_zratio_test_calibrated_eval(SEXP graphsSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP sigmaSEXP, SEXP betaSEXP, SEXP seedSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP freeze_afterSEXP, SEXP slab_cauchySEXP) { +Rcpp::List zratio_test_calibrated_eval(Rcpp::List graphs, arma::imat edges, arma::vec addc, arma::vec tg, arma::vec ihat, arma::vec ghat, arma::vec wt, double psi0, double delta, double eta, int seed, int n_sweep, int burn, int freeze_after, bool slab_cauchy); +RcppExport SEXP _bgms_zratio_test_calibrated_eval(SEXP graphsSEXP, SEXP edgesSEXP, SEXP addcSEXP, SEXP tgSEXP, SEXP ihatSEXP, SEXP ghatSEXP, SEXP wtSEXP, SEXP psi0SEXP, SEXP deltaSEXP, SEXP etaSEXP, SEXP seedSEXP, SEXP n_sweepSEXP, SEXP burnSEXP, SEXP freeze_afterSEXP, SEXP slab_cauchySEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -850,14 +849,13 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::vec >::type wt(wtSEXP); Rcpp::traits::input_parameter< double >::type psi0(psi0SEXP); Rcpp::traits::input_parameter< double >::type delta(deltaSEXP); - Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP); - Rcpp::traits::input_parameter< double >::type beta(betaSEXP); + Rcpp::traits::input_parameter< double >::type eta(etaSEXP); Rcpp::traits::input_parameter< int >::type seed(seedSEXP); Rcpp::traits::input_parameter< int >::type n_sweep(n_sweepSEXP); Rcpp::traits::input_parameter< int >::type burn(burnSEXP); Rcpp::traits::input_parameter< int >::type freeze_after(freeze_afterSEXP); Rcpp::traits::input_parameter< bool >::type slab_cauchy(slab_cauchySEXP); - rcpp_result_gen = Rcpp::wrap(zratio_test_calibrated_eval(graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, sigma, beta, seed, n_sweep, burn, freeze_after, slab_cauchy)); + rcpp_result_gen = Rcpp::wrap(zratio_test_calibrated_eval(graphs, edges, addc, tg, ihat, ghat, wt, psi0, delta, eta, seed, n_sweep, burn, freeze_after, slab_cauchy)); return rcpp_result_gen; END_RCPP } @@ -923,10 +921,10 @@ static const R_CallMethodDef CallEntries[] = { {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4}, {"_bgms_test_warmup_schedule", (DL_FUNC) &_bgms_test_warmup_schedule, 5}, {"_bgms_zratio_scan_graph", (DL_FUNC) &_bgms_zratio_scan_graph, 7}, - {"_bgms_zratio_audit_edges", (DL_FUNC) &_bgms_zratio_audit_edges, 15}, + {"_bgms_zratio_audit_edges", (DL_FUNC) &_bgms_zratio_audit_edges, 14}, {"_bgms_zratio_test_eval", (DL_FUNC) &_bgms_zratio_test_eval, 8}, {"_bgms_zratio_test_saddle", (DL_FUNC) &_bgms_zratio_test_saddle, 8}, - {"_bgms_zratio_test_calibrated_eval", (DL_FUNC) &_bgms_zratio_test_calibrated_eval, 16}, + {"_bgms_zratio_test_calibrated_eval", (DL_FUNC) &_bgms_zratio_test_calibrated_eval, 15}, {"_bgms_zratio_test_precompute", (DL_FUNC) &_bgms_zratio_test_precompute, 10}, {NULL, NULL, 0} }; diff --git a/src/models/ggm/zratio_engine.cpp b/src/models/ggm/zratio_engine.cpp index c270d02b..a7a7db7f 100644 --- a/src/models/ggm/zratio_engine.cpp +++ b/src/models/ggm/zratio_engine.cpp @@ -245,15 +245,15 @@ bool ZRatioEngine::audit_edge(const arma::imat& G, int i, int j, return true; } -void ZRatioEngine::enable_calibration(double delta, double sigma, double beta, - SafeRNG* rng, int n_sweep, int burn, +void ZRatioEngine::enable_calibration(double delta, double eta, SafeRNG* rng, + int n_sweep, int burn, double maha_thresh, int min_anchors, bool slab_cauchy) { calibration_enabled_ = true; frozen_ = false; delta_ = delta; - sigma_ = sigma; - beta_ = beta; + sigma_ = 1.0; + beta_ = eta; rng_ = rng; n_sweep_ = n_sweep; burn_ = burn; diff --git a/src/models/ggm/zratio_engine.h b/src/models/ggm/zratio_engine.h index 02f8cf8d..90886bf5 100644 --- a/src/models/ggm/zratio_engine.h +++ b/src/models/ggm/zratio_engine.h @@ -118,18 +118,19 @@ class ZRatioEngine { /** * Set the standardized-cell prior constants and RNG the block-Gibbs - * oracle samples under, without entering calibration mode. rng must - * outlive the engine. slab_cauchy selects the Cauchy slab family: the - * block couplings then run omega-augmented (scale-mixture of normals) - * and the endpoint legs mix per sweep, matching the marginal-Cauchy - * normalizer the tables integrate. + * oracle samples under, without entering calibration mode. The frame is + * standardized (unit slab scale), so only the diagonal rate eta is free; + * sigma is fixed to 1. rng must outlive the engine. slab_cauchy selects + * the Cauchy slab family: the block couplings then run omega-augmented + * (scale-mixture of normals) and the endpoint legs mix per sweep, + * matching the marginal-Cauchy normalizer the tables integrate. */ - void set_oracle_params(double delta, double sigma, double beta, - SafeRNG* rng, int n_sweep = 300, int burn = 30, + void set_oracle_params(double delta, double eta, SafeRNG* rng, + int n_sweep = 300, int burn = 30, bool slab_cauchy = false) { delta_ = delta; - sigma_ = sigma; - beta_ = beta; + sigma_ = 1.0; + beta_ = eta; rng_ = rng; n_sweep_ = n_sweep; burn_ = burn; @@ -147,13 +148,13 @@ class ZRatioEngine { * box into the addc layout, after which the engine behaves exactly * like one constructed with a full 23-slot constant block. * - * (delta, sigma, beta) are the standardized-cell prior constants the - * oracle samples under; rng must outlive the engine (the model's - * chain RNG). slab_cauchy selects the Cauchy slab family for the - * oracle (see set_oracle_params). + * (delta, eta) are the standardized-cell prior constants the oracle + * samples under (unit slab scale, diagonal rate eta); rng must outlive + * the engine (the model's chain RNG). slab_cauchy selects the Cauchy + * slab family for the oracle (see set_oracle_params). */ - void enable_calibration(double delta, double sigma, double beta, - SafeRNG* rng, int n_sweep = 300, int burn = 30, + void enable_calibration(double delta, double eta, SafeRNG* rng, + int n_sweep = 300, int burn = 30, double maha_thresh = 9.0, int min_anchors = 6, bool slab_cauchy = false); diff --git a/src/sample_ggm.cpp b/src/sample_ggm.cpp index 7b465c86..ad159fd9 100644 --- a/src/sample_ggm.cpp +++ b/src/sample_ggm.cpp @@ -129,8 +129,7 @@ Rcpp::List sample_ggm( // The rng pointer is rebound per chain clone by GGMModel. engine->enable_calibration( Rcpp::as(zs["delta"]), - Rcpp::as(zs["sigma"]), - Rcpp::as(zs["beta"]), nullptr, 300, 30, 9.0, 6, + Rcpp::as(zs["eta"]), nullptr, 300, 30, 9.0, 6, zr_cauchy); } model.set_zratio_engine(std::move(engine)); diff --git a/src/sample_mixed.cpp b/src/sample_mixed.cpp index 10c51d6a..ce93015e 100644 --- a/src/sample_mixed.cpp +++ b/src/sample_mixed.cpp @@ -187,8 +187,7 @@ Rcpp::List sample_mixed_mrf( // The rng pointer is rebound per chain clone by MixedMRFModel. engine->enable_calibration( Rcpp::as(zs["delta"]), - Rcpp::as(zs["sigma"]), - Rcpp::as(zs["beta"]), nullptr, 300, 30, 9.0, 6, + Rcpp::as(zs["eta"]), nullptr, 300, 30, 9.0, 6, zr_cauchy); } model.set_zratio_engine(std::move(engine)); diff --git a/src/zratio_audit_interface.cpp b/src/zratio_audit_interface.cpp index 8402b44c..2ac67deb 100644 --- a/src/zratio_audit_interface.cpp +++ b/src/zratio_audit_interface.cpp @@ -72,9 +72,9 @@ arma::mat zratio_scan_graph( // Measurement-only oracle audit of picked edges on one graph. Per edge: // the deployed correction under `addc` versus the block-Gibbs local // oracle discrepancy log(saddle on oracle moments) - log(additive), and -// the audit score |pred - oracle|. (delta, sigma, beta) are the -// standardized-cell prior constants; slab_cauchy selects the Cauchy -// slab family for the oracle. Rows with ok = 0 had a one-sided block, +// the audit score |pred - oracle|. (delta, eta) are the standardized-cell +// prior constants (unit slab scale); slab_cauchy selects the Cauchy slab +// family for the oracle. Rows with ok = 0 had a one-sided block, // non-positive additive moments, or an oracle with no finite sweep; // their scores are NA. // ----------------------------------------------------------------------------- @@ -90,8 +90,7 @@ Rcpp::List zratio_audit_edges( arma::vec wt, double psi0, double delta, - double sigma, - double beta, + double eta, int n_sweep, int burn, int seed, @@ -99,8 +98,7 @@ Rcpp::List zratio_audit_edges( ) { ZRatioEngine engine(addc, tg, ihat, ghat, wt, psi0); SafeRNG rng(seed); - engine.set_oracle_params(delta, sigma, beta, &rng, n_sweep, burn, - slab_cauchy); + engine.set_oracle_params(delta, eta, &rng, n_sweep, burn, slab_cauchy); const arma::uword n = edges.n_rows; arma::vec err(n), pred(n), oracle(n); arma::ivec ok(n); diff --git a/src/zratio_test_interface.cpp b/src/zratio_test_interface.cpp index 64208cfb..f61cdd7f 100644 --- a/src/zratio_test_interface.cpp +++ b/src/zratio_test_interface.cpp @@ -90,8 +90,7 @@ Rcpp::List zratio_test_calibrated_eval( arma::vec wt, double psi0, double delta, - double sigma, - double beta, + double eta, int seed, int n_sweep, int burn, @@ -100,7 +99,7 @@ Rcpp::List zratio_test_calibrated_eval( ) { ZRatioEngine engine(addc, tg, ihat, ghat, wt, psi0); SafeRNG rng(seed); - engine.enable_calibration(delta, sigma, beta, &rng, n_sweep, burn, 9.0, 6, + engine.enable_calibration(delta, eta, &rng, n_sweep, burn, 9.0, 6, slab_cauchy); arma::vec out(edges.n_rows); for (arma::uword e = 0; e < edges.n_rows; ++e) { diff --git a/tests/testthat/test-hier-zratio-identity.R b/tests/testthat/test-hier-zratio-identity.R index 78008f3f..e0f92f07 100644 --- a/tests/testthat/test-hier-zratio-identity.R +++ b/tests/testthat/test-hier-zratio-identity.R @@ -69,8 +69,7 @@ test_that("Z-ratio constants build in the standardized cell", { ) expect_identical(a, b) expect_identical(a, e) - expect_identical(a$sigma, 1) - expect_identical(a$beta, 1) + expect_identical(a$eta, 1) }) test_that("hierarchical graph marginal holds at a non-unit slab scale", { diff --git a/tests/testthat/test-zratio-alarms.R b/tests/testthat/test-zratio-alarms.R index ee2bea7a..2723182f 100644 --- a/tests/testthat/test-zratio-alarms.R +++ b/tests/testthat/test-zratio-alarms.R @@ -70,7 +70,7 @@ test_that("the audit flags a defective frozen kernel", { skip_on_cran() p = 10 delta = 0.5 * log(p) - zc = bgms:::zratio_constants(delta, sigma = 1, beta = 1) + zc = bgms:::zratio_constants(delta, eta = 1) # Pack a fit whose correction is a constant 0.05 on every coupled-bridge # block, with a hull box wide enough that the clamp never engages. The # oracle discrepancy on healthy blocks is O(0.001), so every targeted @@ -114,7 +114,7 @@ test_that("the audit flags a defective frozen kernel", { ) spec = list( tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, wt = zc$wt, psi0 = zc$psi0, - delta = delta, sigma = 1, beta = 1 + delta = delta, eta = 1 ) old = options(bgms.verbose = TRUE) on.exit(options(old), add = TRUE) diff --git a/tests/testthat/test-zratio-cauchy.R b/tests/testthat/test-zratio-cauchy.R index d28e2911..2c36b95b 100644 --- a/tests/testthat/test-zratio-cauchy.R +++ b/tests/testthat/test-zratio-cauchy.R @@ -107,7 +107,7 @@ test_that("Cauchy block-Gibbs oracle tracks the additive prediction", { audit = zratio_audit_edges( G, matrix(c(1L, 2L), 1, 2), zc$addc, zc$tg, zc$ihat, zc$ghat, zc$wt, zc$psi0, - zc$delta, zc$sigma, zc$beta, 2000L, 50L, 42L, + zc$delta, zc$eta, 2000L, 50L, 42L, identical(slab, "cauchy") ) expect_identical(as.integer(audit$ok), 1L, label = slab) diff --git a/tests/testthat/test-zratio-engine.R b/tests/testthat/test-zratio-engine.R index d9f9c1d5..248bb612 100644 --- a/tests/testthat/test-zratio-engine.R +++ b/tests/testthat/test-zratio-engine.R @@ -40,8 +40,10 @@ test_that("engine reproduces the reference log Z-ratios (all cells/variants)", { test_that("fit-time constant builders match the reference builders", { skip_if_not(file.exists(fixture_path), "zratio_reference.rds not generated") fx = readRDS(fixture_path) - for(cell in fx[c(1, 3)]) { # two cells keep the runtime modest - zc = bgms:::zratio_constants(cell$delta, cell$sigma, cell$beta) + # The builders run in the standardized cell (sigma = 1), so match against + # the two sigma = 1 reference cells; eta = sigma * beta = beta there. + for(cell in fx[c(1, 4)]) { + zc = bgms:::zratio_constants(cell$delta, cell$sigma * cell$beta) expect_equal(zc$addc, cell$addc6, tolerance = 1e-8) expect_equal(zc$psi0, cell$psi0, tolerance = 1e-8) expect_equal(zc$tg, cell$tg, tolerance = 1e-12) @@ -120,9 +122,8 @@ test_that("online calibrator: oracle taper, freeze pack, accuracy gain", { skip_if_not(file.exists(fixture_path), "zratio_reference.rds not generated") q = 14L dlt = 0.5 * log(30) - sg = 2 - bt = 0.5 - zc = bgms:::zratio_constants(dlt, sg, bt) + et = 1 + zc = bgms:::zratio_constants(dlt, et) set.seed(4) graphs = list() @@ -156,7 +157,7 @@ test_that("online calibrator: oracle taper, freeze pack, accuracy gain", { truth = vapply(seq_len(n), function(e) { bgms:::zratio_test_calibrated_eval( graphs[e], edges[e, , drop = FALSE], zc$addc, zc$tg, zc$ihat, - zc$ghat, zc$wt, zc$psi0, dlt, sg, bt, + zc$ghat, zc$wt, zc$psi0, dlt, et, seed = 500L + e, n_sweep = 2000L, burn = 50L, freeze_after = 0L )$log_zratio[1] }, 0.0) @@ -164,7 +165,7 @@ test_that("online calibrator: oracle taper, freeze pack, accuracy gain", { fz = ceiling(0.6 * n) res = bgms:::zratio_test_calibrated_eval( graphs, edges, zc$addc, zc$tg, zc$ihat, zc$ghat, zc$wt, zc$psi0, - dlt, sg, bt, + dlt, et, seed = 7L, n_sweep = 300L, burn = 30L, freeze_after = fz ) add = vapply(seq_len(n), function(e) { From f122e56fbea1e6b10c6a666aedcdff3e3d344c5e Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Tue, 7 Jul 2026 08:58:38 +0200 Subject: [PATCH 4/4] Add the dense high-q Cauchy graph-law identity gate Slow-gated standing check that the hierarchical Cauchy spec reproduces the edge prior in the coupling regime (q = 15, 20; p_inc = 0.5, 0.7; both update methods), where the single-edge memo flags the mixture slab as coupling-sensitive. The deployed system (additive saddle + warm-up OLS correction, default window engaging at p >= 15) holds the marginal to max |dev| = 0.0033. Asserts the graph-law marginal; the per-block alarm verdict is expected to flag at dense high q (it does for the Normal slab too) and is not asserted. --- tests/testthat/test-hier-zratio-identity.R | 34 ++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/tests/testthat/test-hier-zratio-identity.R b/tests/testthat/test-hier-zratio-identity.R index e0f92f07..d4b747a1 100644 --- a/tests/testthat/test-hier-zratio-identity.R +++ b/tests/testthat/test-hier-zratio-identity.R @@ -104,6 +104,40 @@ test_that("hierarchical graph marginal holds for the Cauchy slab", { } }) +test_that("Cauchy graph law holds at dense high q (coupling regime)", { + skip_on_cran() + skip_if( + !identical(Sys.getenv("BGMS_RUN_SLOW_TESTS"), "true"), + "Set BGMS_RUN_SLOW_TESTS=true to run the dense high-q Cauchy identity" + ) + # The single-edge memo flags the Gaussian-mixture slab as a coupling- + # sensitive case (the additive moment approximation weakens on dense + # blocks). At dense high q (q = 15, 20; p_inc = 0.5, 0.7) the deployed + # system -- additive saddle + warm-up OLS correction on maxbd >= 2 blocks + # (the default window engages at p >= 15) -- must still reproduce the edge + # prior. The marginal is the release-relevant statistic; the per-block + # alarm verdict is expected to flag here (it does for the Normal slab too) + # and is not asserted. + suppressMessages(library(parallel)) + cells = expand.grid( + q = c(15L, 20L), p_inc = c(0.5, 0.7), + um = c("adaptive-metropolis", "gibbs"), stringsAsFactors = FALSE + ) + devs = unlist(mclapply(seq_len(nrow(cells)), function(r) { + cl = cells[r, ] + d = sample_ggm_prior( + p = cl$q, n_samples = 6000L, n_warmup = 2000L, + interaction_prior = cauchy_prior(scale = 0.5), + precision_scale_prior = gamma_prior(shape = 1, rate = 2), + spec = "hierarchical", edge_inclusion_prob = cl$p_inc, + update_method = cl$um, delta = 0.5 * log(cl$q), + seed = 4000L + r, verbose = FALSE + ) + mean(d$edge_indicators) - cl$p_inc + }, mc.cores = min(8L, nrow(cells)))) + expect_lt(max(abs(devs)), 0.02) +}) + test_that("hierarchical BB identity: theta ~ Beta(a, b), PIP = a/(a+b)", { skip_on_cran() d = hier_prior_run(