Convergence Diagnostics

This page documents the C++ implementation of convergence diagnostics in bgms. The ESS and R-hat computations are implemented in C++ for speed, enabling fast diagnostics even with long chains and many parameters. For how to interpret these diagnostics as a user, see MCMC Diagnostics in the Guide; the user-facing accessors are extract_ess() and extract_rhat().

Overview

After sampling, bgms computes convergence diagnostics for all parameters:

  • Effective sample size (ESS) — equivalent number of independent draws
  • Split-R-hat — between-chain vs within-chain variance ratio, computed on half-chains (Vehtari et al., 2021)
  • Monte Carlo standard error (MCSE) — precision of the posterior mean estimate

The C++ implementation matches coda::effectiveSize for ESS and computes the classic split-R-hat for convergence, running in parallel across parameters for speed; the R layer splits each chain in half before calling the R-hat kernel.

C++ implementation

The diagnostic computations are in src/mcmc_diagnostics.cpp. The main entry points are:

ESS via AR spectral density

// [[Rcpp::export(.compute_ess_cpp)]]
Rcpp::NumericVector compute_ess_cpp(Rcpp::NumericVector array3d);

Computes ESS for a 3D array [niter x nchains x nparam]. Multi-chain ESS is the sum of per-chain ESS values.

The per-chain algorithm follows Plummer et al. (2006):

  1. Compute biased autocovariance \(c[0..L]\) where \(L = \lfloor 10 \log_{10}(n) \rfloor\)
  2. Fit autoregressive (AR) models of orders \(1..L\) via the Levinson-Durbin recursion; \(v_k\) denotes the innovation (residual) variance of the order-\(k\) fit
  3. AIC order selection: \(\arg\min_k \{ n \log(v_k) + 2k \}\)
  4. Spectral density at 0: \(\text{spec}_0 = v_{\text{best}} / (1 - \sum_j a_{\text{best},j})^2\), where \(a_{\text{best},j}\) are the AR coefficients of the selected order
  5. ESS = \(n \cdot \text{var}_{\text{unbiased}} / \text{spec}_0\)

This matches coda::spectrum0.ar + coda::effectiveSize.

Split-R-hat

// [[Rcpp::export(.compute_rhat_cpp)]]
Rcpp::NumericVector compute_rhat_cpp(Rcpp::NumericVector array3d);

Computes R-hat for a 3D array [niter x nchains x nparam]. Returns NA for single-chain input.

The C++ kernel computes the classic statistic (Gelman & Rubin, 1992):

  1. Compute per-chain means \(\bar{x}_c\) and variances \(s^2_c\)
  2. Within-chain variance: \(W = \frac{1}{m} \sum_c s^2_c\)
  3. Between-chain variance: \(B = \frac{n}{m-1} \sum_c (\bar{x}_c - \bar{x})^2\)
  4. Pooled variance estimate: \(\hat{V} = \frac{n-1}{n} W + \frac{1}{n} B\)
  5. R-hat = \(\sqrt{\hat{V} / W}\)

The Brooks–Gelman degrees-of-freedom adjustment that coda::gelman.diag applies is deliberately omitted: on a nearly-saturated binary indicator chain (one half-chain carries a brief excursion, the rest stay constant) the adjustment’s degrees of freedom degenerate and R-hat converges to \(\sqrt{5/3} \approx 1.29\) regardless of the data, making decisive edges look unconverged. Degenerate inputs are handled explicitly instead: all sub-chains constant and equal returns NA, while constant-but-unequal sub-chains (chains stuck in different states) return +Inf, so a real failure is loud rather than silent.

The R layer feeds this kernel split chains: split_chains() in R/mcmc_summary.R halves each chain into two sub-chains (dropping the middle draw of an odd-length chain so the halves match), so the statistic reported everywhere — summary() tables and extract_rhat() — is the split-R-hat of Vehtari et al. (2021). Splitting makes within-chain drift visible: a chain whose first and second half disagree inflates R-hat even when whole-chain means happen to agree across chains.

Mixture ESS for indicators

