In the Bayesian paradigm, it is often the case that a posterior distirbution cannot be found analytically, such as wehn analysis is conjugate. However, this is need not be the end. Markov chain Monte Carlo (MCMC) is a common method used to sample from a posterior distribution when the analysis is not conjugate. The premise of the method is to construct a Markov chain whose stationary distribution is the posterior density. Once the Markov chain has converged to the posterior, any sample generated by the chain will be a realisation from the posterior distribution and posterior beliefs can be inpsected by inspection of the posterior samples.

In this blog common methods to sample from the posterior distribtion are introduced and practical applications of how to asses posterior draws and implement MCMC methods are discussed. It is assumed the reader is familiar with the Bayesian paradigm and somewhat familiar with Markov chains, the two topics are discussed here and here respectively.

Gibbs sampler

The Gibbs sampler, Algorithm 1, has become one of the most common methods to construct an MCMC scheme. It is generally less computationally expensive than the other proposed algorithms. However, it imposes some restrictions on the models which can be used. The algorithm requires analytic expressions for the full-conditional distribution (FCD) of each parameter. An FCD is the distribution of a single parameter, $\theta_i$ (the $i$-th element of $\boldsymbol{\theta}$), conditional on every other entity (parameters and data), $$p(\theta_i |\theta_1, \dots, \theta_{i-1}, \theta_{i+1}, \dots, \theta_m, \boldsymbol{x}).$$

Suppose we wish to infer a set of unknown parameters in a model, whose posterior, $p(\boldsymbol{\theta}|\boldsymbol{x})$, cannot be found analytically nor can it be sampled from. However, the FCD for each element of $\boldsymbol{\theta}$ can be found analytically and sampled from. Then, a Markov chain can be constructed by sequential sampling of the FCDs of each element of $\boldsymbol{\theta}$, and its stationary distribution is the joint-posterior, $p(\boldsymbol{\theta}| \boldsymbol{x})$. The Gibbs sampling algorithm is described in Algorithm 1.

\[ \begin{aligned} &\text{1. Initialise the state of the chain } \boldsymbol{\theta}^{(0)} = (\theta_1^{(0)}, \dots, \theta_m^{(0)})^\text{T} \text{ and set } j=1. \\ &\text{2. Generate } \boldsymbol{\theta}^{(j)} \text{ by sequential realisations from full-conditionals } \\ &\hspace{1cm} \theta_1^{(j)} \sim p(\theta_1 \mid \theta_2^{(j-1)}, \dots, \theta_m^{(j-1)}, \boldsymbol{x}) \\ &\hspace{1cm} \theta_2^{(j)} \sim p(\theta_2 \mid \theta_1^{(j)}, \theta_3^{(j-1)}, \dots, \theta_m^{(j-1)}, \boldsymbol{x}) \\ &\hspace{2cm} \vdots \\ &\hspace{1cm} \theta_m^{(j)} \sim p(\theta_m \mid \theta_1^{(j)}, \dots, \theta_{m-1}^{(j)}, \boldsymbol{x}) \\ &\text{3. Increment } j \text{ by 1 and return to step 2.} \end{aligned} \]

The Gibbs sampler produces an $m$-dimensional Markov chain which will eventually, after a burn-in period, samples from the joint-posterior density. However, the time taken to reach the stationary distribution is not known, and the quality of the posterior draw produced is not guaranteed. How to assess both aspects of inference is discussed later.

Metropolis-Hastings

The Metropolis-Hastings algorithm can be used when sampling from full-conditionals is not available, i.e. the FCD can not be sampled from directly. Unlike the Gibbs sampler, the method requires new $\boldsymbol{\theta}$ values to be actively proposed before being accepted or rejected. The acceptance probability combines information from the proposal distribution and posterior density. The proposal distribution, commonly denoted $q(\cdot \mid \cdot)$, often depends on the previously accepted value of $\boldsymbol{\theta}$. Following this notation, the algorithm is described in Algorithm 2.

