MCMC Diagnostics

PINTS provides a number of functions to diagnose MCMC progress and convergence.

class pints.MCMCSummary(chains, time=None, parameter_names=None)

Calculates and prints key summaries of posterior samples and diagnostic quantities from MCMC chains.

These include the posterior mean, standard deviation, quantiles, rhat, effective sample size and (if running time is supplied) effective samples per second.

Parameters:

• chains – An array or list of chains returned by an MCMC sampler.

• time (float) – The time taken for the run, in seconds (optional).

• parameter_names (sequence) – A list of parameter names (optional).

References

[1]

“Inference from iterative simulation using multiple sequences”, A Gelman and D Rubin, 1992, Statistical Science.

[2] (1,2)

“Bayesian data analysis”, 3rd edition, CRC Press., A Gelman et al., 2014.

chains()

Returns posterior samples from all chains separately.

ess()

Return the effective sample size for each parameter as defined in [2].

ess_per_second()

Return the effective sample size (as defined in [2]) per second of run time for each parameter.

This is only defined if a run time was passed in at construction time, if no run time is known None is returned.

mean()

Return the posterior means of all parameters.

quantiles()

Return the 2.5%, 25%, 50%, 75% and 97.5% posterior quantiles.

rhat()

Return Gelman and Rubin’s rhat value as defined in [1]. If a single chain is used, the chain is split into two halves and rhat is calculated using these two parts.

std()

Return the posterior standard deviation of all parameters.

summary()

Return a list of the parameter name, posterior mean, posterior std deviation, the 2.5%, 25%, 50%, 75% and 97.5% posterior quantiles, rhat, effective sample size (ess) and ess per second of run time.

time()

Return the run time taken for sampling.

pints.rhat(chains, warm_up=0.0)

Returns the convergence measure \(\hat{R}\) for the approximate posterior according to [3].

\(\hat{R}\) diagnoses convergence by checking mixing and stationarity of \(m\) chains (at least two, \(m\geq 2\)). To diminish the influence of starting values, the first portion of each chain can be excluded from the computation. Subsequently, the truncated chains are split in half, resulting in a total number of \(m'=2m\) chains of length \(n'=(1-\text{warm_up})n/2\). The mean of the variances within and between the resulting chains are computed, \(W\) and \(B\) respectively. Based on those variances an estimator of the marginal posterior variance is constructed

\[\widehat{\text{var}}^+ = \frac{n'-1}{n'}W + \frac{1}{n'}B,\]

The estimator overestimates the variance of the marginal posterior if the chains are not well mixed and stationary, but is unbiased if the original chains equal the target distribution. At the same time, the mean within variance \(W\) underestimates the marginal posterior variance for finite \(n\), but converges to the true variance for \(n\rightarrow \infty\). By comparing \(\widehat{\text{var}}^+\) and \(W\) the mixing and stationarity of the chains can be quantified

\[\hat{R} = \sqrt{\frac{\widehat{\text{var}}^+}{W}}.\]

For well mixed and stationary chains \(\hat{R}\) will be close to one.

The mean within \(W\) and mean between \(B\) variance of the \(m'=2m\) chains of length \(n'=(1-\text{warm_up})n/2\) are defined as

\[W = \frac{1}{m'}\sum _{j=1}^{m'}s_j^2\quad \text{where}\quad s_j^2=\frac{1}{n'-1}\sum _{i=1}^{n'}(\psi _{ij} - \bar{\psi} _j)^2,\]

\[B = \frac{n'}{m'-1}\sum _{j=1}^{m'}(\bar{\psi} _j - \bar{\psi})^2.\]

Here, \(\psi _{ij}\) is the jth sample of the ith chain, \(\bar{\psi _j}=\sum _{i=1}^{n'}\psi _{ij}/n'\) is the within chain mean of the parameter \(\psi\) and \(\bar{\psi } = \sum _{j=1}^{m'}\bar{\psi} _{j}/m'\) is the between chain mean of the within chain means.

References

[3]

“Bayesian data analysis”, ch. 11.4 ‘Inference and assessing convergence’, 3rd edition, Gelman et al., 2014.

Parameters:

• chains (np.ndarray of shape (m, n) or (m, n, p)) – A numpy array with \(n\) samples for each of \(m\) chains. Optionally the \(\hat{R}\) for \(p\) parameters can be computed by passing a numpy array with \(m\) chains of length \(n\) for \(p\) parameters.

• warm_up (float) – First portion of each chain that will not be used for the computation of \(\hat{R}\).

Returns:

rhat – \(\hat{R}\) of the posteriors for each parameter.

Return type:

float or np.ndarray of shape (p,)

pints.rhat_all_params(chains)

Deprecated alias of rhat().

pints.effective_sample_size(samples)

Calculates effective sample size (ESS) for a list of n-dimensional samples.

Parameters:

samples – A 2d array of shape (n_samples, n_parameters).