// [[Rcpp::export(.compute_indicator_ess_cpp)]]
Rcpp::NumericMatrix compute_indicator_ess_cpp(Rcpp::NumericVector array3d);

Computes ESS for binary indicator sequences using the two-state Markov chain approach:

  1. Counts transitions: \(n_{01}\) (off→on), \(n_{10}\) (on→off), \(n_{00}\), \(n_{11}\)
  2. Estimates transition probabilities: \(\hat{a} = n_{01} / (n_{00} + n_{01})\), \(\hat{b} = n_{10} / (n_{10} + n_{11})\)
  3. Returns ESS = \(n \cdot (a + b) / (2 - a - b)\), with \(n\) the chain length as above

This formula derives from the integrated autocorrelation time of a two-state Markov chain (van den Bergh et al., 2026): the lag-\(k\) autocorrelation of the indicator sequence is \(\rho(k) = (1-a-b)^k\), giving integrated autocorrelation time \(\tau = (2 - a - b)/(a + b)\) and hence \(\text{ESS} = n / \tau\).

When the total number of transitions (\(n_{01} + n_{10}\)) is below 5, the estimate is unreliable and flagged in the output.

Indicator summaries and the RB chain

When the fit carries Rao-Blackwellized inclusion draws (see GGM Internals), the printed indicator summary estimates the inclusion probability and its MCSE from that continuous \(J\)-chain (summarize_rb_inclusion()), with ESS and R-hat computed on the same chain; the mixture ESS and the directional transition counts (n0->1, n1->0) continue to come from the binary indicator draws. The two ESS columns are read as a pair: the RB n_eff measures precision conditional on exploration and cannot see a stuck sampler, while n_eff_mixt measures the exploration itself. For edges with no transitions at all — or a machine-constant \(J\)-chain — the RB mcse/n_eff/Rhat are masked to NA, since they would otherwise reassure without evidence. The raw-indicator average remains available through extract_posterior_inclusion_probabilities(fit, estimator = "raw").

Parallelization

Diagnostics are computed in parallel across parameters using RcppParallel (TBB). Each parameter’s ESS and R-hat can be computed independently:

struct ESSWorker : public RcppParallel::Worker {
  void operator()(std::size_t begin, std::size_t end) {
    for(std::size_t j = begin; j < end; j++) {
      double total_ess = 0.0;
      for(int c = 0; c < nchains; c++) {
        const double* col = data + c * niter + j * niter * nchains;
        total_ess += compute_column_ess(col, niter, max_order);
      }
      ess[j] = total_ess;
    }
  }
};

R interface

The diagnostics are exposed to R through: - .compute_ess_cpp() — batch ESS for all parameters - .compute_rhat_cpp() — batch R-hat for all parameters - .compute_indicator_ess_cpp() — batch indicator ESS with transition counts

These are called from R/mcmc_summary.R which assembles the summary tables.

User-facing extractors: - extract_ess(fit) — returns ESS for all parameters - extract_rhat(fit) — returns R-hat for all parameters - summary(fit) — includes ESS, R-hat, and MCSE in tabular output

References

  • Plummer et al. (2006) — AR spectral density ESS (coda package)
  • Gelman & Rubin (1992) — original R-hat formulation
  • van den Bergh et al. (2026) — mixture ESS for spike-and-slab models

References

Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472. https://doi.org/10.1214/ss/1177011136
Plummer, M., Best, N., Cowles, K., & Vines, K. (2006). CODA: Convergence diagnosis and output analysis for MCMC. R News, 6(1), 7–11. https://journal.r-project.org/archive/2006-1/RNews_2006-1.pdf
van den Bergh, D., Clyde, M. A., Raftery, A. E., & Marsman, M. (2026). Reversible jump MCMC with no regrets: Bayesian variable selection using mixtures of mutually singular distributions. Manuscript in Preparation.
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved \(\widehat{R}\) for assessing convergence of MCMC (with discussion). Bayesian Analysis, 16(2), 667–718. https://doi.org/10.1214/20-BA1221