Diagnosing MCMC: ESS, R-hat, MCSE, and Plots

Adapted from an appendix of my MS thesis.

Runtime Diagnostics

The following numerical and visual diagnostic tools are illustrated with three synthetic posteriors. The first named good_chains was generated from independent and identically distributed (IID) draws from a $\mathrm{Beta}(2,5)$ . Samples that are IID are ideal for approximating the posterior [1].

The second is bad_chains0 generated by sorting good_chains and then adding Gaussian noise. These values are highly autocorrelated (not independent) and not identically distributed since the flattened and sorted array is being reshaped into a 2D array. The third synthetic posterior called bad_chains1 is generated from good_chains by randomly introducing portions where consecutive samples are highly correlated to each other [1].

Effective Sample Size

When using MCMC sampling methods we may wonder if a particular sample is large enough to confidently compute quantities of interest like a mean or variance. This cannot be answered directly because samples from MCMC methods have some degree of autocorrelation. The amount of information contained in autocorrelated samples will be less than what we would get from an IID sample of the same size. The effective sample size (ESS) is an estimator that takes autocorrelation into account and provides the number of draws we would have if our sample was IID [1].

good_chains	bad_chains0	bad_chains1
4.389e+03	2.436	111.1

The ESS when the count of actual samples in our synthetic posteriors is 4000 [1].

Convergence of Markov chains is not uniform across parameter space. It is easier to get a good approximation from the bulk of a distribution than from the tails simply because the tails are dominated by rare events. As a general rule of thumb, an ESS greater than 400 is recommended, otherwise the estimation is unreliable [1].

Potential Scale Reduction Factor (R-hat)

We need ways to estimate convergence for finite samples. One pervasive idea is to run more than one chain starting from very different points and then check the resulting chains to see if they look similar to each other. This intuitive notion can be formalized into a numerical diagnostic known as $\hat{R}$ [1].

The overall idea is that the $\hat{R}$ for the parameter $\theta$ is the standard deviation of all the samples of $\theta$ , that is including all chains together, divided by the root mean square of the separated within chain standard deviations. Ideally we should get a value of 1, as the variance between chains should be the same as the variance within chains. From a practical point of view, values of $\hat{R} \lesssim 1.01$ are considered safe [1].

good_chains	bad_chains0	bad_chains1
1.000	2.408	1.033

The $\hat{R}$ from our synthetic posteriors [1].

Monte Carlo Standard Error

Using MCMC methods introduces an additional layer of uncertainty as we are approximating the posterior with a finite number of samples. We can estimate the amount of error introduced using the Monte Carlo Standard Error (MCSE), which is based on Markov chain central limit theorem. We should check the MCSE only once we are sure ESS is high enough and $\hat{R}$ is low enough, otherwise MCSE is of no use. As with the ESS the MCSE varies across the parameter space and we may want to evaluate it for different regions like specific quantiles [1].

good_chains	bad_chains0	bad_chains1
0.002381	0.1077	0.01781

The MCSE of our example posteriors [1].

Trace Plots

Trace plots are the most common plots in Bayesian literature. They are the first plot after inference for visual inspection. A trace plot is made by drawing the sampled values at each iteration step. From these plots were are able to see if different changes converged to the same distribution, getting a sense of the degree of autocorrelation [1].

Autocorrelation Plots

As see when discussing effective sample size, autocorrelation decreases the amount of information contained in a sample and thus is something we want to keep to a minimum. An autocorrelation plot lets us qualitatively inspect the results from ESS. The figure is the bar plot of the autocorrelation function, where the gray band represents the 95% confidence interval [1].

Rank Plots

Rank plots are histograms of the ranked samples. The ranks are computed by first combining all chains but then plotting the results separately for each chain. If all the chains are targeting the same distribution, we expect the ranks to have a uniform distribution. Additionally, if rank plots of all chains look similar, this indicates good mixing of the chains. Rank plots can be more sensitive than trace plots, are thus are recommended after refinement of sampling and modeling definitions. Trace plots are especially useful during early stages when we need to explore many different alternatives [1].

Divergences

The above diagnostics are summaries of the generated samples. Another way to perform a diagnostic is through monitoring the sampling method. One prominent example are the divergences present in Hamiltonian Monte Carlo (HMC) methods. Divergences can be visualized in trace plots as vertical bars at the bottom of the KDEs indicating where something went wrong during sampling [1].

References

1. Martin, Osvaldo A., Kumar, Ravin, Lao, Junpeng (2021) Bayesian Modeling and Computation in Python. Chapman and Hall/CRC.