feat(ggm): log Z(G) closed-form primitives for hierarchical-spec MH

.Rbuildignore

@@ -33,8 +33,10 @@
 ^paper$
 ^LICENSE$
 
-# ---- Legacy test fixtures (GitHub CI only, not shipped to CRAN) ----
-^tests/testthat/fixtures$
+# ---- Legacy test fixtures: only the `legacy` subdirectory is excluded
+#      so it does not bloat the tarball; live fixtures used by tests in
+#      R CMD check must ship so the matching tests can run.
+^tests/testthat/fixtures/legacy$
 
 # ---- C/C++ build artifacts (REQUIRED) ----
 ^src/.*\.o$

NAMESPACE

@@ -49,6 +49,7 @@ export(beta_bernoulli_prior)
 export(beta_prime_prior)
 export(bgm)
 export(bgmCompare)
+export(bgms_posterior_mean)
 export(cauchy_prior)
 export(exponential_prior)
 export(extract_arguments)

R/RcppExports.R

@@ -53,6 +53,70 @@ sample_ggm_prior_cpp <- function(p, n_samples, n_warmup = 1000L, pairwise_scale
     .Call(`_bgms_sample_ggm_prior`, p, n_samples, n_warmup, pairwise_scale, interaction_prior_type, scale_prior_type, gamma_shape, gamma_rate, step_size, max_depth, seed, verbose, edge_indicators_nullable, delta)
 }
 
+log_Z_NLO_gamma_cpp <- function(G, alpha, beta, sigma, include_F = FALSE, delta = 0.0) {
+    .Call(`_bgms_log_Z_NLO_gamma_cpp`, G, alpha, beta, sigma, include_F, delta)
+}
+
+log_Z_NLO_gamma_degord_cpp <- function(G, i, j, alpha, beta, sigma, include_F = FALSE, delta = 0.0) {
+    .Call(`_bgms_log_Z_NLO_gamma_degord_cpp`, G, i, j, alpha, beta, sigma, include_F, delta)
+}
+
+log_Z_NLO_gamma_delta_incr_alpha1_cpp <- function(G_before, i, j, beta, sigma, delta, include_F = FALSE) {
+    .Call(`_bgms_log_Z_NLO_gamma_delta_incr_alpha1_cpp`, G_before, i, j, beta, sigma, delta, include_F)
+}
+
+log_Z_NLO_gamma_delta_incr_alphaN_cpp <- function(G_before, i, j, alpha, beta, sigma, delta, include_F = FALSE) {
+    .Call(`_bgms_log_Z_NLO_gamma_delta_incr_alphaN_cpp`, G_before, i, j, alpha, beta, sigma, delta, include_F)
+}
+
+degord_chain_aux_cpp <- function(q, alpha, beta, sigma, delta) {
+    .Call(`_bgms_degord_chain_aux_cpp`, q, alpha, beta, sigma, delta)
+}
+
+degord_pi_aux_cpp <- function(G_pi, alpha, beta, sigma, delta) {
+    .Call(`_bgms_degord_pi_aux_cpp`, G_pi, alpha, beta, sigma, delta)
+}
+
+degord_permute_graph_cpp <- function(G, i, j) {
+    .Call(`_bgms_degord_permute_graph_cpp`, G, i, j)
+}
+
+degord_log_Zhat_pi_from_pool_cpp <- function(noise_pool_t, G_pi, alpha, beta, sigma, delta, slab_tilt_mode = 0L) {
+    .Call(`_bgms_degord_log_Zhat_pi_from_pool_cpp`, noise_pool_t, G_pi, alpha, beta, sigma, delta, slab_tilt_mode)
+}
+
+degord_delta_log_Zhat_pi_toggle_cpp <- function(noise_pool, noise_pool_t, G_curr, i, j, alpha, beta, sigma, delta, slab_tilt_mode = 0L) {
+    .Call(`_bgms_degord_delta_log_Zhat_pi_toggle_cpp`, noise_pool, noise_pool_t, G_curr, i, j, alpha, beta, sigma, delta, slab_tilt_mode)
+}
+
+degord_draw_bartlett_pool_cpp <- function(q, M_inner, seed) {
+    .Call(`_bgms_degord_draw_bartlett_pool_cpp`, q, M_inner, seed)
+}
+
+degord_V_at_Gamma_pi_cpp <- function(K_depth, pools_t, G_pi, alpha, beta, sigma, delta, c_val, rho, slab_tilt_mode = 0L) {
+    .Call(`_bgms_degord_V_at_Gamma_pi_cpp`, K_depth, pools_t, G_pi, alpha, beta, sigma, delta, c_val, rho, slab_tilt_mode)
+}
+
+degord_V_log_at_Gamma_pi_cpp <- function(K_depth, pools_t, G_pi, alpha, beta, sigma, delta, log_c, rho, slab_tilt_mode = 0L) {
+    .Call(`_bgms_degord_V_log_at_Gamma_pi_cpp`, K_depth, pools_t, G_pi, alpha, beta, sigma, delta, log_c, rho, slab_tilt_mode)
+}
+
+degord_V_log_pair_at_Gamma_curr_star_cpp <- function(K_depth, pools_t, G_pi_curr, G_pi_star, alpha, beta, sigma, delta, log_c_curr, log_c_star, rho, slab_tilt_mode = 0L) {
+    .Call(`_bgms_degord_V_log_pair_at_Gamma_curr_star_cpp`, K_depth, pools_t, G_pi_curr, G_pi_star, alpha, beta, sigma, delta, log_c_curr, log_c_star, rho, slab_tilt_mode)
+}
+
+degord_log_Zhat_star_from_cache_cpp <- function(noise_pool_t, G_pi_curr, G_pi_star, alpha, beta, sigma, delta, slab_tilt_mode = 0L) {
+    .Call(`_bgms_degord_log_Zhat_star_from_cache_cpp`, noise_pool_t, G_pi_curr, G_pi_star, alpha, beta, sigma, delta, slab_tilt_mode)
+}
+
+degord_draw_U_rr_cpp <- function(M_inner, q, rho, seed) {
+    .Call(`_bgms_degord_draw_U_rr_cpp`, M_inner, q, rho, seed)
+}
+
+ggm_hierarchical_smoke_cpp <- function(observations, inclusion_prob, interaction_scale, diagonal_shape, diagonal_rate, delta, M_inner, kappa, rho, n_sweeps, seed) {
+    .Call(`_bgms_ggm_hierarchical_smoke_cpp`, observations, inclusion_prob, interaction_scale, diagonal_shape, diagonal_rate, delta, M_inner, kappa, rho, n_sweeps, seed)
+}
+
 .compute_ess_cpp <- function(array3d) {
     .Call(`_bgms_compute_ess_cpp`, array3d)
 }
