diff --git a/DESCRIPTION b/DESCRIPTION
index a5c0ddd0..701d9a2f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -28,14 +28,15 @@ Copyright: Includes datasets 'ADHD' and 'Boredom', which are licensed under CC-B
 License: GPL (>= 2)
 URL: https://Bayesian-Graphical-Modelling-Lab.github.io/bgms/, https://github.com/Bayesian-Graphical-Modelling-Lab/bgms
 BugReports: https://github.com/Bayesian-Graphical-Modelling-Lab/bgms/issues
-Imports: 
-    Rcpp (>= 1.0.7), 
-    RcppParallel, 
-    Rdpack, 
+Imports:
+    Rcpp (>= 1.0.7),
+    RcppParallel,
+    Rdpack,
     S7,
-    methods, 
+    methods,
     lifecycle,
-    stats
+    stats,
+    utils
 RdMacros: Rdpack
 LinkingTo: 
     Rcpp, 
diff --git a/NAMESPACE b/NAMESPACE
index bfb05b15..c7e82c2f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -73,6 +73,8 @@ export(gamma_prior)
 export(mrfSampler)
 export(normal_prior)
 export(sample_ggm_prior)
+export(sample_graph_prior)
+export(sample_sbm_prior)
 export(sbm_prior)
 export(simulate_mrf)
 export(summarize_zratio_diagnostics)
diff --git a/NEWS.md b/NEWS.md
index 0f668e1c..3087e7db 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -12,6 +12,9 @@
 
 ## New features
 
