hopsy.rhat#

class hopsy.rhat(data, series=0, method='rank', dask_kwargs=None)#

Compute estimate of rank normalized splitR-hat for a set of traces.

The rank normalized R-hat diagnostic tests for lack of convergence by comparing the variance between multiple chains to the variance within each chain. If convergence has been achieved, the between-chain and within-chain variances should be identical. To be most effective in detecting evidence for nonconvergence, each chain should have been initialized to starting values that are dispersed relative to the target distribution.

Parameters:
  • data (numpy.ndarray) – MCMC samples with data.shape == (n_chains, n_draws, dim).

  • series (int) – Compute a series of R-hat statistics every series samples, so R-hat will be computed for data[:,:n] for n in range(series, n_draws+1, series). For the default value series==0, R-hat will be computed only once for the whole data.

  • method (str) – Select R-hat method. Valid methods are: - “rank” # recommended by Vehtari et al. (2019) - “split” - “folded” - “z_scale” - “identity”

  • n_procs (int = 1) – In combination with “series”: compute series of R-hat in parallel using n_procs subprocesses.

  • dask_kwargs (dict) – Dask related kwargs passed to wrap_xarray_ufunc().

Returns:

Returns dataset of the potential scale reduction factors, \(\hat{R}\)

Return type:

numpy.ndarray

Notes

The diagnostic is computed by:

\[\hat{R} = \frac{\hat{V}}{W}\]

where \(W\) is the within-chain variance and \(\hat{V}\) is the posterior variance estimate for the pooled rank-traces. This is the potential scale reduction factor, which converges to unity when each of the traces is a sample from the target posterior. Values greater than one indicate that one or more chains have not yet converged.

Rank values are calculated over all the chains with scipy.stats.rankdata. Each chain is split in two and normalized with the z-transform following Vehtari et al. (2019).

Note that if all chains contain the same constant for some parameter (due to equality constraints), rhat will be 1 for this parameter.

References