classDiagram
class SamplerBase {
<<abstract>>
+step(model, iteration) StepResult
+initialize(model)
+has_nuts_diagnostics() bool
}
class MetropolisSampler {
element-wise MH
eager adaptation init
}
class NUTSSampler {
owns NUTSAdaptationController
adaptive tree depth
}
class GibbsSampler {
conjugate row sweep
exact draw, accept_prob 1
}
SamplerBase <|-- MetropolisSampler
SamplerBase <|-- NUTSSampler
SamplerBase <|-- GibbsSampler
Sampler Hierarchy
The sampler classes in src/mcmc/samplers/ are header-only wrappers that own adaptation state and delegate to the stateless algorithm functions in src/mcmc/algorithms/ (or to model-specific update code). The chain runner interacts with samplers through the SamplerBase interface.
Class hierarchy
One sampler instance drives the whole model: the chain runner issues a single sampler->step(model, iter) per iteration, and under NUTS all parameter blocks — including the mixed MRF’s five blocks (thresholds, means, and the three interaction blocks) — are updated jointly in one trajectory (see Model Classes). The models’ own element-wise Metropolis code runs only when the user selects update_method = "adaptive-metropolis", and the GGM’s conjugate sweep only under update_method = "gibbs".
SamplerBase interface
class SamplerBase {
virtual StepResult step(BaseModel& model, int iteration) = 0;
virtual void initialize(BaseModel& model) {}
virtual bool has_nuts_diagnostics() const { return false; }
virtual void set_warm_step_size(double eps) {}
virtual double get_final_step_size() const; // NaN if none
virtual void set_warm_inv_mass(const arma::vec& inv_mass) {}
virtual arma::vec get_final_inv_mass() const; // empty if none
};step() performs one MCMC iteration and returns a StepResult containing the accepted parameter vector, acceptance probability, and optional NUTS diagnostics.
The four warm-start methods carry adapted state between fits: a refit injects the previous fit’s step size and diagonal metric before the loop, and the chain runner reads the adapted values back out afterwards (see Parallel Chains). Non-gradient samplers ignore the setters and return the NaN and empty sentinels.
initialize() is called once by the chain runner before the MCMC loop. For NUTS it runs the step-size heuristic that finds a reasonable initial step size by iteratively doubling or halving until the acceptance probability is near the target (default 0.8). For Metropolis it eagerly sets up the model’s proposal-SD adaptation controllers (init_metropolis_adaptation). For Gibbs it switches the model’s between-model step to the conjugate edge proposal (set_conjugate_edge_proposal(true)).
MetropolisSampler
Each call to step() invokes model.do_one_metropolis_step(iteration), which performs a full sweep of element-wise proposals, and returns the model’s mean acceptance probability over the sweep (last_metropolis_mean_accept_prob()). That per-iteration value is stored as the am_accept_prob__ trace during sampling and summarized in fit$am_diag (see Metropolis Updates). Edge selection is handled by the chain runner before the sampler step, not by MetropolisSampler itself.
Proposal standard deviations are tuned during warmup via Robbins-Monro updates (see Warmup Schedule).
NUTSSampler
Uses the No-U-Turn criterion to adapt trajectory length automatically. The maximum tree depth is set from config.max_tree_depth (default 10). Returns NUTS diagnostics (tree depth, divergence flag, energy, mean acceptance probability) with each step. See NUTS Algorithm for the tree-building details.
NUTSSampler owns a NUTSAdaptationController that manages:
- Dual averaging for step-size adaptation
- Diagonal mass matrix accumulation via Welford’s online algorithm (a numerically stable running-variance estimate), estimated on the full (zero-padded) parameter layout (see NUTS — Mass matrix)
Each call to step():
- Runs one NUTS trajectory in the active (theta-space) parameterization — graph constraints are enforced by the models’ null-space coordinates, so there is no separate constrained path
- Passes the acceptance probability to the adaptation controller
- If the mass matrix was just updated (end of a warmup window), reinitializes the step size with the heuristic; at the stage-3c boundary (edge selection activating), restarts dual averaging
After warmup, the adaptation controller freezes and returns the smoothed step size from dual averaging.
GibbsSampler
A thin wrapper for the GGM’s conjugate updates (update_method = "gibbs"; see GGM Internals for the row-block draw). Each call to step() invokes model.do_one_gibbs_step(iteration) — an exact conjugate row sweep — and reports accept_prob = 1 (an exact draw always “accepts”). With edge selection, the between-model step uses the full-conditional edge birth/death proposal, which needs no tuning. Two gates keep this sampler on models whose full conditionals are conjugate: the R layer rejects update_method = "gibbs" on anything but all-continuous data, and the GGM entry point rejects an unsupported prior family (row_block_gibbs_eligible() requires a Normal or Cauchy slab with a Gamma diagonal).
Warmup staging under Gibbs is minimal: a short full-model settle window, then selection-active warmup for the remainder (see Warmup Schedule).
Factory dispatch
Dispatch happens in two steps in chain_runner.cpp. The R layer passes the update_method string through unchanged; resolve_sampler_spec() maps it to a typed SamplerSpec (any other string is an error):
SamplerSpec resolve_sampler_spec(const std::string& sampler_type) {
if (sampler_type == "nuts") {
return SamplerSpec{SamplerKind::NUTS, /*learn_sd=*/true,
/*nuts_diag=*/true, /*am_diag=*/false};
} else if (sampler_type == "adaptive-metropolis") {
return SamplerSpec{SamplerKind::AdaptiveMetropolis, /*learn_sd=*/false,
/*nuts_diag=*/false, /*am_diag=*/true};
} else if (sampler_type == "gibbs") {
return SamplerSpec{SamplerKind::Gibbs, /*learn_sd=*/false,
/*nuts_diag=*/false, /*am_diag=*/false};
} else {
// std::runtime_error rather than Rcpp::stop: this runs on worker
// threads, where constructing an Rcpp exception is not safe.
throw std::runtime_error("Unknown sampler_type: '" + sampler_type + "'");
}
}create_sampler() then switches on the enum:
unique_ptr<SamplerBase> create_sampler(SamplerKind kind,
const SamplerConfig& config,
WarmupSchedule& schedule) {
switch (kind) {
case SamplerKind::NUTS:
return make_unique<NUTSSampler>(config, schedule);
case SamplerKind::AdaptiveMetropolis:
return make_unique<MetropolisSampler>(config, schedule);
case SamplerKind::Gibbs:
return make_unique<GibbsSampler>(config, schedule);
}
throw std::runtime_error("Unhandled SamplerKind");
}The SamplerSpec records which diagnostics the sampler produces (nuts_diag for NUTS, am_diag for adaptive Metropolis — the per-iteration mean-acceptance trace) and whether proposal-SD learning applies (learn_sd), which the warmup schedule consumes. The spec also determines the warmup staging host: the schedule is built with select_during_warmup = (kind == Gibbs), which switches it to the Gibbs settle-then-select staging (see Warmup Schedule).
All three model types (GGM, OMRF, mixed MRF) default to NUTS. The user can select Metropolis via update_method = "adaptive-metropolis" for any model, and the conjugate Gibbs sweep via update_method = "gibbs" for all-continuous (GGM) data.
SamplerConfig
SamplerConfig (mcmc/execution/sampler_config.h) is the second argument to create_sampler() above: one plain struct carrying the run settings the chain runner and the samplers need. Each of the three bgm() entry points (sample_ggm.cpp, sample_omrf.cpp, sample_mixed.cpp) fills one from its R arguments and passes it by const reference from there on, so no sampler reads the R list and nothing downstream can disagree about a setting. bgmCompare() does not use it at all: it runs its own sampler in models/bgmCompare/bgmCompare_sampler.cpp and never constructs a SamplerConfig.
The struct holds eleven fields in four groups: the sampler selection (sampler_type, the same string resolve_sampler_spec() consumes); the run length (no_iter, no_warmup); the gradient-sampler tuning (max_tree_depth, initial_step_size, target_acceptance, learn_mass_matrix); and the per-run switches (edge_selection, na_impute, zratio_gauge_sweeps, seed).
Only NUTSSampler reads any of it. MetropolisSampler and GibbsSampler both take the config and discard it with (void)config: Metropolis gets its adaptation from the model and the schedule, and the Gibbs draw is exact, so neither has a tuning parameter to receive.
Two fields are easy to misread. seed is the base seed, not the seed a chain uses: the runner gives chain c the seed config.seed + c (see Parallel Chains). And zratio_gauge_sweeps is the number of post-sampling assessment sweeps the in-chain trust gauge runs, with 0 switching it off; the reference draws per pair and the per-sweep referenced-pair cap are not settings here, but fixed design parameters owned by ZRatioGauge itself. It is filled by the GGM and mixed entry points only, since an ordinal model has no continuous block to gauge.
Every field carries a default in the struct declaration, and those declared values are not the defaults a user sees. sampler_type declares "adaptive-metropolis" while the R layer resolves NUTS for every model, and each entry point overwrites it before the struct is used. initial_step_size is the one field no production entry point sets, so NUTS always starts from the declared 0.1 and adapts from there; it is overridden only in the gradient test interface. Read the declared values as fallbacks, and R Scaffolding for the defaults that actually reach a fit.
Step-size and mass-matrix adaptation
The full adaptation schedule for the NUTSAdaptationController introduced above — stage boundaries, doubling windows, blending formula, and dual-averaging parameters — is documented in Warmup Schedule.
MetropolisSampler uses a separate MetropolisAdaptationController that tunes proposal standard deviations via Robbins–Monro updates targeting an acceptance rate of 0.44, the componentwise random-walk optimum. Note that the same 0.44 target also governs the models’ edge-move proposals during NUTS runs; the user’s target_accept argument applies to the Metropolis tuner only under update_method = "adaptive-metropolis" (under NUTS it sets the step-size adaptation target instead, and under Gibbs it is unused — exact draws have no acceptance target).
GibbsSampler performs no adaptation at all.