@@ -141,8 +205,8 @@ ggm_test_logp_and_gradient_full_prior <- function(x, suf_stat, n, edge_indicator
     .Call(`_bgms_ggm_test_logp_and_gradient_full_prior`, x, suf_stat, n, edge_indicators, interaction_prior_type, interaction_scale, interaction_alpha, interaction_beta, diagonal_prior_type, diagonal_shape, diagonal_rate, inv_mass_diag)
 }
 
-sample_ggm <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, target_acceptance = 0.8, max_tree_depth = 10L, na_impute = FALSE, missing_index_nullable = NULL, delta = 0.0) {
-    .Call(`_bgms_sample_ggm`, inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, target_acceptance, max_tree_depth, na_impute, missing_index_nullable, delta)
+sample_ggm <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, target_acceptance = 0.8, max_tree_depth = 10L, na_impute = FALSE, missing_index_nullable = NULL, delta = 0.0, graph_prior_spec = "joint", z_ratio_M_inner = 100L, z_ratio_kappa = 1.0, z_ratio_rho = 0.5) {
+    .Call(`_bgms_sample_ggm`, inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, target_acceptance, max_tree_depth, na_impute, missing_index_nullable, delta, graph_prior_spec, z_ratio_M_inner, z_ratio_kappa, z_ratio_rho)
 }
 
 sample_mixed_mrf <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, sampler_type = "mh", target_acceptance = 0.80, max_tree_depth = 10L, na_impute = FALSE, missing_index_discrete_nullable = NULL, missing_index_continuous_nullable = NULL, delta = 0.0) {

R/bgm.R

@@ -111,6 +111,30 @@
 #'   apply the tilt. Not allowed for pure ordinal models (no precision
 #'   matrix to tilt).
 #'
+#' @param graph_prior_spec Character; one of \code{"joint"} (default)
+#'   or \code{"hierarchical"}. Controls the marginal prior on the graph
+#'   indicators \eqn{\Gamma}. Under \code{"joint"} the implicit
+#'   \eqn{\Gamma}-marginal is \eqn{\pi(\Gamma) \cdot Z(\Gamma)}, where
+#'   \eqn{Z(\Gamma)} is the normalising constant of the precision-matrix
+#'   prior conditional on the graph. Under \code{"hierarchical"} the
+#'   chain compensates with an unbiased estimator of
+#'   \eqn{Z(\Gamma_\text{curr}) / Z(\Gamma_\text{star})}, recovering
+#'   \eqn{\pi(\Gamma)} as the \eqn{\Gamma}-marginal. Only supported when
+#'   the interaction prior is \code{normal_prior(...)} and the precision-
+#'   scale prior is \code{gamma_prior(...)} (the prior families for which
+#'   the closed-form Laplace-NLO normaliser approximation is implemented).
+#'   Default: \code{"joint"}.
+#'
+#' @param z_ratio_tuning Named list with components \code{M_inner}
+#'   (positive integer, default 100), \code{kappa} (positive numeric,
+#'   default 1.0), and \code{rho} (numeric in (0, 1), default 0.5).
+#'   Tuning knobs for the V/Russian-Roulette estimator used when
+#'   \code{graph_prior_spec = "hierarchical"}; ignored otherwise.
+#'   \code{M_inner} is the number of inner Bartlett-Cholesky importance
+#'   samples per Russian-Roulette pool, \code{kappa} sets the analytic
+#'   centring \eqn{c = \kappa \exp(\log Z_\text{NLO}(\Gamma))}, and
+#'   \code{rho} is the geometric-truncation continuation probability.
+#'
 #' @param pairwise_scale `r lifecycle::badge("deprecated")` Double.
 #'   Scale of the Cauchy prior for pairwise
 #'   interaction parameters. Use \code{interaction_prior} instead.
@@ -343,6 +367,8 @@ bgm = function(
   means_prior = normal_prior(scale = 1),
   precision_scale_prior = gamma_prior(shape = 1, rate = 1),
   delta = NULL,
+  graph_prior_spec = c("joint", "hierarchical"),
+  z_ratio_tuning = list(M_inner = 100L, kappa = 1.0, rho = 0.5),
   edge_selection = TRUE,
   edge_prior = bernoulli_prior(0.5),
   na_action = c("listwise", "impute"),
@@ -511,6 +537,8 @@ bgm = function(
     scale_shape = sp$scale_shape,
     scale_rate = sp$scale_rate,
     delta = delta,
+    graph_prior_spec = graph_prior_spec,
+    z_ratio_tuning = z_ratio_tuning,
     standardize = standardize,
     edge_selection = edge_selection,
     edge_prior = edge_prior,

R/bgm_spec.R

@@ -266,6 +266,10 @@ bgm_spec = function(x,
                     scale_shape = 1,
                     scale_rate = 1,
                     delta = NULL,
+                    graph_prior_spec = c("joint", "hierarchical"),
+                    z_ratio_tuning = list(M_inner = 100L,
+                                          kappa   = 1.0,
+                                          rho     = 0.5),
                     standardize = FALSE,
                     edge_selection = TRUE,
                     edge_prior = bernoulli_prior(0.5),
@@ -364,6 +368,61 @@ bgm_spec = function(x,
      !is.finite(delta) || delta < 0) {
     stop("'delta' must be a single finite non-negative numeric, or NULL.")
   }
+  # Validate hierarchical-spec args (only meaningful for ggm/mixed_mrf).
+  graph_prior_spec = if(is.character(graph_prior_spec) &&
+                        length(graph_prior_spec) > 1L) {
+    match.arg(graph_prior_spec)
+  } else {
+    if(!(length(graph_prior_spec) == 1L &&
+         is.character(graph_prior_spec) &&
+         graph_prior_spec %in% c("joint", "hierarchical"))) {
+      stop("'graph_prior_spec' must be \"joint\" or \"hierarchical\".")
+    }
+    graph_prior_spec
+  }
+  if(graph_prior_spec == "hierarchical" &&
+     !(model_type %in% c("ggm", "mixed_mrf"))) {
+    stop(
+      "'graph_prior_spec = \"hierarchical\"' requires continuous data; ",
+      "the current model_type is '", model_type, "', which has no ",
+      "continuous precision block. Use \"joint\" or supply continuous data."
+    )
+  }
+  if(graph_prior_spec == "hierarchical" &&
+     interaction_prior_type != "normal") {
+    stop(
+      "'graph_prior_spec = \"hierarchical\"' requires a Normal slab ",
+      "prior (interaction_prior_type = \"normal\"). Re-fit with ",
+      "interaction_prior = normal_prior(scale = ...)."
+    )
+  }
+  if(graph_prior_spec == "hierarchical" &&
+     scale_prior_type != "gamma") {
+    stop(
+      "'graph_prior_spec = \"hierarchical\"' requires a Gamma diagonal ",
+      "prior (scale_prior_type = \"gamma\")."
+    )
+  }
+  # Validate z_ratio_tuning shape (only enforced if hierarchical; for joint
+  # the defaults pass through unused).
+  if(!is.list(z_ratio_tuning))
+    stop("'z_ratio_tuning' must be a list with components M_inner, kappa, rho.")
+  zrt_M_inner = z_ratio_tuning$M_inner %||% 100L
+  zrt_kappa   = z_ratio_tuning$kappa   %||% 1.0
+  zrt_rho     = z_ratio_tuning$rho     %||% 0.5
+  if(!is.numeric(zrt_M_inner) || length(zrt_M_inner) != 1L ||
+     !is.finite(zrt_M_inner) || zrt_M_inner < 1L)
+    stop("'z_ratio_tuning$M_inner' must be a positive integer.")
+  if(!is.numeric(zrt_kappa) || length(zrt_kappa) != 1L ||
+     !is.finite(zrt_kappa) || zrt_kappa <= 0)
+    stop("'z_ratio_tuning$kappa' must be a positive number.")
+  if(!is.numeric(zrt_rho) || length(zrt_rho) != 1L ||
+     !is.finite(zrt_rho) || zrt_rho <= 0 || zrt_rho >= 1)
+    stop("'z_ratio_tuning$rho' must be in (0, 1).")
+  z_ratio_tuning = list(M_inner = as.integer(zrt_M_inner),
+                        kappa   = as.numeric(zrt_kappa),
+                        rho     = as.numeric(zrt_rho))
+
   if(delta > 0 && model_type %in% c("omrf", "compare")) {
     stop(
       "'delta' (determinant tilt) requires continuous variables; the ",
@@ -444,6 +503,8 @@ bgm_spec = function(x,
       scale_shape = scale_shape,
       scale_rate = scale_rate,
       delta = delta,
+      graph_prior_spec = graph_prior_spec,
+      z_ratio_tuning = z_ratio_tuning,
       edge_prior_flat = ep_flat
     )
   } else if(model_type == "mixed_mrf") {
@@ -536,6 +597,10 @@ build_spec_ggm = function(x, data_columnnames, num_variables,
                           interaction_alpha, interaction_beta,
                           scale_prior_type, scale_shape, scale_rate,
                           delta = 0,
+                          graph_prior_spec = "joint",
+                          z_ratio_tuning = list(M_inner = 100L,
+                                                kappa   = 1.0,
+                                                rho     = 0.5),
                           edge_prior_flat) {
   # Missing data
   md = validate_missing_data(
@@ -577,6 +642,8 @@ build_spec_ggm = function(x, data_columnnames, num_variables,
       scale_shape = scale_shape,
       scale_rate = scale_rate,
       delta = delta,
+      graph_prior_spec = graph_prior_spec,
+      z_ratio_tuning = z_ratio_tuning,
       edge_selection = ep$edge_selection,
       edge_prior = ep$edge_prior,
       inclusion_probability = ep$inclusion_probability,

R/bgms_posterior_mean.R

@@ -0,0 +1,125 @@
+#' Sign-corrected posterior means for a bgms fit
+#'
+#' Computes posterior means using Lyne (2015)'s sign-corrected ergodic
+#' estimator. For a functional \eqn{f(K)} of the precision matrix,
+#'
+#' \deqn{E[f(K) \mid Y] \;\approx\; \frac{\sum_i \mathrm{sign}_i \cdot f(K_i)}{\sum_i \mathrm{sign}_i}}
+#'
+#' where \eqn{\mathrm{sign}_i} is the sign of the Russian-Roulette
+#' estimator of \eqn{1/Z(\Gamma)} at iteration \eqn{i}. The correction
+#' is required for the ergodic averages to converge to the true
+#' posterior expectation when running the hierarchical-spec GGM
+#' (\code{graph_prior_spec = "hierarchical"}) under settings that
+#' permit signed V values.
+#'
+#' Under operational tunings (low \eqn{\delta}, reasonably large
+#' \eqn{\kappa}) sign is identically \eqn{+1} and the correction
+#' collapses to the plain posterior mean. Sign correction matters
+#' primarily at high \eqn{\delta} or aggressive \eqn{\kappa}. The
+#' diagnostic vector lives at \code{fit$v_ratio_diagnostics$sign}.
+#'
+#' For joint-spec fits, or for any fit that has no
+#' \code{v_ratio_diagnostics} field, this function returns the
+#' \code{posterior_mean_*} fields unchanged.
+#'
+#' @param fit  A bgm() fit object.
+#' @return A list with components:
+#'   \describe{
+#'     \item{\code{main}}{Sign-corrected posterior mean of main effects.
+#'       \code{NULL} for GGM (which has no main effects).}
+#'     \item{\code{pairwise}}{Sign-corrected posterior mean of pairwise
+#'       associations, as a symmetric \eqn{p \times p} matrix. For GGM
+#'       this is on the association scale \eqn{-K_{ij}/2}.}
+#'     \item{\code{indicator}}{Sign-corrected posterior mean of edge
+#'       indicators \eqn{\gamma_{ij}}, as a symmetric \eqn{p \times p}
+#'       matrix. Present only when the fit used edge selection.}
+#'     \item{\code{residual_variance}}{Sign-corrected posterior mean
+#'       of the residual variance \eqn{1/K_{ii}}. Present only for GGM
+#'       fits.}
+#'   }
+#' @export
+bgms_posterior_mean = function(fit) {
+  diag = fit$v_ratio_diagnostics
+
+  # No sign data → identity fall-through to existing posterior means.
+  # (Joint-spec fits, or hierarchical fits before F2 wiring landed.)
+  if(is.null(diag)) {
+    return(list(
+      main              = fit$posterior_mean_main,
+      pairwise          = fit$posterior_mean_pairwise,
+      indicator         = fit$posterior_mean_indicator,
+      residual_variance = fit$posterior_mean_residual_variance
+    ))
+  }
+
+  raw = fit$raw_samples
+  signs = unlist(diag$sign)
+  sum_signs = sum(signs)
+  if(sum_signs == 0) {
+    stop(
+      "All sign(V) entries sum to zero; the sign-corrected posterior ",
+      "mean is undefined. This usually indicates that the V/RR estimator ",
+      "is operating outside its convergence regime - try refitting with ",
+      "a larger kappa (passed via z_ratio_tuning).",
+      call. = FALSE
+    )
+  }
+
+  # Pool per-chain samples and compute sign-weighted column means in one shot.
+  # Chain order in `signs` must match chain order in `raw$*` (both unlisted
+  # / rbind'd over the same lapply order).
+  weighted_col_means = function(samples_per_chain) {
+    if(is.null(samples_per_chain)) return(NULL)
+    pooled = do.call(rbind, samples_per_chain)
+    colSums(pooled * signs) / sum_signs
+  }
+
+  # Whether this is a GGM fit (the only model class that produces sign data
+  # currently). Heuristic: GGM has residual_variance, OMRF/mixed do not.
+  is_continuous = !is.null(fit$posterior_mean_residual_variance)
+  num_variables = nrow(fit$posterior_mean_pairwise)
+
+  out = list()
+
+  # Pairwise → symmetric matrix in association scale (GGM: precision * -0.5).
+  pair_means = weighted_col_means(raw$pairwise)
+  pairwise_mat = matrix(0, num_variables, num_variables,
+                        dimnames = dimnames(fit$posterior_mean_pairwise))
+  pairwise_mat[lower.tri(pairwise_mat)] = pair_means
+  pairwise_mat = pairwise_mat + t(pairwise_mat)
+  if(is_continuous) pairwise_mat = -0.5 * pairwise_mat
+  out$pairwise = pairwise_mat
+
+  # Main: GGM has no main effects; OMRF/mixed don't (yet) support
+  # hierarchical V-correction. Carry the existing field through for shape
+  # parity with the field of the same name on the fit object. Use
+  # `out["main"] = list(NULL)` so the slot name survives even when the
+  # value is NULL (plain `out$main = NULL` would drop the slot).
+  out["main"] = list(fit$posterior_mean_main)
+
+  # Indicator (when edge selection was active).
+  if(!is.null(raw$indicator)) {
+    ind_means = weighted_col_means(raw$indicator)
+    indicator_mat = matrix(0, num_variables, num_variables,
+                           dimnames = dimnames(fit$posterior_mean_indicator))
+    indicator_mat[lower.tri(indicator_mat)] = ind_means
+    indicator_mat = indicator_mat + t(indicator_mat)
+    out$indicator = indicator_mat
+  } else {
+    out["indicator"] = list(NULL)
+  }
+
+  # Residual variance for GGM: sign-correct on the per-sample functional
+  # 1/K_ii (matches the per-sample-inversion form used by the plain field
+  # to avoid Jensen bias).
+  if(is_continuous && !is.null(raw$main)) {
+    inv_main_per_chain = lapply(raw$main, function(m) 1 / m)
+    rv = weighted_col_means(inv_main_per_chain)
+    names(rv) = names(fit$posterior_mean_residual_variance)
+    out$residual_variance = rv
+  } else {
+    out["residual_variance"] = list(NULL)
+  }
+
+  out
+}