Accept a Cauchy slab under the hierarchical spec; rename graph_prior_spec

NEWS.md

@@ -27,7 +27,7 @@
 * 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(graph_prior_spec = "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()` 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 NUTS diagnostics summary now prints the `warmup_incomplete` flag (energy not stationary) it already computed.
 * `extract_prior_inclusion_probabilities()`: prior edge-inclusion probabilities in the same matrix layout as `extract_posterior_inclusion_probabilities()`, for prior/posterior inclusion-odds computations. Under the joint spike-and-slab prior on a continuous block the graph marginal is reweighted by the per-graph normalizer (positive-definite-cone mass shaped by the determinant tilt), so continuous-continuous edges do not keep the edge-prior marginal at any `delta` — with a uniform Beta-Bernoulli hyperprior at 3 variables the prior edge probability is about 0.37 rather than 0.5, and a fixed `bernoulli_prior(0.5)` at `delta = 0` lands near 0.27. The values are read from the cached correction table (`bernoulli_prior()`, `beta_bernoulli_prior()`) or estimated by a prior-only chain with the fit's own correction settings (`sbm_prior()`, cached on the fit). Mixed models report per-edge-class values (discrete-discrete, continuous-continuous, cross); ordinal models use the analytic edge-prior marginals, including the exchangeable partition mixture for `sbm_prior()`.

R/bgm.R

@@ -168,7 +168,7 @@
 #'   sampled inclusion probability is returned per chain in
 #'   \code{fit$inclusion_parameter_samples}.
 #'
-#' @param graph_prior_spec Character. How the precision prior composes with
+#' @param precision_graph_prior Character. How the precision prior composes with
 #'   the edge prior under edge selection for continuous (GGM) data:
 #'   \describe{
 #'     \item{"joint"}{(default) The un-normalised joint specification
@@ -189,10 +189,10 @@
 #'       sampling (\code{\link{summarize_zratio_diagnostics}}; the summary
 #'       is returned as \code{fit$zratio_diag} and issues print like other
 #'       sampler warnings). Requires \code{edge_selection = TRUE}, a
-#'       \code{normal_prior()} interaction prior, a shape-1
-#'       \code{gamma_prior()} (or \code{exponential_prior()}) precision
-#'       scale prior, and continuous data — either all-continuous (GGM) or
-#'       mixed with at least two continuous variables. On mixed data the
+#'       \code{normal_prior()} or \code{cauchy_prior()} interaction prior, a
+#'       shape-1 \code{gamma_prior()} (or \code{exponential_prior()})
+#'       precision scale prior, and continuous data — either all-continuous
+#'       (GGM) or mixed with at least two continuous variables. On mixed data the
 #'       normalizer lives on the continuous block \eqn{K_{yy}}, so the
 #'       Z-ratio enters the continuous-continuous edge moves only, with the
 #'       mediating-block counts read off the continuous subgraph; discrete
@@ -201,7 +201,7 @@
 #'   Default: \code{"joint"}.
 #'
 #' @param calibration_window Non-negative integer or \code{NULL} (default).
-#'   Only for \code{graph_prior_spec = "hierarchical"}: length of the
+#'   Only for \code{precision_graph_prior = "hierarchical"}: length of the
 #'   appended warm-up window in which the Z-ratio correction is calibrated
 #'   online against a block-Gibbs oracle and then frozen (with its hull
 #'   clamp) before sampling. \code{NULL} resolves to no window for
@@ -407,7 +407,7 @@ bgm = function(
   delta = NULL,
   edge_selection = TRUE,
   edge_prior = bernoulli_prior(0.5),
-  graph_prior_spec = c("joint", "hierarchical"),
+  precision_graph_prior = c("joint", "hierarchical"),
   calibration_window = NULL,
   na_action = c("listwise", "impute"),
   update_method = c("nuts", "adaptive-metropolis", "gibbs"),
@@ -577,7 +577,7 @@ bgm = function(
     delta = delta,
     edge_selection = edge_selection,
     edge_prior = edge_prior,
-    graph_prior_spec = graph_prior_spec,
+    precision_graph_prior = precision_graph_prior,
     calibration_window = calibration_window,
     update_method = update_method,
     target_accept = if(hasArg(target_accept)) target_accept else NULL,

R/bgm_spec.R

@@ -268,7 +268,7 @@ bgm_spec = function(x,
                     delta = NULL,
                     edge_selection = TRUE,
                     edge_prior = bernoulli_prior(0.5),
-                    graph_prior_spec = c("joint", "hierarchical"),
+                    precision_graph_prior = c("joint", "hierarchical"),
                     calibration_window = NULL,
                     # Legacy edge prior params (accepted for backward compat)
                     inclusion_probability = 0.5,
@@ -379,46 +379,52 @@ bgm_spec = function(x,
   # The Z-ratio constants are derived for the Normal slab with an exponential
   # (shape-1 Gamma) diagonal, on the continuous precision matrix, under edge
   # selection. Anything else keeps the joint specification.
-  graph_prior_spec = match.arg(graph_prior_spec)
-  if(graph_prior_spec == "hierarchical") {
+  precision_graph_prior = match.arg(precision_graph_prior)
+  if(precision_graph_prior == "hierarchical") {
     if(!model_type %in% c("ggm", "mixed_mrf")) {
       stop(
-        "graph_prior_spec = \"hierarchical\" needs a continuous precision ",
+        "precision_graph_prior = \"hierarchical\" needs a continuous precision ",
         "block to normalize; the current model_type is '", model_type,
         "'. Use the joint specification, or data with continuous variables."
       )
     }
     if(model_type == "mixed_mrf" && sum(variable_type == "continuous") < 2) {
       stop(
-        "graph_prior_spec = \"hierarchical\" on mixed data needs at least ",
+        "precision_graph_prior = \"hierarchical\" on mixed data needs at least ",
         "two continuous variables (the specification normalizes the ",
         "continuous-block prior across its graphs). Use the joint ",
         "specification."
       )
     }
     if(!edge_selection) {
       stop(
-        "graph_prior_spec = \"hierarchical\" normalizes p(K | Gamma) ",
+        "precision_graph_prior = \"hierarchical\" normalizes p(K | Gamma) ",
         "across graphs and needs edge_selection = TRUE; with a fixed ",
         "graph the specifications coincide."
       )
     }
-    if(!identical(interaction_prior_type, "normal")) {
-      stop(
-        "graph_prior_spec = \"hierarchical\" requires a normal ",
-        "interaction (slab) prior; the Z-ratio constants are derived for ",
-        "the Normal slab. Use interaction_prior = normal_prior(), or the ",
-        "joint specification."
-      )
+    if(!interaction_prior_type %in% c("normal", "cauchy")) {
+      stop(sprintf(
+        paste0(
+          "precision_graph_prior = \"hierarchical\" supports a normal or Cauchy ",
+          "interaction (slab) prior. Got %s_prior(). Use interaction_prior = ",
+          "normal_prior() or cauchy_prior(), or keep precision_graph_prior = ",
+          "\"joint\"."
+        ),
+        interaction_prior_type
+      ))
     }
     if(abs(scale_shape - 1) > 1e-12) {
-      stop(
-        "graph_prior_spec = \"hierarchical\" requires shape = 1 on the ",
-        "precision scale prior (gamma_prior(shape = 1) or ",
-        "exponential_prior()); the Z-ratio constants are derived for the ",
-        "exponential diagonal. Adjust the prior, or use the joint ",
-        "specification."
-      )
+      stop(sprintf(
+        paste0(
+          "precision_graph_prior = \"hierarchical\" requires shape = 1 on the ",
+          "precision scale prior; the Z-ratio normalizer is derived for ",
+          "the exponential diagonal. Got shape = %s. Use ",
+          "gamma_prior(shape = 1) or exponential_prior(), or keep ",
+          "precision_graph_prior = \"joint\"."
+        ),
+        format(scale_shape)
+      ))
     }
   }
 
@@ -495,7 +501,7 @@ bgm_spec = function(x,
       scale_rate = scale_rate,
       scale_eta = scale_eta,
       delta = delta,
-      graph_prior_spec = graph_prior_spec,
+      precision_graph_prior = precision_graph_prior,
       calibration_window = calibration_window,
       edge_prior_flat = ep_flat
     )
@@ -522,7 +528,7 @@ bgm_spec = function(x,
       scale_rate = scale_rate,
       scale_eta = scale_eta,
       delta = delta,
-      graph_prior_spec = graph_prior_spec,
+      precision_graph_prior = precision_graph_prior,
       calibration_window = calibration_window,
       edge_prior_flat = ep_flat
     )