rhat#

cultivars.bayes.chains.rhat(*chains)[source]#

Rank-normalized split-\(\widehat{R}\) with folding, per quantity.

The potential scale reduction of Gelman and Rubin (1992) compares the spread of chain means with the spread within chains: for \(M\) chains of \(N\) draws with within-chain variance \(W\) and between-chain variance \(B\),

\[\widehat{R} = \sqrt{\frac{\widehat{\operatorname{var}}^{+}}{W}}, \qquad \widehat{\operatorname{var}}^{+} = \frac{N - 1}{N}\, W + \frac{1}{N}\, B .\]

Three modifications from Vehtari et al. (2021) make the statistic usable as a stopping rule. Each chain is split in half before the comparison, so \(M = 2C\) for \(C\) supplied chains and a single run still yields a between-chain term, which is what detects a chain that is still drifting. The draws are replaced by their normal scores, \(z = \Phi^{-1}\!\bigl((r - 3/8) / (S + 1/4)\bigr)\) for the pooled rank \(r\) among all \(S\) draws, which makes the statistic invariant to monotone transformations and defined for quantities with no finite variance. And the statistic is computed a second time on the scores of \(|x - \operatorname{median}(x)|\); this folded version is sensitive to chains that agree in location but not in scale, which the plain version misses. The larger of the two is reported.

Parameters:

*chains (ArrayLike) – One or more chains of the same quantity from independent runs, each (S,) for a scalar or (S, d) for d quantities drawn jointly. All chains must have the same shape, with S at least eight and finite throughout.

Return type:

float | NDArray[double]

Returns:

A float for (S,) chains, else (d,). nan for a quantity whose draws never moved; the report built by convergence() lists such quantities as degenerate rather than counting them for or against the verdict. Values above about 1.01 indicate the chains, or the halves of a chain, disagree.

Raises:

Note

With a single chain the between-chain term compares its two halves, which detects drift but not a chain confined to one of several posterior modes. The two half-chains are also short: an independent sequence of 200 draws exceeds 1.01 about a fifth of the time, one of 500 about a fiftieth. A value near the threshold from one short chain is a reason to run longer or add a second seed, not a verdict.

See also

ess_bulk(): The companion size statistic on the same

rank-normalized split chains.

convergence(): The full report, which applies the threshold.

References

Gelman, A., & Rubin, D. B. (1992). Inference from iterative

simulation using multiple sequences. Statistical Science, 7(4), 457-472.

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 that target the same distribution agree:

>>> import numpy as np
>>> rng = np.random.default_rng(1)
>>> float(rhat(rng.standard_normal(1000), rng.standard_normal(1000))) < 1.02
True

Two runs stuck in different places do not:

>>> float(rhat(rng.standard_normal(1000), 3.0 + rng.standard_normal(1000))) > 1.5
True

Rank normalization makes the statistic indifferent to heavy tails, so a Cauchy quantity is assessed like any other:

>>> float(rhat(rng.standard_cauchy(2000), rng.standard_cauchy(2000))) < 1.02
True
On this page