\[ \begin{aligned} &\text{1. Initialise the state of the chain } \boldsymbol{\theta}^{(0)} = (\theta_1^{(0)}, \dots, \theta_m^{(0)})^\text{T} \text{ and set } j=1. \\ &\text{2. Propose a new parameter value} \\ &\hspace{1cm} \boldsymbol{\theta}^{*} \sim q(\cdot \mid \boldsymbol{\theta}^{(j-1)}). \\ &\text{3. Calculate the acceptance probability } \\ &\hspace{1cm} \alpha(\boldsymbol{\theta}^{(j-1)}, \boldsymbol{\theta}^{*}) = \min \left\lbrace 1, \frac{p(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{p(\boldsymbol{\theta}^{(j-1)}\mid\boldsymbol{x})} \frac{q(\boldsymbol{\theta}^{(j-1)}\mid\boldsymbol{\theta}^{*})}{q(\boldsymbol{\theta}^{*} \mid \boldsymbol{\theta}^{(j-1)})} \right\rbrace \\ &\hspace{2.8cm} = \min \left\lbrace 1, \frac{p(\boldsymbol{\theta}^{*})L(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{p(\boldsymbol{\theta}^{(j-1)})L(\boldsymbol{\theta}^{(j-1)}\mid\boldsymbol{x})} \frac{q(\boldsymbol{\theta}^{(j-1)}\mid\boldsymbol{\theta}^{*})}{q(\boldsymbol{\theta}^{*} \mid \boldsymbol{\theta}^{(j-1)})} \right\rbrace. \\ &\text{4. Set } \boldsymbol{\theta}^{(j)} = \boldsymbol{\theta}^{*} \text{ with probability } \alpha(\boldsymbol{\theta}^{(j-1)}, \boldsymbol{\theta}^{*}) \text{ otherwise set } \boldsymbol{\theta}^{(j)} = \boldsymbol{\theta}^{(j-1)} \\ &\text{5. Increment } j \text{ by 1 and return to step 2.} \end{aligned} \]

Note that the acceptance probability contains a ratio of posterior densities, the ratio cancels out the normalising constant, and thus is equal to a ratio of prior multiplied by likelihood. The form of the acceptance probability in Algorithm 2 ensures that the Markov chain’s stationary distribution is the posterior density of interest, with proposed values in areas of higher posterior density having a higher acceptance probability.

The proposal distribution, $q(\cdot|\boldsymbol{\theta})$, is of particular importance, as it determines the algorithm’s efficiency. Ideally, the proposal distribution closely resembles the stationary distribution (posterior distribution). If this is not the case, the chain may have inefficient exploration of the parameter space, which can lead to long burn-in periods and increased computational cost. Common proposal distributions are discussed in next. Similarly to the Gibbs sampler, no guarantees are made about the quality of posterior draws or burn-in length.

Independent proposal

An independent proposal is one that does not rely on the current state of the Markov chain. The proposal density can then be written $q(\boldsymbol{\theta}^{*})$, without the $\boldsymbol{\theta}^{(j-1)}$ dependence, and the acceptance probability becomes

\[ \begin{aligned} \alpha(\boldsymbol{\theta}, \boldsymbol{\theta}^{*}) &= \min \left\lbrace 1, \frac{p(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{p(\boldsymbol{\theta}\mid\boldsymbol{x})} \frac{q(\boldsymbol{\theta})}{q(\boldsymbol{\theta}^{*})} \right\rbrace, \\ &= \min \left\lbrace 1, \frac{p(\boldsymbol{\theta}^{*})L(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{p(\boldsymbol{\theta})L(\boldsymbol{\theta}\mid\boldsymbol{x})} \frac{q(\boldsymbol{\theta})}{q(\boldsymbol{\theta}^{*})} \right\rbrace. \end{aligned} \]

Note that the superscript of the previously accepted value of $\boldsymbol{\theta}^{(j-1)}$ has been dropped for brevity. A special case of the independent proposal is when the prior distribution is used to propose values. In this case, the acceptance ratio simplifies to become the ratio of the likelihoods,

\[ \begin{aligned} \alpha(\boldsymbol{\theta}, \boldsymbol{\theta}^{*}) &= \min \left\lbrace 1, \frac{p(\boldsymbol{\theta}^{*})L(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{p(\boldsymbol{\theta})L(\boldsymbol{\theta}\mid\boldsymbol{x})} \frac{p(\boldsymbol{\theta})}{p(\boldsymbol{\theta}^{*})} \right\rbrace, \\ &= \min \left\lbrace 1, \frac{L(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{L(\boldsymbol{\theta}\mid\boldsymbol{x}) } \right\rbrace. \end{aligned} \]

A prior proposal distribution is efficient if the prior closely resembles the posterior. Often, this is not the case, and the algorithm leads to a low acceptance probability and slow exploration of the parameter space.

Symmetric proposal

A proposal is considered symmetric if $q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{}) = q(\boldsymbol{\theta}^{}\mid \boldsymbol{\theta})$ for all $\boldsymbol{\theta}, \boldsymbol{\theta}^{*} \in \Theta$. When this is true, the acceptance probability simplifies to a ratio of the stationary distributions, as such

\[ \begin{aligned} \alpha(\boldsymbol{\theta}, \boldsymbol{\theta}^{*}) &= \min \left\lbrace 1, \frac{p(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{p(\boldsymbol{\theta}\mid\boldsymbol{x})} \right\rbrace \\ &= \min \left\lbrace 1, \frac{p(\boldsymbol{\theta}^{*})L(\boldsymbol{\theta}^{*}\mid\boldsymbol{x})}{p(\boldsymbol{\theta})L(\boldsymbol{\theta}\mid\boldsymbol{x})} \right\rbrace. \end{aligned} \]

Symmetric proposals are common when constructing MCMC schemes, the most common of which is the random walk proposal.

Random walk proposal

A special case of the symmetric distribution is a random walk, which adds independent and identically distributed random noise to the previously accepted parameter values. That is,

\[ \boldsymbol{\theta}^{*} = \boldsymbol{\theta}^{(j-1)} + \boldsymbol{\varepsilon}_j, \]

where \(\boldsymbol{\varepsilon}_j\) are independent and identically distributed random variables. Usually, the noise is normally distributed and centred around \(\boldsymbol{0}\) (\(m\)-dimensional zero vector). Let

\[ \boldsymbol{\varepsilon}_j \sim \text{N}\left(\boldsymbol{0}, \Sigma\right), \]

for all \(j\), then the a proposed value of \(\boldsymbol{\theta}\) is drawn from

\[ \boldsymbol{\theta}^{*}\mid\boldsymbol{\theta}^{(j-1)} \sim \text{N}\left(\boldsymbol{\theta}^{(j-1)}, \Sigma\right). \]

It remains to choose a covariance matrix, $\Sigma$, the choice of which is important for the efficiency of the inference scheme. Similarly to before, a quick exploration of parameter space is preferred. To achieve this, a covariance matrix is needed that shows similar correlations to the posterior distribution and marginal variances which are not too big nor too small. Small marginal variances lead to slow exploration of the parameter space with many proposed values being accepted. Large marginal variances will lead to too few proposed values being accepted and slow exploration.

A common rule of thumb is that the covariance matrix should depend on the posterior covariance and the number of parameters being inferred, such that $$ \text{Var} (\boldsymbol{\varepsilon}_j) = \frac{2.38^2}{m}\widehat{\text{Var}} (\boldsymbol{\theta}|\boldsymbol{x}). $$ Here, $\widehat{\text{Var}}(\boldsymbol{\theta}|\boldsymbol{x})$ denotes a sample variance taken from an initial run of the inference scheme. Initial runs of inference would add additional computational cost to the scheme. However, for complex models with complex posterior distributions, the benefits of a more efficient sampling algorithm would generally outweigh the additional cost.

Metropolis-within-Gibbs

Hybrid MCMC schemes are capable of implementing Metropolis-Hastings updates to a subset of parameters whose FCDs are not tractable (solvable) while allowing component-wise Gibbs updates for those that are.

The Metropolis-within-Gibbs algorithm is a combination of Algorithms 1 and 2. Prior to inference, a proposal distribution for the Metropolis-Hastings update(s) must be chosen. The chain is initialised before new values are drawn for the parameters with tractable FCDs. The Metropolis-Hastings update follows Algorithm 2; new values are proposed, an acceptance probability is calculated, and the chain is updated accordingly. A proposal distribution is only required for parameters whose FCDs can not be sampled from. Note that the updates can occur in any order, and the Metropolis-Hastings update can be separated into several distinct updates, each with its own proposal distribution and acceptance probability. Such schemes are beneficial as high dimensional parameter spaces are often more difficult to explore, with increased difficulty in finding appropriate proposal distributions.

Convergence and autocorrelation

The period before the chain reaches the stationary distribution is referred to as the ‘burn-in’ period, and these samples are removed before analysing the posteriors. Often, multiple chains are executed to ensure they converge to the same distribution, and checks are implemented to ensure this.

Convergence

Once a suitable inference scheme is constructed, a key question remains to be answered; has the chain reached the target distribution? Despite guarantees that the chain will reach a stationary distribution, as the number of iterations approaches infinity, it is not practical to run an MCMC scheme for an infinite number of iterations. Therefore, practical applications of how to answer the question must be considered. Not providing a suitable answer may have dramatic consequences on the results of analysis, skewing posterior distributions, predictions and any conclusions drawn from them.

An initial convergence check is usually done by eye. A converged chain should sample from a steady distribution. Conversely, a chain that does not maintain some equilibrium shows signs of non-convergence and may still be in its burn-in period. Whether or not the chain is moving centred upon some steady state should be visible in a trace plot, a time series line plot of the chain’s value at each iteration. Also informal, although more convincing, is the inspection of multiple chains, initialised at a range of values. All chains should reach the same stationary distribution, and overlayed trace plots are an easy method to highlight when this is or is not the case. When the chains do not overlap, this could indicate a lack of convergence. Multiple chains have the additional benefit that the posterior draws can be combined, provided checks have been made that they come from the same distribution, essentially parallelising inference.

Several formal convergence diagnostics have been proposed. The most widely used method is a likely to be a statistic which compares the within and between chain variation of multiple chains. The statistic, $\hat{\text{R}}$, tends to 1.0 as the number of iterations increases. Many other convergence diagnostics have been suggested, however, they will not be discussed in detail here.

Autocorrelation

Autocorrelation refers to the correlation between realisations of the Markov chain at varying lags, the number of iterations between two realisations, and consequently draws from the posterior of interest. Ideally, all posterior draws would be independent but this is generally not possible. The autocorrelation can be reduced by only saving every $k$-th realisation from the chain, and deleting the rest, referred to as thinning. However, this ignores the information that could be gained from the deleted realisation. A better method is to increase the number of iterations to gain a better understanding of the posterior distribution and reduce the uncertainty in posterior statistics, e.g. mean and variance. Thinning the output can be reserved for when an extremely large number of iterations are needed and the associated storage costs become too high.

The amount of autocorrelation can observed informally by plotting the autocorrelation at a range of lags. More formally, the effective sample size (ESS) can be calculated. This is the number of independent realisations which have the same estimation power as the correlated realisations from the Markov chain. When autocorrelation is low, ESS and the number of realisations will be very close in value. When autocorrelation is high the ESS can be drastically smaller. The ESS gives a good indication of the accuracy of moments calculated from the posterior draws. Single and multi-variate ESS calculations exist and both can used together.

Statistical software

Constructing an appropriate inference scheme can be complicated, and finding appropriate proposal distributions may be difficult for high-dimensional problems. To overcome this hurdle, software exists that can implement a Bayesian inference scheme (semi-)automatically. The two main pieces of software are STAN and JAGS. Both are available within R and Python and they both have advantages and disadvantages. The largest disadvantage for both is the initial time spent learning how to write a model in their respective probabilistic programming languages. Although how to use them is not within the scope of this blog post, a number of high quality resources exist and a quick google search should yield a number of how-to’s in your chosen programming language.

STAN implements another Bayesian inference algorithm called Hamiltonian Monte Carlo (HMC), which generally offers better exploration of the parameter space. This allows STAN to quickly find the main areas of the target support, and once there, the samples generated are almost uncorrelated. However, due to the specific algorithm STAN uses, categorical variables cannot directly be a part of the model. For many models this is not an issue and a ‘workaround’ can be implemented. Again, this will not be discussed within this blog but a number of free online resources cover this. JAGS is not as efficient as STAN in its exploration of the parameter space and its posteriors often produce more correlated samples. However, the initial learning curve is less steep and it allows categorical variables within the model.

Both STAN and JAGS are likely faster than hard-coded inference schemes, both using optimised C code. However, it is not always possible to use them as models becoome more complex.

Final remarks

This introduction into MCMC is fairly brief and omits a lot of the detailed arguments and background which would be a part of any university’s ‘Introduction to Bayesian Methods’ course. If such mathematical and technical detail is needed many texts books (or other online resources) are available. For example, the book Markov Chain Monte Carlo Stochastic Simulation for Bayesian Inference by Gamerman and Lopes provides a comprehensive introduction to MCMC methods, a PDF of which can be found here.