sample_ggm_prior()

Draws precision matrices \(K\) from the prior of a Gaussian graphical model. The likelihood is omitted (\(n = 0\), \(S = 0\)), so the chain targets the prior alone.

Description

Three prior specifications are supported via the spec argument:

"conditional" (default): fix a graph \(\Gamma\) and sample \(K \mid \Gamma\) using the same theta-space NUTS sampler that drives bgm() for continuous data. Off-diagonals at excluded positions are held exactly at zero throughout the chain.
"joint": sample \((K, \Gamma)\) jointly from the un-normalised joint prior, using the adaptive-Metropolis (or Gibbs) chain from bgm() with edge selection on and the likelihood off. The marginal on \(\Gamma\) is then the edge prior times a graph-dependent normalizing constant, \(\pi(\Gamma) \cdot Z(\Gamma)\) — the joint specification, not the hierarchical one. Useful for simulation-based calibration of bgm()’s default sampler.
"hierarchical": sample \((K, \Gamma)\) from the hierarchical specification \(p(\Gamma)\,p(K \mid \Gamma)\) with \(p(K \mid \Gamma)\) normalized per graph, so the marginal on \(\Gamma\) is exactly the edge prior \(\pi(\Gamma)\). The per-graph normalizer ratio in each between-edge move is evaluated by the deterministic local Z-ratio approximation. Requires normal_prior() or cauchy_prior() interactions.

Priors are specified on the partial-association scale \(K_{yy} = -K/2\): interaction_prior acts on \(K_{yy,ij} = -K_{ij}/2\), so a normal_prior(scale = s) on the association scale is equivalent to a \(\mathrm{Normal}(0, 2s)\) prior on \(K_{ij}\) itself. Likewise precision_scale_prior acts on \(K_{ii}/2\): the default exponential_prior(eta = 1) with the default normal_prior(scale = 1) interaction prior resolves to \(K_{ii}/2 \sim \mathrm{Exp}(1)\) and therefore \(K_{ii} \sim \mathrm{Exp}(1/2)\) (mean \(2\)). The same convention is used by bgm() and the continuous block of the mixed-MRF model, so a prior argument passed here means the same distribution it would mean there. Output samples are reported as entries of \(K\); convert with \(K_{yy} = -K/2\) if you want them on the partial-association scale.

Usage

sample_ggm_prior(
  p,
  n_samples,
  n_warmup              = 2000,
  interaction_prior     = normal_prior(scale = 1),
  precision_scale_prior = exponential_prior(eta = 1),
  step_size             = 0.1,
  max_depth             = 10L,
  seed                  = 1L,
  verbose               = TRUE,
  edge_indicators       = NULL,
  delta                 = NULL,
  spec                  = c("conditional", "joint", "hierarchical"),
  edge_inclusion_prob   = 0.5,
  update_method         = c("adaptive-metropolis", "gibbs"),
  edge_prior            = NULL,
  apply_correction      = TRUE,
  zratio_diagnostics    = TRUE
)

Arguments

