Gibbs Sampling by Example

In order to grasp a new technique it’s often useful to see how it is used in real-world examples. Here we will examine a particular implementation of a Gibbs sampler which pertains to a model we discussed in a previous series of articles.


As was discussed in the previous article, the Gibbs sampling algorithm is a Markov Chain Monte Carlo (MCMC) technique that can be used to simulate outcomes from a multivariate distribution.

As with all MCMC techniques, simulations are produced by generating a realisation of a random walk and recording the position of the walk at the end of consecutive time intervals. MCMC techniques are designed such that, no matter where it starts, the positions produced by a long enough walk will resemble a sample from the target distribution.

The Gibbs sampler produces a random walk by updating each coordinate associated with the multivariate distribution sequentially. Once each coordinate has been updated the new values are recorded and another iteration occurs.

The key to the Gibbs sampler is in how each of the coordinates are updated. In particular, updating a particular variable involves simulating an outcome from its conditional distribution, given the latest values of all the other variables.

The animation below illustrates how this process unfolds in the case of a bivariate Normal distribution.

Here you can see that the random walk moves around the space by first taking a horizontal step and then a vertical step.

In the animation I’ve also plotted the positions of the random walk on top of some bivariate Normal contours. One can see that after a large number of steps the distribution of the dots resembles the target distribution.

In this article we will show how Gibbs sampling can be applied to another more specialised problem.

The Model

In the Incomplete Data series I introduced a model that could be used to describe the timing of individual claim payments made by an insurer.

From the data we made two observations regarding the timing of payments; payments occurred sporadically and eventually they ceased.

These phenomena were assumed to be caused by the state of each claim. In particular, we assumed that the claim could be in one of three states: an open state in which no payment occurs; a closed state, also with no payment; and an open state in which a payment does occur. In addition, it was assumed that if a claim closed it remained closed.

The probability of a claim moving from an open state to the closed state (  p ) and, given it was open, giving rise to a payment (  q ) were both assumed to be constant and applicable to all claims.

To understand more about the claim payment process, therefore, it’s useful to investigate what can be learned about the unknown factors in our model based on the observed data. In this case, the unknowns are the transition probability parameters (  p and   q ) and if and when any of the claims closed.

As was discussed in the Incomplete Data series, this problem can be tackled by looking at the posterior distribution of the unknowns. Specifically, assuming a uniform prior on   p and   q , the posterior density takes the following form:

    \begin{equation*} \begin{split} P(p,q,\mathbf{c},& \mathbf{z^{open}}) \propto \\ & \prod\limits_{i=1}^{m}{(1-p)^{z^{open}_{i} + n_{i} + d_{i}} \cdot p^{c_{i}} \cdot q^{d_{i}} \cdot (1-q)^{z^{open}_{i}+n_{i}} } \end{split} \end{equation*}

Here, the event that a given claim is closed is represented by   c_{i} ,   z^{open}_{i} represents the number of trailing periods with no payments that each claim remained open for, and   p and   q are the parameters we discussed previously.

The observed data is represented by the number of claims (  m), the number of periods in which payments occurred (  d_{i} ), the number of periods prior to the last payment in which no payments occurred (  n_{i} ), and the number of observed periods after the last payment (  z_{i} ).

In this article we will use a Gibbs sampler to simulate from this distribution. However, as discussed in the previous article, to do this we need to know what each of the unknowns’ conditional distributions are.

The Conditional Distributions

The purpose of this article isn’t to dwell on how we arrive at the conditional distribution of each unknown. After all, in your problem the posterior and conditional distributions will undoubtedly be different. However, as they are integral to performing Gibbs sampling it’s worth getting a feel for how you might derive them.

In this case you can get most of the way to deriving the conditional distributions by “inspection”. This means considering the posterior density as a function of each unknown and seeing whether it looks like any of the standard distributions you already know (see Wikipedia for a list of standard distributions).

For example, as a function of (  p) the posterior looks like:

    \begin{equation*} P(p|q,\mathbf{c},\mathbf{z^{open}}) \propto (1-p)^{\sum\limits_{i=1}^{m}{z^{open}_{i}+n_{i}+d_{i}}} \cdot p^{\sum\limits_{i=1}^{m}{c_{i}}} \end{equation*}

This density is of the same form as a beta distribution.

By a similar albeit more long-winded approach it’s possible to see that the conditional distributions of each of the unknowns take the form of the following standard distributions:


As you can see the conditional distributions each take a standard form. The exception to this being   z^{open}_{i} –which is either a constant or follows a truncated Geomtric distribution–nevertheless it is still easy to simulate from.

With these conditional distributions, we have all that we need to implement a Gibbs sampler.

The Gibbs Sampler

At this point it may be obvious to you how we would set up a Gibbs sampler using the conditional distributions above. However, I expect others might feel uncomfortable translating what we’ve discussed into an actual implementation.

With this in mind, below is the Gibbs sampler relevant to this article written in “pseudo-code”:


Hence, in the programming language of your choice, one iteration of the Gibbs sampler involves: for each claim simulating whether it is closed or not (step 7); simulating how many periods each claim has been open for (step 9); simulating new parameter values for   p and   q (steps 12 and 13).

Final Word

In this article we discussed a specialised example of a Gibbs sampling algorithm. In particular, we examined how to implement a Gibbs sampler for a model that could be used to explain the timing of payments on insurance claims.

I hope that by working through a real-world example applying Gibbs sampling to your problems will seem less daunting.

If you have any questions on this article or how it pertains to your problem feel free to leave a comment below.

Leave a Reply

Your email address will not be published. Required fields are marked *