Model Classes
This page has not yet passed technical or readability review.
All graphical models in bgms inherit from a single abstract base class, BaseModel, defined in src/models/base_model.h. The chain runner interacts with models exclusively through this interface, so adding a new model type requires implementing the relevant virtual methods without modifying the MCMC engine.
BaseModel interface
BaseModel declares over 30 virtual methods. The most important groups are:
Capability queries
virtual bool has_edge_selection() const;
virtual bool has_missing_data() const;
virtual bool gauge_available() const;The interface does not ask models what sampler they support: the sampler choice is configuration-driven. The R layer validates the combination of update_method, variable types, and priors (validate_sampler()), and the chain runner dispatches on the resulting SamplerSpec (see Sampler Hierarchy). What remains are queries about model features: whether spike-and-slab edge selection applies, whether missing-data imputation is active, and whether a Z-ratio trust gauge is attached (hierarchical graph prior only).
Sampling methods
virtual pair<double, vec> logp_and_gradient(const vec& parameters);
virtual void do_one_metropolis_step(int iteration) = 0;
virtual void do_one_gibbs_step(int iteration);
virtual void set_conjugate_edge_proposal(bool enable);
virtual double last_metropolis_mean_accept_prob() const;
virtual void set_metropolis_target_accept(double target);Every model must implement do_one_metropolis_step(), a full sweep of element-wise proposals. Models that support NUTS override logp_and_gradient(); the joint function is more efficient than separate calls because the log-posterior and gradient share intermediate computations. do_one_gibbs_step() defaults to a throw: only models with an exact conjugate full-conditional sweep override it (currently the GGM row-block Gibbs), and set_conjugate_edge_proposal() switches the GGM’s between-model step from the random-walk proposal to the full-conditional edge birth/death used by the Gibbs sampler.
last_metropolis_mean_accept_prob() reports the mean acceptance probability of the most recent Metropolis sweep (NaN for models without one); the adaptive-Metropolis sampler stores it per iteration as the am_accept_prob__ trace. set_metropolis_target_accept() receives the resolved acceptance target before the MCMC loop — the user’s target_accept under update_method = "adaptive-metropolis", or the fixed 0.44 edge-move target under NUTS.
Lifecycle hooks
virtual void prepare_iteration();
virtual void collect_chain_diagnostics(ChainResult& chain_result) const;prepare_iteration() runs at the start of every iteration (e.g., to shuffle the edge-update order). collect_chain_diagnostics() runs once at the end of the chain and copies run-level diagnostic state (e.g., the Z-ratio engine’s counters and calibration anchors) into the ChainResult.
Trust-gauge hooks
virtual bool gauge_available() const;
virtual void set_gauge_active(bool on, int n_draws, int cap);
virtual void gauge_begin_sweep();
virtual void gauge_end_sweep();Under the hierarchical graph prior, the sampler decides edge moves with a fast normalizer-ratio approximation, and an optional trust gauge audits those decisions after sampling (see summarize_zratio_gauge()). While the gauge is active, update_edge_indicators() re-decides each non-trivial edge move (one with a non-empty mediating block) and also evaluates the exact block-local reference for it, so the fast decision and the exact one can be compared; the begin/end hooks bracket each assessment sweep. All four methods are no-ops for models without a Z-ratio engine.
Parameter access
virtual vec get_vectorized_parameters() const = 0;
virtual void set_vectorized_parameters(const vec& params);
virtual vec get_full_vectorized_parameters() const = 0;
virtual vec get_storage_vectorized_parameters() const;
virtual size_t parameter_dimension() const = 0;
virtual size_t full_parameter_dimension() const;
virtual size_t storage_dimension() const;Three vector layouts serve three consumers:
- Active (
get_vectorized_parameters(),parameter_dimension()) — the free coordinates the sampler updates. Under edge selection this dimension changes when the graph changes. - Full (
get_full_vectorized_parameters(),full_parameter_dimension()) — the zero-padded layout in which every possible coordinate keeps a fixed slot. The NUTS adaptation controller estimates the mass matrix on this layout so entries survive active-set changes (see NUTS — Mass matrix). - Storage (
get_storage_vectorized_parameters(),storage_dimension()) — the fixed-size layout saved to output each retained iteration, chosen to match what the R layer expects.
For the GGM the three are: active \(p + |E|\) (log-diagonals plus null-space coordinates of included edges), full \(p + p(p-1)/2\), and storage \(p(p+1)/2\) — the upper triangle of the precision matrix \(\boldsymbol{\Theta}\) itself, not the theta coordinates.
Edge selection
virtual void update_edge_indicators() = 0;
virtual void set_edge_selection_active(bool active);
virtual const arma::imat& get_edge_indicators() const = 0;
virtual arma::mat& get_inclusion_probability() = 0;When edge selection is active, update_edge_indicators() performs a full scan over all edges, proposing add/delete moves (Gottardo & Raftery, 2008; van den Bergh et al., 2026). Excluded edges are represented as parameters fixed at zero. The chain runner calls this before the sampler step on iterations where the warmup schedule enables selection, and passes the indicator matrix and inclusion probabilities to the edge-prior update afterwards.
Adaptation support
virtual void init_metropolis_adaptation(const WarmupSchedule& schedule);
virtual void tune_proposal_sd(int iteration, const WarmupSchedule& schedule);Models with element-wise Metropolis proposals implement these to tune proposal standard deviations during warmup stage 3b, using the Robbins-Monro weight supplied by the schedule’s rm_weight_for_proposal_sd() (see Warmup Schedule).
Mass matrix access
virtual void set_inv_mass(const arma::vec& inv_mass);
virtual const arma::vec& get_inv_mass() const;
virtual arma::vec get_active_inv_mass() const;The stored inverse-mass diagonal lives on the full layout; get_active_inv_mass() gathers the entries for the currently active coordinates, which is what the integrator receives.
Infrastructure
virtual void set_seed(int seed) = 0;
virtual unique_ptr<BaseModel> clone() const = 0;
virtual SafeRNG& get_rng() = 0;clone() produces an independent deep copy for parallel chain execution. Each clone gets its own RNG state via set_seed().
Graph constraints
There is no constraint-projection interface. Models with continuous variables enforce excluded edges and positive definiteness through the constrained Cholesky parameterization: the active parameter vector contains only free coordinates, every value of which maps to a valid precision matrix. The sampler never sees a constraint — NUTS runs its standard unconstrained integrator on the active vector for every model type.
GGMModel
Defined in src/models/ggm/ggm_model.h and ggm_model.cpp.
Parameterization. The GGM is parameterized by the precision matrix \(\boldsymbol{\Theta}\) (inverse covariance). The model maintains an upper-triangular Cholesky factor \(\boldsymbol{\Phi}\) such that \(\boldsymbol{\Phi}^\top\boldsymbol{\Phi} = \boldsymbol{\Theta}\) (see GGM Internals and Constrained Cholesky Parameterization, which use the same convention).
Sampling. The GGM supports three update methods. Under NUTS, the model works in theta space: log-scale diagonals for positivity and per-column null-space coordinates for the off-diagonals, so excluded edges stay exactly zero without projection; logp_and_gradient() computes the gradient in this parameterization. Under element-wise Metropolis-Hastings, each off-diagonal element is proposed from a Gaussian centered on the current value, and the Cholesky factor is updated via a rank-1 Givens rotation (or hyperbolic rotation for downdates), avoiding a full \(O(p^3)\) decomposition. Under Gibbs (update_method = "gibbs"), do_one_gibbs_step() performs an exact conjugate row-block sweep, and the between-model step uses the full-conditional edge birth/death proposal (see GGM Internals).
Sufficient statistics. For complete data, the GGM only needs the sufficient statistic \(\mathbf{X}^\top\mathbf{X}\), not the raw data matrix. When missing data imputation is active, the sufficient statistics are recomputed after each imputation step.
Spike-and-slab. Edge indicators \(\gamma_{ij} \in \{0, 1\}\) control whether edge \((i,j)\) is included. When \(\gamma_{ij} = 0\), the corresponding precision element is exactly zero by construction of the null-space coordinates.
OMRFModel
Defined in src/models/omrf/omrf_model.h and omrf_model.cpp.
Parameterization. The ordinal MRF has two parameter types: main effects (category thresholds) and pairwise interactions. For ordinal variables, there is one threshold per category. For Blume-Capel variables, there are two parameters: a linear effect \(\alpha\) and a quadratic effect \(\beta\), with the category score \(\mu_c = \alpha(c - b) + \beta(c - b)^2\) where \(b\) is the reference category.
Gradient. The OMRF supports logp_and_gradient(), enabling NUTS sampling. The gradient involves category probabilities computed by the variable helpers, which use a FAST/SAFE dual-path strategy for numerical stability.
Residual matrix. The model maintains a precomputed residual matrix \(r_{ni} = 2\sum_{j \neq i} x_{nj} \omega_{ij}\) that caches the contribution of pairwise interactions to each variable’s conditional distribution. This matrix is updated incrementally when a single interaction changes.
Edge selection. Edge indicators are updated via add/delete moves (Gottardo & Raftery, 2008; van den Bergh et al., 2026), independent of the parameter sampler. Because all parameters are unconstrained reals, the OMRF handles an excluded edge by simply dropping it from the active parameter vector. The OMRF also provides an element-wise adaptive Metropolis fallback for parameter updates (used during warmup proposal-SD tuning and as a user-selectable alternative to NUTS).
MixedMRFModel
Defined in src/models/mixed/mixed_mrf_model.h and four .cpp files.
Parameterization. The mixed MRF handles \(p\) discrete and \(q\) continuous variables jointly. The log-density has the form:
\[ \log f(\mathbf{x}, \mathbf{y}) = \sum_i \mu_i(x_i) + \mathbf{x}^\top \boldsymbol{\Omega}_{xx} \mathbf{x} + \mathbf{y}^\top \boldsymbol{\Omega}_{yy} \mathbf{y} + 2\mathbf{x}^\top \boldsymbol{\Omega}_{xy} \mathbf{y} \]
where \(\boldsymbol{\Omega}_{xx}\), \(\boldsymbol{\Omega}_{yy}\), and \(\boldsymbol{\Omega}_{xy}\) are the discrete-discrete, continuous-continuous, and cross-type pairwise interaction matrices (see the User’s Guide for their statistical interpretation).
Sampling. All five parameter blocks are updated jointly by NUTS. The continuous precision block uses the same constrained Cholesky parameterization as the GGM — log-scale diagonals plus null-space coordinates for the off-diagonals — so positive definiteness and excluded continuous-continuous edges are handled by construction; the gradient includes the Jacobian of the theta-to-precision mapping. Discrete-discrete and cross edges are plain coordinates that are dropped from the active vector when excluded. The model also provides a component-wise adaptive Metropolis fallback (do_one_metropolis_step()) used when the sampler is set to "adaptive-metropolis".
Three dimension concepts.
parameter_dimension()— Active parameters: thresholds, means, included discrete and cross edges, plus the continuous block’s \(q + |E_{yy}|\) theta coordinatesfull_parameter_dimension()— The zero-padded layout with a slot for every possible parameterstorage_dimension()— The fixed output layout (raw precision entries for the continuous block)
See Mixed MRF Internals for the block layout and the discrete pseudolikelihood.