Argument	Description
`p`	Integer. Dimension of the precision matrix (\(p \ge 2\)).
`n_samples`	Integer. Number of post-warmup draws to keep.
`n_warmup`	Integer. NUTS warmup iterations. Default: `2000`.
`interaction_prior`	A `bgms_parameter_prior` for the partial-association off-diagonals \(K_{yy,ij} = -K_{ij}/2\). Allowed: `normal_prior()`, `cauchy_prior()`. `beta_prime_prior()` is not supported here. Default: `normal_prior(scale = 1)`, matching `bgm()`.
`precision_scale_prior`	A `bgms_scale_prior` for \(K_{ii}/2\). Allowed: `exponential_prior()`, `gamma_prior()`. Both accept the rate in the raw frame (`rate`) or the standardized frame (`eta`; the raw rate is derived as `eta / s` for interaction-prior scale `s`). Default: `exponential_prior(eta = 1)`.
`step_size`	Positive numeric. Initial NUTS step size used to seed dual-averaging adaptation (see NUTS). Default: `0.1`. Used only for `spec = "conditional"` (NUTS path); ignored otherwise.
`max_depth`	Integer. Maximum NUTS tree depth. Default: `10`. Used only for `spec = "conditional"`.
`seed`	Integer. RNG seed for the chain. Default: `1`.
`verbose`	Logical. If `TRUE` (default), print a progress bar.
`edge_indicators`	Optional integer \(p \times p\) matrix with `1` = edge included, `0` = excluded. Must be symmetric with `1`s on the diagonal. Default: full graph (all edges included). Used only for `spec = "conditional"` (the chain samples \(K \mid \Gamma\)); ignored when the chain samples \(\Gamma\).
`delta`	Non-negative numeric, or `NULL` for the dimension-adaptive default. Determinant-tilt exponent: multiplies the prior by \(\|K\|^{\delta}\), softly repelling the chain from the boundary of the positive-definite cone. `NULL` (default) auto-resolves to \(0.5 \log(p)\). Pass `delta = 0` for the untilted prior.
`spec`	One of `"conditional"` (default: sample \(K \mid \Gamma\) at fixed \(\Gamma\)), `"joint"` (sample \((K, \Gamma)\) jointly from the un-normalised joint prior), or `"hierarchical"` (sample \((K, \Gamma)\) from the per-graph normalized specification via the Z-ratio approximation).
`edge_inclusion_prob`	Probability in \((0, 1)\) for the Bernoulli edge prior used when `spec = "joint"`. Default: `0.5`. Ignored when `spec = "conditional"`.
`update_method`	One of `"adaptive-metropolis"` (default) or `"gibbs"`. Sampler driving the `spec = "joint"` chain; the Gibbs chain uses the conjugate row and edge updates and needs no proposal tuning. Ignored when `spec = "conditional"` (NUTS).
`edge_prior`	A `bgms_indicator_prior` object (`bernoulli_prior()`, `beta_bernoulli_prior()`, or `sbm_prior()`), or `NULL` (default) for a Bernoulli prior with probability `edge_inclusion_prob`. Only for `spec = "joint"`.
`apply_correction`	Logical. For the hierarchical edge priors (`beta_bernoulli_prior()`, `sbm_prior()`), apply the normalizing-constant correction to the hyperparameter updates (default `TRUE`; the correction table is built from the tilted prior sampler and cached across calls). With `FALSE` the plain conjugate updates are used, whose hyperparameter marginals do not match the hyperpriors under the determinant tilt.
`zratio_diagnostics`	Logical (default `TRUE`). Only for `spec = "hierarchical"`: run the trust gauge (`summarize_zratio_gauge()`) on the returned chain and attach the result; detected issues are printed when `verbose`. The gauge redoes a subset of the chain’s edge decisions with the exact calculation and reports two alarms: how often the decision outcome differs (`flip_rate`), and the projected distortion of the mean inclusion probability from the measured error under the edge prior’s feedback (`harm_pred`). Evidence-free sampling is the regime where the second alarm matters: a small consistent error can shift the graph marginal without flipping individual decisions.

Value

A list with components:

Component	Description
`K_offdiag`	Numeric matrix `n_samples` \(\times\) \(p(p-1)/2\) of upper-triangle off-diagonal entries of \(K\) for each draw, in row-major order (the upper triangle traversed by row): \((K_{12}, K_{13}, \ldots, K_{1p}, K_{23}, \ldots, K_{2p}, K_{34}, \ldots)\). Under `spec = "conditional"`, excluded edges are returned as `0`; under `spec = "joint"`, off-diagonals at excluded edges are sampled at `0` per the inclusion indicator.
`K_diag`	Numeric matrix `n_samples` \(\times\) \(p\) of diagonal entries \(K_{11}, \ldots, K_{pp}\).
`offdiag_names`	Character vector of length \(p(p-1)/2\) naming the columns of `K_offdiag` (e.g. `"K_1_2"`).
`diag_names`	Character vector of length \(p\) naming the columns of `K_diag`.
`edge_indicators`	Under `spec = "conditional"`, the \(p \times p\) integer matrix of fixed inclusion indicators used (full graph if not supplied). Under `spec = "joint"`, an `n_samples` \(\times\) \(p(p-1)/2\) integer matrix of sampled \(\Gamma_{ij}\) indicators (column order matches `K_offdiag`).
`theta`	Only with `beta_bernoulli_prior()`: numeric vector of length `n_samples` with the sampled inclusion probability.
`allocations`	Only with `sbm_prior()`: integer matrix (`n_samples` \(\times\) \(p\)) of sampled cluster allocations (1-based).
`zratio_diagnostics`	Only with `spec = "hierarchical"` and `zratio_diagnostics = TRUE`: the trust-gauge summary from `summarize_zratio_gauge()`.

Details

When spec = "joint", the chain is initialized from an ancestral draw of the edge prior (hyperparameters from their prior, then indicators given the hyperparameters), keyed to seed. Under a hierarchical edge prior the inclusion parameter and the graph density are coupled, and a full-graph start could otherwise pin both near 1 for a large number of sweeps in zero-evidence chains.

Examples

# Default Normal(0, 1) off-diagonal, Exponential(1) diagonal, p = 4.
draws = sample_ggm_prior(
  p = 4, n_samples = 200, n_warmup = 200,
  verbose = FALSE
)
dim(draws$K_offdiag) # 200 x 6
colnames(draws$K_offdiag) = draws$offdiag_names
head(draws$K_offdiag)

# Sparser graph: drop the (1, 4) edge.
E = matrix(1L, 4, 4)
E[1, 4] = E[4, 1] = 0L
draws = sample_ggm_prior(
  p = 4, n_samples = 200, n_warmup = 200,
  edge_indicators = E, verbose = FALSE
)
colnames(draws$K_offdiag) = draws$offdiag_names
all(draws$K_offdiag[, "K_1_4"] == 0) # TRUE