+* `sample_graph_prior()`: draws edge-inclusion indicators, together with any edge-prior hyperparameters, from the graph level of the spike-and-slab prior. Under the hierarchical specification the draw is ancestral and exact; under the joint specification it runs the zero-data prior chain, so the draws carry the per-graph normalizer tilt. Hyperparameters can be fixed instead of sampled (`theta` for Bernoulli and Beta-Bernoulli priors, `allocations` plus `block_probs` for the Stochastic-Block prior).
+* `sample_sbm_prior()`: ancestral draws of block allocations and pair-inclusion probabilities from the MFM-SBM edge-prior hyperprior.
+* The edge-prior correction table now announces itself only when it actually builds (a cache hit is silent) and states that the build is one-time and cached. In interactive sessions it draws a progress bar for both serial and parallel builds (a parallel build advances the bar as each batch of forked chains completes). The bar follows `display_progress` (the same control as the sampler's bar), not the advisory `bgms.verbose` flag.
 * Gaussian graphical models (GGM): `bgm(x, variable_type = "continuous")` fits a GGM with Bayesian edge selection. Sampling uses NUTS on a free-element Cholesky (theta-space) parameterization of the precision matrix, which keeps the precision matrix positive-definite by construction; adaptive-metropolis is also available.
 * GGM Gibbs sampler: `update_method = "gibbs"` fits a GGM with a conjugate row-by-row update of the precision matrix, with or without edge selection. It needs no step-size or proposal tuning and supports a Normal or Cauchy prior on the edges. Continuous data only.
 * Gibbs warmup staging: with `update_method = "gibbs"` and edge selection, the first 15% of the warmup runs the full model so the precision matrix settles, and edge selection is active for the remaining 85%. Previously edge selection only started at the first retained iteration, so the graph's equilibration happened inside the retained samples. Both windows scale with the warmup budget; the warmup default is unchanged.
diff --git a/R/bgm.R b/R/bgm.R
index 92f5d6de..e00f81e8 100644
--- a/R/bgm.R
+++ b/R/bgm.R
@@ -286,7 +286,9 @@
 #' @param chains Integer. Number of parallel chains to run. Default: \code{4}.
 #'
 #' @param cores Integer. Number of CPU cores for parallel execution.
-#'   Default: \code{parallel::detectCores()}.
+#'   Sampling uses \code{min(cores, chains)}; some computations outside of
+#'   sampling (such as building the edge-selection prior correction table)
+#'   use all \code{cores}. Default: \code{parallel::detectCores()}.
 #'
 #' @param seed Optional integer. Random seed for reproducibility. Must be a
 #'   single non-negative integer.
diff --git a/R/correction_tables.R b/R/correction_tables.R
index 555d9326..18d94800 100644
--- a/R/correction_tables.R
+++ b/R/correction_tables.R
@@ -153,6 +153,92 @@ normalize_builder_cores = function(cores) {
 }
 
 
+# ------------------------------------------------------------------
+# new_correction_progress (internal)
+# ------------------------------------------------------------------
+# In-place progress bar matching the bgms MCMC bar (see the theme and
+# bar-width logic in src/utils/progress_manager.cpp): tortoise-shell
+# brackets, a heavy horizontal rule filled in blue and empty in gray,
+# a sub-cell partial glyph, and a "prefix: <bar> cur/tot (xx.x%)"
+# layout redrawn with a carriage return. Unicode + ANSI colour when the
+# session is UTF-8, an ASCII "[=== ]" fallback otherwise. Returns
+# update(current) and close(); the caller decides whether to draw it.
+# ------------------------------------------------------------------
+new_correction_progress = function(total, prefix = "Correction table") {
+  unicode = isTRUE(l10n_info()[["UTF-8"]])
+  is_rstudio = Sys.getenv("RSTUDIO") == "1"
+
+  # Bar width, mirroring progress_manager.cpp so the bar does not wrap.
+  console_width = if(is_rstudio) {
+    max(0L, as.integer(getOption("width", 80L))) + 3L
+  } else {
+    80L
+  }
+  line_width = if(is_rstudio) {
+    max(10L, min(console_width - 25L, 70L))
+  } else {
+    70L
+  }
+  bar_width = if(line_width <= 5L) {
+    0L
+  } else if(line_width < 20L) {
+    line_width - 10L
+  } else if(line_width < 40L) {
+    line_width - 15L
+  } else if(line_width > 70L) {
+    40L
+  } else {
+    line_width - 30L
+  }
+  if(is_rstudio) {
+    bar_width = if(bar_width > 30L) bar_width - 20L else 10L
+  }
+
+  if(unicode) {
+    # Tortoise-shell brackets (U+2997/U+2998), heavy horizontal rule
+    # (U+2501) filled/empty, heavy sub-cell (U+257A); ANSI 38;5;73 blue
+    # and 37 gray, matching progress_manager.cpp.
+    lhs = "\u2997"
+    rhs = "\u2998"
+    filled = "\u001b[38;5;73m\u2501\u001b[39m"
+    partial_more = filled
+    partial_less = "\u001b[37m\u257a\u001b[39m"
+    empty = "\u001b[37m\u2501\u001b[39m"
+  } else {
+    lhs = "["
+    rhs = "]"
+    filled = "="
+    partial_more = " "
+    partial_less = " "
+    empty = " "
+  }
+
+  draw = function(current) {
+    frac = if(total > 0L) current / total else 1
+    exact = frac * bar_width
+    n_filled = min(as.integer(exact), bar_width)
+    bar = strrep(filled, n_filled)
+    if(n_filled < bar_width) {
+      part = exact - n_filled
+      if(part > 0) {
+        bar = paste0(bar, if(part > 0.5) partial_more else partial_less)
+        n_filled = n_filled + 1L
+      }
+    }
+    if(n_filled < bar_width) {
+      bar = paste0(bar, strrep(empty, bar_width - n_filled))
+    }
+    cat(sprintf(
+      "\r%s: %s%s%s %d/%d (%.1f%%)", prefix, lhs, bar, rhs,
+      current, total, 100 * frac
+    ))
+    utils::flush.console()
+  }
+
+  list(update = draw, close = function() cat("\n"))
+}
+
+
 # ------------------------------------------------------------------
 # sweep_prior_edge_density (internal)
 # ------------------------------------------------------------------
@@ -166,7 +252,8 @@ sweep_prior_edge_density = function(p, theta, delta,
                                     n_samples = 2000L, n_warmup = 500L,
                                     n_seeds = 3L,
                                     update_method = "gibbs",
-                                    cores = 1L, base_seed = 1L) {
+                                    cores = 1L, base_seed = 1L,
+                                    show_progress = FALSE) {
   cores = normalize_builder_cores(cores)
   num_pairs = p * (p - 1) / 2
   cells = expand.grid(theta = theta, seed = seq_len(n_seeds))
@@ -185,13 +272,42 @@ sweep_prior_edge_density = function(p, theta, delta,
     mean(draws$edge_indicators)
   }
 
+  # The bar follows display_progress (like the sampler's own bar), not the
+  # advisory bgms.verbose flag. It redraws in place with a carriage return, so
+  # it is drawn only in an interactive session; batch runs (scripts, R CMD
+  # check, the test suite) rely on the one-line build announcement in
+  # ggm_correction_table(). A parallel sweep forks the cells in batches of
+  # `cores` and advances the bar as each batch of chains completes.
+  n_cells = nrow(cells)
+  draw_bar = show_progress && interactive()
+  pb = NULL
+  if(draw_bar) {
+    pb = new_correction_progress(n_cells)
+    on.exit(pb$close(), add = TRUE)
+    pb$update(0L)
+  }
   edens_cells = if(cores > 1L) {
-    unlist(parallel::mclapply(
-      seq_len(nrow(cells)), one_cell,
-      mc.cores = cores, mc.preschedule = FALSE
-    ))
+    out = numeric(n_cells)
+    done = 0L
+    for(batch in split(seq_len(n_cells), ceiling(seq_len(n_cells) / cores))) {
+      out[batch] = as.numeric(unlist(parallel::mclapply(
+        batch, one_cell,
+        mc.cores = cores, mc.preschedule = FALSE
+      )))
+      done = done + length(batch)
+      if(draw_bar) {
+        pb$update(done)
+      }
+    }
+    out
   } else {
-    vapply(seq_len(nrow(cells)), one_cell, numeric(1))
+    vapply(seq_len(n_cells), function(k) {
+      v = one_cell(k)
+      if(draw_bar) {
+        pb$update(k)
+      }
+      v
+    }, numeric(1))
   }
   if(anyNA(edens_cells)) {
     stop("Correction-table sweep: a prior chain failed.")
@@ -242,7 +358,7 @@ build_ggm_correction_table = function(
   precision_scale_prior = gamma_prior(shape = 1, eta = 1),
   n_grid = 120L, n_samples = 2000L, n_warmup = 500L, n_seeds = 3L,
   update_method = c("gibbs", "adaptive-metropolis"),
-  cores = 1L, base_seed = 1L
+  cores = 1L, base_seed = 1L, show_progress = FALSE
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -255,7 +371,8 @@ build_ggm_correction_table = function(
     interaction_prior = interaction_prior,
     precision_scale_prior = precision_scale_prior,
     n_samples = n_samples, n_warmup = n_warmup, n_seeds = n_seeds,
-    update_method = update_method, cores = cores, base_seed = base_seed
+    update_method = update_method, cores = cores, base_seed = base_seed,
+    show_progress = show_progress
   )
 
   table = correction_table_from_edens(
@@ -371,18 +488,20 @@ ggm_edge_prior_correction = function(prior, sampler, num_variables,
   interaction_prior = correction_interaction_prior(prior)
   precision_scale_prior = correction_scale_prior(prior)
 
-  if(isTRUE(sampler$verbose)) {
-    message(
-      "Edge-prior correction: building or loading the normalizing-constant ",
-      "table for this model (cached across fits)."
-    )
-  }
+  # The build announcement follows the advisory bgms.verbose flag; the progress
+  # bar follows display_progress (progress_type 0 is "none"), matching the
+  # sampler's own bar.
+  show_progress = is.null(sampler$progress_type) ||
+    !identical(as.integer(sampler$progress_type), 0L)
+
   table = ggm_correction_table(
     p = num_continuous, delta = prior$delta,
     interaction_prior = interaction_prior,
     precision_scale_prior = precision_scale_prior,
     update_method = "gibbs",
-    cores = sampler$cores
+    cores = sampler$cores,
+    verbose = isTRUE(sampler$verbose),
+    show_progress = show_progress
   )
   if(identical(prior$edge_prior, "Stochastic-Block") &&
     is.null(table$fprime)) {
@@ -420,7 +539,8 @@ ggm_correction_table = function(
   precision_scale_prior = gamma_prior(shape = 1, eta = 1),
   n_grid = 120L, n_samples = 2000L, n_warmup = 500L, n_seeds = 3L,
   update_method = c("gibbs", "adaptive-metropolis"),
-  cores = 1L, base_seed = 1L, refresh = FALSE
+  cores = 1L, base_seed = 1L, refresh = FALSE, verbose = FALSE,
+  show_progress = FALSE
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -451,13 +571,19 @@ ggm_correction_table = function(
     }
   }
 
+  if(verbose) {
+    message(
+      "Building the edge-selection prior correction table (one-time for ",
+      "this model size and prior; cached for later fits)."
+    )
+  }
   table = build_ggm_correction_table(
     p = p, delta = delta,
     interaction_prior = interaction_prior,
     precision_scale_prior = precision_scale_prior,
     n_grid = n_grid, n_samples = n_samples, n_warmup = n_warmup,
     n_seeds = n_seeds, update_method = update_method,
-    cores = cores, base_seed = base_seed
+    cores = cores, base_seed = base_seed, show_progress = show_progress
   )
 
   if(use_cache) {
diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R
index 18cffb4e..7311abd6 100644
--- a/R/sample_ggm_prior.R
+++ b/R/sample_ggm_prior.R
@@ -354,7 +354,8 @@ sample_ggm_prior = function(
       p = p, delta = delta,
       interaction_prior = interaction_prior,
       precision_scale_prior = precision_scale_prior,
-      update_method = "gibbs"
+      update_method = "gibbs",
+      verbose = isTRUE(verbose)
     )
     correction = correction_list_from_table(table, ep$edge_prior)
   }
diff --git a/R/sample_graph_prior.R b/R/sample_graph_prior.R
new file mode 100644
index 00000000..b15dae19
--- /dev/null
+++ b/R/sample_graph_prior.R
@@ -0,0 +1,427 @@
+# ==============================================================================
+# Graph-level prior samplers
+# ==============================================================================
+#
+# sample_graph_prior(): draws (hyperparameters, edge indicators) from the
+# graph level of the spike-and-slab prior, under either specification of
+# how the precision prior composes with the graph:
+#
+#   hierarchical: p(hyper) . pi(Gamma | hyper) -- ancestral and exact
+#   joint:        p(hyper) . q(Gamma | hyper), with
+#                 q(Gamma | hyper) proportional to Z(Gamma) pi(Gamma | hyper)
+#                 -- the determinant-tilted law, sampled by the zero-data
+#                 (K, Gamma) chain with K discarded from the output
+#
+# sample_sbm_prior(): ancestral draws of (allocations, pair probabilities)
+# from the MFM-SBM edge-prior hyperprior.
+#
+# Edge indicators are returned in row-major upper-triangle order (i < j),
+# matching the K_offdiag column order of sample_ggm_prior().
+# ==============================================================================
+
+
+# Run fn() inside an RNG scope keyed to `seed`, restoring the caller's RNG
+# state on exit (same idiom as ggm_prior_ancestral_indicators).
+with_graph_prior_seed = function(seed, fn) {
+  has_seed = exists(".Random.seed", envir = globalenv(), inherits = FALSE)
+  if(has_seed) {
+    old_seed = get(".Random.seed", envir = globalenv(), inherits = FALSE)
+    on.exit(assign(".Random.seed", old_seed, envir = globalenv()), add = TRUE)
+  } else {
+    on.exit(
+      if(exists(".Random.seed", envir = globalenv(), inherits = FALSE)) {
+        rm(list = ".Random.seed", envir = globalenv())
+      },
+      add = TRUE
+    )
+  }
+  set.seed(seed)
+  fn()
+}
+
+
+# Row-major upper-triangle pair indices (i < j): a (E x 2) matrix with
+# E = p(p-1)/2 rows, ordered (1,2), (1,3), ..., (1,p), (2,3), ...
+graph_pair_indices = function(p) {
+  ii = rep(seq_len(p - 1L), times = (p - 1L):1L)
+  jj = unlist(lapply(seq_len(p - 1L), function(i) (i + 1L):p))
+  cbind(ii, jj)
+}
+
+
+# Resolve the optional conditioning arguments to a fixed p x p pair
+# probability matrix, or NULL when no conditioning was requested.
+graph_prior_conditioning = function(ep, p, theta, allocations, block_probs) {
+  if(!is.null(theta) && (!is.null(allocations) || !is.null(block_probs))) {
+    stop(
+      "Supply either 'theta' or ('allocations', 'block_probs') to condition ",
+      "the graph prior, not both."
+    )
+  }
+  if(!is.null(theta)) {
+    if(identical(ep$edge_prior, "Stochastic-Block")) {
+      stop(
+        "Condition the Stochastic-Block prior with 'allocations' and ",
+        "'block_probs'; 'theta' conditions the Bernoulli and Beta-Bernoulli ",
+        "priors."
+      )
+    }
+    if(!is.numeric(theta) || length(theta) != 1L || is.na(theta) ||
+      theta <= 0 || theta >= 1) {
+      stop("'theta' must be a single numeric in (0, 1).")
+    }
+    return(matrix(theta, p, p))
+  }
+  if(!is.null(allocations) || !is.null(block_probs)) {
+    if(is.null(allocations) || is.null(block_probs)) {
+      stop(
+        "Conditioning the Stochastic-Block prior needs both 'allocations' ",
+        "and 'block_probs'."
+      )
+    }
+    if(!identical(ep$edge_prior, "Stochastic-Block")) {
+      stop(
+        "'allocations'/'block_probs' condition the Stochastic-Block prior; ",
+        "use edge_prior = sbm_prior()."
+      )
+    }
+    allocations = as.integer(allocations)
+    if(length(allocations) != p || anyNA(allocations) ||
+      any(allocations < 1L)) {
+      stop("'allocations' must be a length-p vector of 1-based block labels.")
+    }
+    if(!is.matrix(block_probs) || nrow(block_probs) != ncol(block_probs) ||
+      !isSymmetric(unname(block_probs)) ||
+      any(block_probs <= 0) || any(block_probs >= 1)) {
+      stop(
+        "'block_probs' must be a symmetric matrix with entries in (0, 1)."
+      )
+    }
+    if(max(allocations) > nrow(block_probs)) {
+      stop("'allocations' labels exceed the size of 'block_probs'.")
+    }
+    prob = matrix(0.5, p, p)
+    for(i in seq_len(p - 1L)) {
+      for(j in (i + 1L):p) {
+        prob[i, j] = block_probs[allocations[i], allocations[j]]
+        prob[j, i] = prob[i, j]
+      }
+    }
+    return(prob)
+  }
+  NULL
+}
+
+
+#' @title Sample from the Graph Prior
+#'
+#' @description
+#' Draws edge-inclusion indicators, together with any edge-prior
+#' hyperparameters, from the graph level of the spike-and-slab prior used by
+#' \code{\link{bgm}} for models with continuous variables. The \code{spec}
+#' argument selects how the precision prior composes with the graph:
+#' \itemize{
+#'   \item \code{"hierarchical"} (default): the graph marginal is exactly the
+#'     edge prior, \eqn{p(\mathrm{hyper}) \, \pi(\Gamma \mid \mathrm{hyper})}.
+#'     Sampling is ancestral and exact: hyperparameters from their prior,
+#'     then independent pair flips.
+#'   \item \code{"joint"}: the graph marginal is reweighted by the per-graph
+#'     normalizer of the determinant-tilted precision prior,
+#'     \eqn{q(\Gamma \mid \mathrm{hyper}) \propto Z(\Gamma) \,
+#'     \pi(\Gamma \mid \mathrm{hyper})}. Sampling runs the zero-data
+#'     \eqn{(K, \Gamma)} chain of \code{\link{sample_ggm_prior}} and discards
+#'     \eqn{K}; with a Beta-Bernoulli or Stochastic-Block prior the
+#'     hyperparameter updates apply the normalizing-constant correction, so
+#'     the first call for a model cell may build the correction table
+#'     (cached across fits).
+#' }
+#'
+#' @details
+#' The optional conditioning arguments fix the edge-prior hyperparameters
+#' instead of sampling them: \code{theta} fixes the inclusion probability of
+#' a Bernoulli or Beta-Bernoulli prior, and \code{allocations} plus
+#' \code{block_probs} fix the block structure of a Stochastic-Block prior
+#' (reducing it to independent pair flips at the given block probabilities).
+#'
+#' Under \code{spec = "joint"} the tilted graph law depends on the precision
+#' prior through its normalizer, so \code{interaction_prior},
+#' \code{precision_scale_prior}, and \code{delta} are part of the graph law;
+#' they are ignored under \code{spec = "hierarchical"}.
+#'
+#' @param p Integer. Number of nodes (\eqn{p \ge 2}).
+#' @param n_samples Integer. Number of prior draws.
+#' @param edge_prior An edge prior specification object:
+#'   \code{\link{bernoulli_prior}()}, \code{\link{beta_bernoulli_prior}()},
+#'   or \code{\link{sbm_prior}()}. Default \code{bernoulli_prior(0.5)}.
+#' @param spec One of \code{"hierarchical"} (default) or \code{"joint"}.
+#' @param interaction_prior A \code{\link{cauchy_prior}()} or
+#'   \code{\link{normal_prior}()} for the pairwise (slab) part of the
+#'   precision prior. Used only when \code{spec = "joint"}.
+#' @param precision_scale_prior A \code{\link{gamma_prior}()} or
+#'   \code{\link{exponential_prior}()} for the precision diagonal. Used only
+#'   when \code{spec = "joint"}.
+#' @param delta Non-negative numeric or \code{NULL} (default): determinant
+#'   tilt exponent; \code{NULL} resolves to \eqn{0.5 \log(p)}. Used only
+#'   when \code{spec = "joint"}.
+#' @param theta Optional numeric in (0, 1): fix the inclusion probability of
+#'   a Bernoulli or Beta-Bernoulli edge prior instead of sampling it.
+#' @param allocations Optional integer vector of length \code{p} with 1-based
+#'   block labels: fix the Stochastic-Block allocation instead of sampling it.
+#'   Requires \code{block_probs}.
+#' @param block_probs Optional symmetric matrix with entries in (0, 1): the
+#'   block-pair inclusion probabilities that go with \code{allocations}.
+#' @param n_warmup Integer. Warmup iterations of the zero-data chain. Used
+#'   only when \code{spec = "joint"}. Default \code{2e3}.
+#' @param seed Integer. Seed for the draw; the caller's RNG state is
+#'   restored on exit.
+#' @param verbose Logical. Print progress of the zero-data chain and of a
+#'   correction-table build. Default \code{TRUE}.
+#'
+#' @return A list with elements:
+#'   \describe{
+#'     \item{\code{edge_indicators}}{Integer matrix
+#'       (\code{n_samples x p(p-1)/2}) of edge-inclusion indicators, columns
+#'       in row-major upper-triangle order (matching
+#'       \code{sample_ggm_prior()}'s \code{K_offdiag}).}
+#'     \item{\code{pair_names}}{Character vector labeling the columns as
+#'       \code{"i-j"}.}
+#'     \item{\code{theta}}{Only with an unconditioned
+#'       \code{beta_bernoulli_prior()}: numeric vector of sampled inclusion
+#'       probabilities.}
+#'     \item{\code{allocations}}{Only with an unconditioned
+#'       \code{sbm_prior()}: integer matrix (\code{n_samples x p}) of sampled
+#'       block allocations.}
+#'     \item{\code{spec}, \code{edge_prior}, \code{p}}{The specification,
+#'       edge-prior family, and node count of the draw.}
+#'   }
+#'
+#' @examples
+#' # Hierarchical spec: the graph marginal is exactly the edge prior.
+#' g = sample_graph_prior(
+#'   p = 6, n_samples = 200,
+#'   edge_prior = bernoulli_prior(0.3), seed = 11
+#' )
+#' mean(g$edge_indicators) # about 0.3
+#'
+#' # Beta-Bernoulli: inclusion probabilities are sampled alongside.
+#' g = sample_graph_prior(
+#'   p = 6, n_samples = 200,
+#'   edge_prior = beta_bernoulli_prior(2, 4), seed = 11
+#' )
+#' mean(g$theta) # about 1/3
+#'
+#' \donttest{
+#' # Joint spec: the graph law carries the per-graph normalizer Z(Gamma).
+#' g = sample_graph_prior(
+#'   p = 6, n_samples = 500,
+#'   edge_prior = bernoulli_prior(0.3), spec = "joint",
+#'   interaction_prior = normal_prior(scale = 0.5),
+#'   precision_scale_prior = gamma_prior(shape = 1, rate = 2),
+#'   seed = 11, verbose = FALSE
+#' )
+#' mean(g$edge_indicators) # shifted away from 0.3 by the Z(Gamma) tilt
+#' }
+#' @seealso \code{\link{sample_ggm_prior}}, \code{\link{sample_sbm_prior}},
+#'   \code{\link{bernoulli_prior}}, \code{\link{beta_bernoulli_prior}},
+#'   \code{\link{sbm_prior}}, \code{\link{bgm}}
+#'
+#' @export
+sample_graph_prior = function(
+  p,
+  n_samples,
+  edge_prior = bernoulli_prior(0.5),
+  spec = c("hierarchical", "joint"),
+  interaction_prior = cauchy_prior(scale = 2.5),
+  precision_scale_prior = gamma_prior(shape = 1, eta = 1),
+  delta = NULL,
+  theta = NULL,
+  allocations = NULL,
+  block_probs = NULL,
+  n_warmup = 2e3,
+  seed = 1L,
+  verbose = TRUE
+) {
+  spec = match.arg(spec)
+  if(!is.numeric(p) || length(p) != 1L || is.na(p) || p < 2 || p != round(p)) {
+    stop("'p' must be a single integer >= 2.")
+  }
+  p = as.integer(p)
+  if(!is.numeric(n_samples) || length(n_samples) != 1L || is.na(n_samples) ||
+    n_samples < 1 || n_samples != round(n_samples)) {
+    stop("'n_samples' must be a single positive integer.")
+  }
+  n_samples = as.integer(n_samples)
+
+  ep = unpack_indicator_prior(edge_prior, num_variables = p)
+  cond_prob = graph_prior_conditioning(ep, p, theta, allocations, block_probs)
+
+  pairs = graph_pair_indices(p)
+  pair_names = paste0(pairs[, 1L], "-", pairs[, 2L])
+
+  if(spec == "joint") {
+    eff_edge_prior = if(is.null(cond_prob)) {
+      edge_prior
+    } else {
+      bernoulli_prior(inclusion_probability = cond_prob)
+    }
+    draws = sample_ggm_prior(
+      p = p, n_samples = n_samples, n_warmup = as.integer(n_warmup),
+      interaction_prior = interaction_prior,
+      precision_scale_prior = precision_scale_prior,
+      delta = delta,
+      spec = "joint",
+      edge_prior = eff_edge_prior,
+      update_method = "gibbs",
+      seed = as.integer(seed),
+      verbose = verbose
+    )
+    out = list(
+      edge_indicators = draws$edge_indicators,
+      pair_names = pair_names,
+      theta = draws$theta,
+      allocations = draws$allocations,
+      spec = spec,
+      edge_prior = ep$edge_prior,
+      p = p
+    )
+    return(out)
+  }
+
+  # Hierarchical spec: ancestral. Hyperparameters from their prior, then
+  # independent pair flips given the implied pair probabilities.
+  res = with_graph_prior_seed(seed, function() {
+    n_edges = nrow(pairs)
+    theta_out = NULL
+    alloc_out = NULL
+
+    if(!is.null(cond_prob) || identical(ep$edge_prior, "Bernoulli")) {
+      prob = if(is.null(cond_prob)) ep$inclusion_probability else cond_prob
+      pv = prob[pairs]
+      ru = matrix(runif(n_samples * n_edges), n_samples, n_edges)
+      indicators = 1L * (ru < matrix(pv, n_samples, n_edges, byrow = TRUE))
+    } else if(identical(ep$edge_prior, "Beta-Bernoulli")) {
+      theta_out = rbeta(
+        n_samples, ep$beta_bernoulli_alpha, ep$beta_bernoulli_beta
+      )
+      ru = matrix(runif(n_samples * n_edges), n_samples, n_edges)
+      indicators = 1L * (ru < theta_out)
+    } else {
+      alloc_out = matrix(0L, n_samples, p)
+      indicators = matrix(0L, n_samples, n_edges)
+      for(s in seq_len(n_samples)) {
+        z = ancestral_mfm_sbm_partition(p, ep$lambda, ep$dirichlet_alpha)
+        prob = ancestral_sbm_pair_probabilities(
+          z,
+          ep$beta_bernoulli_alpha, ep$beta_bernoulli_beta,
+          ep$beta_bernoulli_alpha_between, ep$beta_bernoulli_beta_between
+        )
+        alloc_out[s, ] = as.integer(z)
+        indicators[s, ] = as.integer(runif(n_edges) < prob[pairs])
+      }
+    }
+    storage.mode(indicators) = "integer"
+    list(indicators = indicators, theta = theta_out, allocations = alloc_out)
+  })
+
+  list(
+    edge_indicators = res$indicators,
+    pair_names = pair_names,
+    theta = res$theta,
+    allocations = res$allocations,
+    spec = spec,
+    edge_prior = ep$edge_prior,
+    p = p
+  )
+}
+
+
+#' @title Sample from the Stochastic-Block Edge-Prior Hyperprior
+#'
+#' @description
+#' Ancestral draws from the MFM-SBM hyperprior of
+#' \code{\link{sbm_prior}()}: a partition of the nodes into blocks
+#' (shifted-Poisson number of components, Dirichlet-weighted allocation) and
+#' Beta-distributed within- and between-block edge-inclusion probabilities.
+#' These are the hyperparameters that \code{\link{sample_graph_prior}} and
+#' \code{\link{bgm}} integrate over when the edge prior is a Stochastic-Block
+#' prior.
+#'
+#' @param p Integer. Number of nodes (\eqn{p \ge 2}).
+#' @param n_samples Integer. Number of prior draws.
+#' @param edge_prior An \code{\link{sbm_prior}()} object. Default
+#'   \code{sbm_prior()}.
+#' @param seed Integer. Seed for the draw; the caller's RNG state is
+#'   restored on exit.
+#'
+#' @return A list with elements:
+#'   \describe{
+#'     \item{\code{allocations}}{Integer matrix (\code{n_samples x p}) of
+#'       block labels.}
+#'     \item{\code{pair_probability}}{Numeric matrix
+#'       (\code{n_samples x p(p-1)/2}) of implied pair-inclusion
+#'       probabilities, columns in row-major upper-triangle order.}
+#'     \item{\code{num_blocks}}{Integer vector: number of occupied blocks
+#'       per draw.}
+#'     \item{\code{pair_names}}{Character vector labeling the pair columns
+#'       as \code{"i-j"}.}
+#'     \item{\code{p}}{The node count.}
+#'   }
+#'
+#' @examples
+#' draws = sample_sbm_prior(p = 8, n_samples = 100, seed = 4)
+#' table(draws$num_blocks)
+#' range(draws$pair_probability)
+#' @seealso \code{\link{sbm_prior}}, \code{\link{sample_graph_prior}},
+#'   \code{\link{bgm}}
+#'
+#' @export
+sample_sbm_prior = function(p, n_samples, edge_prior = sbm_prior(),
+                            seed = 1L) {
+  if(!is.numeric(p) || length(p) != 1L || is.na(p) || p < 2 || p != round(p)) {
+    stop("'p' must be a single integer >= 2.")
+  }
+  p = as.integer(p)
+  if(!is.numeric(n_samples) || length(n_samples) != 1L || is.na(n_samples) ||
+    n_samples < 1 || n_samples != round(n_samples)) {
+    stop("'n_samples' must be a single positive integer.")
+  }
+  n_samples = as.integer(n_samples)
+
+  ep = unpack_indicator_prior(edge_prior, num_variables = p)
+  if(!identical(ep$edge_prior, "Stochastic-Block")) {
+    stop("'edge_prior' must be an sbm_prior() object.")
+  }
+
+  pairs = graph_pair_indices(p)
+  pair_names = paste0(pairs[, 1L], "-", pairs[, 2L])
+
+  res = with_graph_prior_seed(seed, function() {
+    allocations = matrix(0L, n_samples, p)
+    pair_probability = matrix(0, n_samples, nrow(pairs))
+    num_blocks = integer(n_samples)
+    for(s in seq_len(n_samples)) {
+      z = ancestral_mfm_sbm_partition(p, ep$lambda, ep$dirichlet_alpha)
+      prob = ancestral_sbm_pair_probabilities(
+        z,
+        ep$beta_bernoulli_alpha, ep$beta_bernoulli_beta,
+        ep$beta_bernoulli_alpha_between, ep$beta_bernoulli_beta_between
+      )
+      allocations[s, ] = as.integer(z)
+      pair_probability[s, ] = prob[pairs]
+      num_blocks[s] = length(unique(z))
+    }
+    list(
+      allocations = allocations,
+      pair_probability = pair_probability,
+      num_blocks = num_blocks
+    )
+  })
+
+  list(
+    allocations = res$allocations,
+    pair_probability = res$pair_probability,
+    num_blocks = res$num_blocks,
+    pair_names = pair_names,
+    p = p
+  )
+}
diff --git a/man/bgm.Rd b/man/bgm.Rd
index 9b9ea71f..c2d70311 100644
--- a/man/bgm.Rd
+++ b/man/bgm.Rd
@@ -259,7 +259,9 @@ matrix. Default: \code{TRUE}.}
 \item{chains}{Integer. Number of parallel chains to run. Default: \code{4}.}
 
 \item{cores}{Integer. Number of CPU cores for parallel execution.
-Default: \code{parallel::detectCores()}.}
+Sampling uses \code{min(cores, chains)}; some computations outside of
+sampling (such as building the edge-selection prior correction table)
+use all \code{cores}. Default: \code{parallel::detectCores()}.}
 
 \item{display_progress}{Character. Controls progress reporting during
 sampling. Options: \code{"per-chain"} (separate bar per chain),
diff --git a/man/sample_graph_prior.Rd b/man/sample_graph_prior.Rd
new file mode 100644
index 00000000..fce4a00c
--- /dev/null
+++ b/man/sample_graph_prior.Rd
@@ -0,0 +1,148 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sample_graph_prior.R
+\name{sample_graph_prior}
+\alias{sample_graph_prior}
+\title{Sample from the Graph Prior}
+\usage{
+sample_graph_prior(
+  p,
+  n_samples,
+  edge_prior = bernoulli_prior(0.5),
+  spec = c("hierarchical", "joint"),
+  interaction_prior = cauchy_prior(scale = 2.5),
+  precision_scale_prior = gamma_prior(shape = 1, eta = 1),
+  delta = NULL,
+  theta = NULL,
+  allocations = NULL,
+  block_probs = NULL,
+  n_warmup = 2000,
+  seed = 1L,
+  verbose = TRUE
+)
+}
+\arguments{
+\item{p}{Integer. Number of nodes (\eqn{p \ge 2}).}
+
+\item{n_samples}{Integer. Number of prior draws.}
+
+\item{edge_prior}{An edge prior specification object:
+\code{\link{bernoulli_prior}()}, \code{\link{beta_bernoulli_prior}()},
+or \code{\link{sbm_prior}()}. Default \code{bernoulli_prior(0.5)}.}
+
+\item{spec}{One of \code{"hierarchical"} (default) or \code{"joint"}.}
+
+\item{interaction_prior}{A \code{\link{cauchy_prior}()} or
+\code{\link{normal_prior}()} for the pairwise (slab) part of the
+precision prior. Used only when \code{spec = "joint"}.}
+
+\item{precision_scale_prior}{A \code{\link{gamma_prior}()} or
+\code{\link{exponential_prior}()} for the precision diagonal. Used only
+when \code{spec = "joint"}.}
+
+\item{delta}{Non-negative numeric or \code{NULL} (default): determinant
+tilt exponent; \code{NULL} resolves to \eqn{0.5 \log(p)}. Used only
+when \code{spec = "joint"}.}
+
+\item{theta}{Optional numeric in (0, 1): fix the inclusion probability of
+a Bernoulli or Beta-Bernoulli edge prior instead of sampling it.}
+
+\item{allocations}{Optional integer vector of length \code{p} with 1-based
+block labels: fix the Stochastic-Block allocation instead of sampling it.
+Requires \code{block_probs}.}
+
+\item{block_probs}{Optional symmetric matrix with entries in (0, 1): the
+block-pair inclusion probabilities that go with \code{allocations}.}
+
+\item{n_warmup}{Integer. Warmup iterations of the zero-data chain. Used
+only when \code{spec = "joint"}. Default \code{2e3}.}
+
+\item{seed}{Integer. Seed for the draw; the caller's RNG state is
+restored on exit.}
+
+\item{verbose}{Logical. Print progress of the zero-data chain and of a
+correction-table build. Default \code{TRUE}.}
+}
+\value{
+A list with elements:
+\describe{
+\item{\code{edge_indicators}}{Integer matrix
+(\code{n_samples x p(p-1)/2}) of edge-inclusion indicators, columns
+in row-major upper-triangle order (matching
+\code{sample_ggm_prior()}'s \code{K_offdiag}).}
+\item{\code{pair_names}}{Character vector labeling the columns as
+\code{"i-j"}.}
+\item{\code{theta}}{Only with an unconditioned
+\code{beta_bernoulli_prior()}: numeric vector of sampled inclusion
+probabilities.}
+\item{\code{allocations}}{Only with an unconditioned
+\code{sbm_prior()}: integer matrix (\code{n_samples x p}) of sampled
+block allocations.}
+\item{\code{spec}, \code{edge_prior}, \code{p}}{The specification,
+edge-prior family, and node count of the draw.}
+}
+}
+\description{
+Draws edge-inclusion indicators, together with any edge-prior
+hyperparameters, from the graph level of the spike-and-slab prior used by
+\code{\link{bgm}} for models with continuous variables. The \code{spec}
+argument selects how the precision prior composes with the graph:
+\itemize{
+\item \code{"hierarchical"} (default): the graph marginal is exactly the
+edge prior, \eqn{p(\mathrm{hyper}) \, \pi(\Gamma \mid \mathrm{hyper})}.
+Sampling is ancestral and exact: hyperparameters from their prior,
+then independent pair flips.
+\item \code{"joint"}: the graph marginal is reweighted by the per-graph
+normalizer of the determinant-tilted precision prior,
+\eqn{q(\Gamma \mid \mathrm{hyper}) \propto Z(\Gamma) \,
+    \pi(\Gamma \mid \mathrm{hyper})}. Sampling runs the zero-data
+\eqn{(K, \Gamma)} chain of \code{\link{sample_ggm_prior}} and discards
+\eqn{K}; with a Beta-Bernoulli or Stochastic-Block prior the
+hyperparameter updates apply the normalizing-constant correction, so
+the first call for a model cell may build the correction table
+(cached across fits).
+}
+}
+\details{
+The optional conditioning arguments fix the edge-prior hyperparameters
+instead of sampling them: \code{theta} fixes the inclusion probability of
+a Bernoulli or Beta-Bernoulli prior, and \code{allocations} plus
+\code{block_probs} fix the block structure of a Stochastic-Block prior
+(reducing it to independent pair flips at the given block probabilities).
+
+Under \code{spec = "joint"} the tilted graph law depends on the precision
+prior through its normalizer, so \code{interaction_prior},
+\code{precision_scale_prior}, and \code{delta} are part of the graph law;
+they are ignored under \code{spec = "hierarchical"}.
+}
+\examples{
+# Hierarchical spec: the graph marginal is exactly the edge prior.
+g = sample_graph_prior(
+  p = 6, n_samples = 200,
+  edge_prior = bernoulli_prior(0.3), seed = 11
+)
+mean(g$edge_indicators) # about 0.3
+
+# Beta-Bernoulli: inclusion probabilities are sampled alongside.
+g = sample_graph_prior(
+  p = 6, n_samples = 200,
+  edge_prior = beta_bernoulli_prior(2, 4), seed = 11
+)
+mean(g$theta) # about 1/3
+
+\donttest{
+# Joint spec: the graph law carries the per-graph normalizer Z(Gamma).
+g = sample_graph_prior(
+  p = 6, n_samples = 500,
+  edge_prior = bernoulli_prior(0.3), spec = "joint",
+  interaction_prior = normal_prior(scale = 0.5),
+  precision_scale_prior = gamma_prior(shape = 1, rate = 2),
+  seed = 11, verbose = FALSE
+)
+mean(g$edge_indicators) # shifted away from 0.3 by the Z(Gamma) tilt
+}
+}
+\seealso{
+\code{\link{sample_ggm_prior}}, \code{\link{sample_sbm_prior}},
+\code{\link{bernoulli_prior}}, \code{\link{beta_bernoulli_prior}},
+\code{\link{sbm_prior}}, \code{\link{bgm}}
+}
diff --git a/man/sample_sbm_prior.Rd b/man/sample_sbm_prior.Rd
new file mode 100644
index 00000000..85a2cce8
--- /dev/null
+++ b/man/sample_sbm_prior.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sample_graph_prior.R
+\name{sample_sbm_prior}
+\alias{sample_sbm_prior}
+\title{Sample from the Stochastic-Block Edge-Prior Hyperprior}
+\usage{
+sample_sbm_prior(p, n_samples, edge_prior = sbm_prior(), seed = 1L)
+}
+\arguments{
+\item{p}{Integer. Number of nodes (\eqn{p \ge 2}).}
+
+\item{n_samples}{Integer. Number of prior draws.}
+
+\item{edge_prior}{An \code{\link{sbm_prior}()} object. Default
+\code{sbm_prior()}.}
+
+\item{seed}{Integer. Seed for the draw; the caller's RNG state is
+restored on exit.}
+}
+\value{
+A list with elements:
+\describe{
+\item{\code{allocations}}{Integer matrix (\code{n_samples x p}) of
+block labels.}
+\item{\code{pair_probability}}{Numeric matrix
+(\code{n_samples x p(p-1)/2}) of implied pair-inclusion
+probabilities, columns in row-major upper-triangle order.}
+\item{\code{num_blocks}}{Integer vector: number of occupied blocks
+per draw.}
+\item{\code{pair_names}}{Character vector labeling the pair columns
+as \code{"i-j"}.}
+\item{\code{p}}{The node count.}
+}
+}
+\description{
+Ancestral draws from the MFM-SBM hyperprior of
+\code{\link{sbm_prior}()}: a partition of the nodes into blocks
+(shifted-Poisson number of components, Dirichlet-weighted allocation) and
+Beta-distributed within- and between-block edge-inclusion probabilities.
+These are the hyperparameters that \code{\link{sample_graph_prior}} and
+\code{\link{bgm}} integrate over when the edge prior is a Stochastic-Block
+prior.
+}
+\examples{
+draws = sample_sbm_prior(p = 8, n_samples = 100, seed = 4)
+table(draws$num_blocks)
+range(draws$pair_probability)
+}
+\seealso{
+\code{\link{sbm_prior}}, \code{\link{sample_graph_prior}},
+\code{\link{bgm}}
+}
diff --git a/tests/testthat/helper-fixtures.R b/tests/testthat/helper-fixtures.R
index 82f9ba73..f35dff31 100644
--- a/tests/testthat/helper-fixtures.R
+++ b/tests/testthat/helper-fixtures.R
@@ -49,8 +49,10 @@
 # Ensure bgms package is loaded
 library(bgms)
 
-# Suppress informational messages during tests
-options(bgms.verbose = FALSE)
+# Advisory output is quieted for the test run in setup.R, not here. Helper files
+# are sourced by devtools::load_all() (setup files are not), so setting the
+# option here would leak bgms.verbose = FALSE into interactive development
+# sessions and silence fit-time messages and progress bars.
 
 # ------------------------------------------------------------------------------
 # 1. Session-Cached Model Fixtures
diff --git a/tests/testthat/test-correction-tables.R b/tests/testthat/test-correction-tables.R
index a9361da2..8813b46a 100644
--- a/tests/testthat/test-correction-tables.R
+++ b/tests/testthat/test-correction-tables.R
@@ -108,6 +108,54 @@ test_that("table build with cache round-trips and reuses the file", {
   expect_equal(tab1$cell$delta, 0.5 * log(4))
 })
 
+test_that("the progress bar renders the label, counts, and percentage", {
+  pb = new_correction_progress(120L, prefix = "Correction table")
+  mid = paste(capture.output(pb$update(70L)), collapse = "")
+  expect_match(mid, "Correction table:", fixed = TRUE)
+  expect_match(mid, "70/120", fixed = TRUE)
+  expect_match(mid, "58.3%", fixed = TRUE)
+  full = paste(capture.output(pb$update(120L)), collapse = "")
+  expect_match(full, "120/120 (100.0%)", fixed = TRUE)
+})
+
+test_that("the parallel sweep matches the serial sweep cell for cell", {
+  skip_on_cran()
+  skip_on_os("windows")
+
+  args = list(
+    p = 4, theta = c(0.2, 0.5, 0.8), delta = 0.5 * log(4),
+    interaction_prior = cauchy_prior(scale = 2.5),
+    precision_scale_prior = gamma_prior(shape = 1, eta = 1),
+    n_samples = 200L, n_warmup = 100L, n_seeds = 2L, update_method = "gibbs"
+  )
+  serial = do.call(sweep_prior_edge_density, c(args, cores = 1L))
+  parallel = do.call(sweep_prior_edge_density, c(args, cores = 2L))
+
+  expect_equal(parallel$edens_raw, serial$edens_raw)
+})
+
+test_that("the build announces itself once; a cache hit is silent", {
+  cache_dir = file.path(tempdir(), "bgms-ctable-msg-test")
+  unlink(cache_dir, recursive = TRUE)
+  old = options(bgms.correction_cache_dir = cache_dir)
+  on.exit(options(old), add = TRUE)
+
+  build_args = list(
+    p = 4, n_grid = 12L, n_samples = 100L, n_warmup = 100L, n_seeds = 1L,
+    update_method = "gibbs", verbose = TRUE
+  )
+  # capture.output silences the progress bar (stdout); messages pass through.
+  capture.output(
+    expect_message(
+      do.call(ggm_correction_table, build_args),
+      "Building the edge-selection prior correction table"
+    )
+  )
+  capture.output(
+    expect_no_message(do.call(ggm_correction_table, build_args))
+  )
+})
+
 test_that("gibbs and adaptive-metropolis sweeps agree on edge density", {
   skip_on_cran()
 
diff --git a/tests/testthat/test-mixed-correction.R b/tests/testthat/test-mixed-correction.R
index 00d0248e..c426418c 100644
--- a/tests/testthat/test-mixed-correction.R
+++ b/tests/testthat/test-mixed-correction.R
@@ -242,19 +242,28 @@ test_that("bgm() applies the correction to a mixed beta-bernoulli fit", {
   skip_on_cran()
   old = mixed_correction_cache()
   on.exit(options(old), add = TRUE)
+  # A fresh cache forces an actual build, the only case that announces
+  # itself; a cache hit is silent.
+  unlink(getOption("bgms.correction_cache_dir"), recursive = TRUE)
 
   d = mixed_smoke_data()
-  msgs = capture.output(
-    fit <- bgm(d$x,
-      variable_type = d$variable_type,
-      edge_prior = beta_bernoulli_prior(),
-      iter = 300, warmup = 300, chains = 1, cores = 1,
-      display_progress = "none", verbose = TRUE
+  # Capture stdout too, so an interactive test run does not leak the
+  # build progress bar; the announcement is a message (stderr).
+  msgs = NULL
+  capture.output(
+    msgs <- capture.output(
+      fit <- bgm(d$x,
+        variable_type = d$variable_type,
+        edge_prior = beta_bernoulli_prior(),
+        iter = 300, warmup = 300, chains = 1, cores = 1,
+        display_progress = "none", verbose = TRUE
+      ),
+      type = "message"
     ),
-    type = "message"
+    type = "output"
   )
 
-  expect_true(any(grepl("Edge-prior correction", msgs)))
+  expect_true(any(grepl("correction table", msgs)))
   expect_s3_class(fit, "bgms")
   th = fit$inclusion_parameter_samples
   expect_length(th, 1L)
@@ -294,7 +303,7 @@ test_that("bgm() skips the correction with one continuous variable", {
     type = "message"
   )
 
-  expect_false(any(grepl("Edge-prior correction", msgs)))
+  expect_false(any(grepl("correction table", msgs)))
   expect_s3_class(fit, "bgms")
 })
 
diff --git a/tests/testthat/test-sample-graph-prior.R b/tests/testthat/test-sample-graph-prior.R
new file mode 100644
index 00000000..08d377bc
--- /dev/null
+++ b/tests/testthat/test-sample-graph-prior.R
@@ -0,0 +1,152 @@
+# Tests for sample_graph_prior(): the hierarchical spec's ancestral graph
+# law, hyperparameter conditioning, and the joint spec's Z(Gamma)-tilted
+# negative control (mirroring test-hier-zratio-identity).
+
+test_that("hierarchical Bernoulli graph law matches its prior", {
+  g = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = bernoulli_prior(0.3), seed = 11
+  )
+  expect_identical(dim(g$edge_indicators), c(4000L, 15L))
+  expect_true(all(g$edge_indicators %in% 0:1))
+  expect_lt(abs(mean(g$edge_indicators) - 0.3), 0.02)
+  expect_null(g$theta)
+  expect_null(g$allocations)
+  expect_identical(g$spec, "hierarchical")
+  expect_identical(g$edge_prior, "Bernoulli")
+  expect_identical(g$pair_names[1:2], c("1-2", "1-3"))
+})
+
+test_that("hierarchical Beta-Bernoulli samples theta and coheres with it", {
+  a = 2
+  b = 4
+  g = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = beta_bernoulli_prior(a, b), seed = 5
+  )
+  expect_length(g$theta, 4000L)
+  expect_true(all(g$theta > 0 & g$theta < 1))
+  expect_lt(abs(mean(g$theta) - a / (a + b)), 0.02)
+  expect_lt(abs(var(g$theta) - a * b / ((a + b)^2 * (a + b + 1))), 0.01)
+  # The graph marginal integrates theta out to a/(a + b), and per draw the
+  # edge mean tracks the sampled theta.
+  expect_lt(abs(mean(g$edge_indicators) - a / (a + b)), 0.02)
+  expect_gt(cor(rowMeans(g$edge_indicators), g$theta), 0.7)
+})
+
+test_that("hierarchical SBM separates within- from between-block edges", {
+  g = sample_graph_prior(
+    p = 8, n_samples = 1000,
+    edge_prior = sbm_prior(
+      alpha = 6, beta = 1, alpha_between = 1, beta_between = 6
+    ),
+    seed = 7
+  )
+  expect_identical(dim(g$allocations), c(1000L, 8L))
+  expect_true(all(g$allocations >= 1L))
+
+  pairs = which(upper.tri(matrix(0, 8, 8)), arr.ind = TRUE)
+  pairs = pairs[order(pairs[, 1L], pairs[, 2L]), , drop = FALSE]
+  within = matrix(
+    g$allocations[, pairs[, 1L]] == g$allocations[, pairs[, 2L]],
+    nrow = 1000L
+  )
+  ind = g$edge_indicators
+  expect_gt(mean(ind[within]), mean(ind[!within]) + 0.3)
+})
+
+test_that("conditioning on theta fixes the graph law", {
+  g = sample_graph_prior(
+    p = 6, n_samples = 2000,
+    edge_prior = beta_bernoulli_prior(1, 1), theta = 0.8, seed = 3
+  )
+  expect_lt(abs(mean(g$edge_indicators) - 0.8), 0.02)
+  expect_null(g$theta)
+})
+
+test_that("conditioning on (allocations, block_probs) fixes the SBM law", {
+  z = c(1L, 1L, 1L, 2L, 2L, 2L)
+  bp = matrix(c(0.9, 0.05, 0.05, 0.9), 2, 2)
+  g = sample_graph_prior(
+    p = 6, n_samples = 2000,
+    edge_prior = sbm_prior(), allocations = z, block_probs = bp, seed = 9
+  )
+  pairs = which(upper.tri(matrix(0, 6, 6)), arr.ind = TRUE)
+  pairs = pairs[order(pairs[, 1L], pairs[, 2L]), , drop = FALSE]
+  within = z[pairs[, 1L]] == z[pairs[, 2L]]
+  expect_lt(abs(mean(g$edge_indicators[, within]) - 0.9), 0.03)
+  expect_lt(abs(mean(g$edge_indicators[, !within]) - 0.05), 0.03)
+  expect_null(g$allocations)
+})
+
+test_that("conditioning arguments are validated", {
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = sbm_prior(), theta = 0.5
+    ),
+    "allocations"
+  )
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = sbm_prior(), allocations = rep(1L, 5)
+    ),
+    "block_probs"
+  )
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = bernoulli_prior(0.5),
+      allocations = rep(1L, 5),
+      block_probs = matrix(0.5, 1, 1)
+    ),
+    "sbm_prior"
+  )
+  expect_error(
+    sample_graph_prior(p = 5, n_samples = 10, theta = 1.5),
+    "theta"
+  )
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = sbm_prior(),
+      allocations = c(1L, 1L, 2L, 2L, 3L),
+      block_probs = matrix(0.5, 2, 2)
+    ),
+    "exceed"
+  )
+})
+
+test_that("the same seed reproduces the draw and restores the RNG state", {
+  set.seed(123)
+  before = runif(1)
+  set.seed(123)
+  g1 = sample_graph_prior(p = 6, n_samples = 50, seed = 42)
+  after = runif(1)
+  g2 = sample_graph_prior(p = 6, n_samples = 50, seed = 42)
+  expect_identical(g1$edge_indicators, g2$edge_indicators)
+  expect_identical(before, after)
+})
+
+test_that("the joint spec produces the Z(Gamma)-tilted law, not the prior", {
+  skip_on_cran()
+  # Mirrors the negative control in test-hier-zratio-identity: at
+  # tau^2 = 2 the joint marginal separates cleanly from Bernoulli(0.3)
+  # (about 0.22), while the hierarchical law sits at 0.3.
+  g = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = bernoulli_prior(0.3), spec = "joint",
+    interaction_prior = normal_prior(scale = 0.5),
+    precision_scale_prior = gamma_prior(shape = 1, rate = 2),
+    n_warmup = 1500, seed = 7, verbose = FALSE
+  )
+  expect_identical(dim(g$edge_indicators), c(4000L, 15L))
+  expect_gt(abs(mean(g$edge_indicators) - 0.3), 0.05)
+
+  h = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = bernoulli_prior(0.3), seed = 7
+  )
+  expect_lt(abs(mean(h$edge_indicators) - 0.3), 0.02)
+})
diff --git a/tests/testthat/test-sample-sbm-prior.R b/tests/testthat/test-sample-sbm-prior.R
new file mode 100644
index 00000000..0ac2124f
--- /dev/null
+++ b/tests/testthat/test-sample-sbm-prior.R
@@ -0,0 +1,59 @@
+# Tests for sample_sbm_prior(): ancestral draws from the MFM-SBM
+# edge-prior hyperprior.
+
+test_that("sample_sbm_prior returns coherent draws", {
+  d = sample_sbm_prior(p = 8, n_samples = 500, seed = 4)
+  expect_identical(dim(d$allocations), c(500L, 8L))
+  expect_identical(dim(d$pair_probability), c(500L, 28L))
+  expect_length(d$num_blocks, 500L)
+  expect_true(all(d$allocations >= 1L))
+  expect_true(all(d$pair_probability > 0 & d$pair_probability < 1))
+  expect_true(all(d$num_blocks >= 1L & d$num_blocks <= 8L))
+  # num_blocks counts the occupied blocks of each allocation.
+  expect_identical(
+    d$num_blocks,
+    apply(d$allocations, 1, function(z) length(unique(z)))
+  )
+})
+
+test_that("pair probabilities follow the block structure", {
+  # Dense within (Beta(8, 1)), sparse between (Beta(1, 8)).
+  d = sample_sbm_prior(
+    p = 8, n_samples = 500,
+    edge_prior = sbm_prior(
+      alpha = 8, beta = 1, alpha_between = 1, beta_between = 8
+    ),
+    seed = 2
+  )
+  pairs = which(upper.tri(matrix(0, 8, 8)), arr.ind = TRUE)
+  pairs = pairs[order(pairs[, 1L], pairs[, 2L]), , drop = FALSE]
+  within = matrix(
+    d$allocations[, pairs[, 1L]] == d$allocations[, pairs[, 2L]],
+    nrow = 500L
+  )
+  expect_gt(mean(d$pair_probability[within]), 0.8)
+  expect_lt(mean(d$pair_probability[!within]), 0.2)
+})
+
+test_that("lambda shifts the number of blocks", {
+  few = sample_sbm_prior(
+    p = 10, n_samples = 400,
+    edge_prior = sbm_prior(lambda = 0.1), seed = 6
+  )
+  many = sample_sbm_prior(
+    p = 10, n_samples = 400,
+    edge_prior = sbm_prior(lambda = 8), seed = 6
+  )
+  expect_lt(mean(few$num_blocks), mean(many$num_blocks))
+})
+
+test_that("input validation and reproducibility", {
+  expect_error(
+    sample_sbm_prior(p = 8, n_samples = 10, edge_prior = bernoulli_prior()),
+    "sbm_prior"
+  )
+  expect_error(sample_sbm_prior(p = 1, n_samples = 10), "'p'")
+  d1 = sample_sbm_prior(p = 6, n_samples = 20, seed = 9)
+  d2 = sample_sbm_prior(p = 6, n_samples = 20, seed = 9)
+  expect_identical(d1, d2)
+})