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/. The chain runner interacts with samplers through the SamplerBase interface.

Class hierarchy

classDiagram
    class SamplerBase {
        <<abstract>>
        +step(model, iteration) StepResult
        +initialize(model)
        +has_nuts_diagnostics() bool
    }
    class MetropolisSampler {
        element-wise MH
    }
    class NUTSSampler {
        owns NUTSAdaptationController
        adaptive tree depth
        -do_unconstrained_step()
        -do_constrained_step()
    }
    SamplerBase <|-- MetropolisSampler
    SamplerBase <|-- NUTSSampler

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".

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; }
};

step() performs one MCMC iteration and returns a StepResult containing the accepted parameter vector, acceptance probability, and optional NUTS diagnostics.

initialize() is called once before the first iteration. For gradient-based samplers, 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).

MetropolisSampler

The simplest sampler. Each call to step() invokes model.do_one_metropolis_step(iteration), which performs a full sweep of element-wise proposals, and returns the result. 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) 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)

Each call to step():

Dispatches inline to do_unconstrained_step() or do_constrained_step() based on model.has_constraints()
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

After warmup, the adaptation controller freezes and returns the smoothed step size from dual averaging.

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, ...};
  } else if (sampler_type == "adaptive-metropolis") {
    return SamplerSpec{SamplerKind::AdaptiveMetropolis, /*learn_sd=*/false, ...};
  } else {
    Rcpp::stop("Unknown sampler_type: '%s'", sampler_type.c_str());
  }
}

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);
  }
}

The SamplerSpec also records which diagnostics the sampler produces and whether proposal-SD learning applies, which the warmup schedule consumes.

All three model types (GGM, OMRF, mixed MRF) default to NUTS. The user can select Metropolis via update_method = "adaptive-metropolis".

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).

Constrained integration dispatch

When model.has_constraints() returns true, NUTSSampler switches from the unconstrained leapfrog to RATTLE integration:

Unconstrained path — Uses get_vectorized_parameters() and logp_and_gradient() on the active parameter subset. Used by the OMRF, which handles edge selection by adjusting its active parameter dimension instead.
Constrained path — Uses get_full_position() and logp_and_gradient_full() on the complete parameter vector. The sampler wraps the model’s project_position and project_momentum methods into lambda callbacks passed to nuts_step().

The RATTLE algorithm itself is documented in Constrained Leapfrog (RATTLE). The model-side projection functions are described in Model Classes — Constraint projection.