PyBNF convergence diagnostics (pybnf.diagnostics)

Convergence diagnostics for MCMC samplers — pure R-hat / ESS math.

Extracted from algorithms/samplers/base.py (M2.2 move 5, ADR-0009): the convergence diagnostics are substantial, already shared by all four samplers, project-specific, and oracle-testable, so their pure numerical core lives here as free functions — a peer of objective.py, navigable, and importable by the benchmark harness (M2.5) without reaching into samplers/. The instance-coupled glue (report_convergence_diagnostics / check_convergence / _write_diagnostics) and the PSet→array bridge (_param_vec) stay on BayesianAlgorithm, which delegates the math here.

Implements the rank-normalized split-R-hat and bulk/tail effective sample size of Vehtari, Gelman, Simpson, Carpenter & Bürkner (2021), Bayesian Analysis (the same conventions as Stan / ArviZ).

The data layout throughout: chain_history is a list of num_parallel per-chain histories, each a list of (n_dim,) parameter vectors in the sampling space (see BayesianAlgorithm._param_vec). The internal chains arrays are (n_chains, n_draws, n_dim) (or (n_chains, n_draws) for the single-parameter ESS kernel).

pybnf.diagnostics.ess(chain_history, num_parallel)

Compute bulk and tail effective sample size per Vehtari et al. (2021).

Bulk ESS: computed on rank-normalized values (same transform as R-hat). Tail ESS: minimum ESS of the 5% and 95% quantile indicators.

Returns (bulk_ess, tail_ess) arrays of shape (n_dim,) or (None, None).

pybnf.diagnostics.ess_from_chains(chains)

Compute effective sample size from an (M, n) array of chains using FFT-based autocovariance and Geyer’s initial positive sequence estimator.

pybnf.diagnostics.rhat(chain_history, num_parallel)

Compute rank-normalized split-R-hat for each parameter (Vehtari, Gelman, Simpson, Carpenter & Burkner, 2021, Bayesian Analysis).

Steps: 1. Split each chain in half (doubles the number of chains, catches within-chain non-stationarity) 2. Rank-normalize across all split chains (replaces values with normal scores of their ranks) 3. Compute R-hat on both the ranked values and folded ranked values (detects scale differences) 4. Return the element-wise maximum

Returns a numpy array of shape (n_dim,) or None if insufficient data.

pybnf.diagnostics.split_chain_rhat(chains)

Compute the potential scale reduction factor R-hat from an array of chains, using the Vehtari et al. (2021) convention R = sqrt(var_plus / W), where var_plus = ((n-1)/n) W + B/n. This is the form paired with the rank normalization and folding done by rhat() (and matches Stan / ArviZ).

It deliberately omits the older Gelman & Rubin (1992) df-style correction sqrt((N+1)/N * (var_plus/W) - (n-1)/(N n)): that factor only inflates R-hat when the chains have not converged (it cancels to ~1 at convergence), making PyBNF report systematically higher R-hat than reference tools on the same chains – harmful for cross-tool comparison – without changing any convergence decision.

chains: (N, n, d) array Returns: (d,) array of R-hat values

pybnf.diagnostics.split_chains(chain_history, num_parallel)

Build the split-chains array used for R-hat and ESS.

Uses the last 50% of each chain, then splits that window in two (doubling the chain count, catching within-chain non-stationarity). Returns an (2 * num_parallel, half, n_dim) array, or None if there is too little history (fewer than 20 recorded steps, or fewer than 5 per split half).