cultivars.bayes.chains#
Convergence diagnostics for Markov chain Monte Carlo output.
A sampler’s draws are only a posterior once the chain has forgotten where it started and has visited the distribution in proportion, and neither is observable directly. The diagnostics here are the standard indirect evidence, following Vehtari, Gelman, Simpson, Carpenter, and Bürkner (2021): rank-normalized split-\(\widehat{R}\) with folding, which compares chains, or the two halves of one chain, on location and on scale; bulk and tail effective sample sizes by Geyer’s initial monotone sequence, which say how many independent draws the correlated ones are worth for a mean and for an interval endpoint; the Monte Carlo standard error of the mean; and Geweke’s early-versus-late score, which asks whether the burn-in was long enough. For \(S\) pooled draws with integrated autocorrelation time \(\tau\), the effective size is \(S / \tau\), and the working thresholds are \(\widehat{R} \le 1.01\) and at least 100 effective draws per chain.
Every sampled result in the package – the stochastic-volatility
posteriors, the Gibbs, Student-t, stochastic-volatility, and hierarchical
BVARs, the time-varying-parameter models, the factor and structural
volatility models, the particle-chain DSGE posterior – carries a
convergence() method that assesses its retained draws and returns a
ConvergenceTest. This module is the same machinery for draws that
did not come from the package: a chain from another sampler, or a
function of the package’s draws (an impulse response at a horizon, a
variance ratio, a ratio of two parameters) whose own convergence is the
question. Each function takes one or more chains of the same quantity: a
chain is (S,) for a scalar or (S, d) for d quantities drawn
jointly, and several chains are several positional arguments from
independent runs.
Two commitments shape the surface. First, a single chain is assessed by
splitting it in half, which detects drift but cannot detect a chain
confined to one of several posterior modes; a second run from a different
seed is what detects that, and passing one is the recommended practice
rather than an option, so every function takes *chains. Second, a
quantity whose draws never moved is reported as degenerate and is neither
passed nor failed: the diagnostics cannot tell a parameter fixed by
construction from a sampler that is stuck, and the report says which
quantities it declined to judge rather than guessing.
Layout. The six public functions are thin: each stacks its chains through
_core’s _stack and applies a per-quantity statistic from
_core._estimators – _rhat, _ess_bulk, _ess_tail,
_mcse_mean, _geweke – which in turn share
_potential_scale_reduction, _effective_sample_size, and the
FFT _autocovariance; the rank normalization and the split are
_core._transforms. convergence() hands the stack to
ConvergenceTest._assess(), the classmethod on the frozen record in
_internals, which computes every column and applies the thresholds;
the convergence() method of every sampled result reaches the same
classmethod through _ConvergenceMixin, so
a report built here and one built by a result are the same object and
render the same table. The thresholds _RHAT_TOL and
_MIN_ESS_PER_CHAIN, and the eight-draw minimum chain length
_MIN_CHAIN_DRAWS, are _core defaults.
References
Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457-472.
Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In J. M. Bernardo, J. O. Berger, A. P. Dawid, & A. F. M. Smith (Eds.), Bayesian Statistics 4 (pp. 169-193). Oxford University Press.
Geyer, C. J. (1992). Practical Markov chain Monte Carlo. Statistical Science, 7(4), 473-483.
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. Bayesian Analysis, 16(2), 667-718.
Example
Two runs of two quantities, assessed together:
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> chain_a = rng.standard_normal((2000, 2))
>>> chain_b = rng.standard_normal((2000, 2))
>>> report = convergence(chain_a, chain_b, names=("alpha", "beta"))
>>> report.converged
True
>>> report.n_chains, report.n_draws
(2, 2000)
Two runs that settled in different places are caught by \(\widehat{R}\), which is what the second seed is for:
>>> stuck = np.stack([chain_a[:, 0], chain_b[:, 0] + 4.0])
>>> rhat(stuck[0], stuck[1]) > 1.5
True