2 Week 7: Introduction to MCMC

Markov Chain Monte Carlo, often referred to simply as MCMC, has probably been the single most important advance in statistics over the last 25 years or so. MCMC has its origins in Metropolis et al. (1953) and Hastings (1970). However, it was only in the late 1980’s with the seminal paper Gelfand and Smith (1990) that the MCMC idea began to invade mainstream statistics and a considerable amount of research has been devoted to it ever since.

MCMC utilises a Monte Carlo method based upon Markov chains (as the name suggests!) MCMC is primarily used in Bayesian statistics to provide samples from the posterior distribution.

We cover three preliminary topics before discussing MCMC. These are:-

1. 1.

   Monte Carlo methods;

2. 2.

   Markov chains

3. 3.

   Bayesian Statistics (very brief overview – details in MATH553).

Monte Carlo methods basically means simulation. These are used primarily to evaluate integrals (for random variables) where analytical solutions are either not possible or extremely laborious.

The simplest example is the case of the expectation

$$𝔼\left[X\right]=\int xf(x)𝑑x$$

of a random variable $X$ which has no analytical form. We can estimate the mean of $X$, $𝔼\left[X\right]$ by taking a sample from the distribution $X$ and using the sample mean as an estimate of the theoretical mean. That is:-

1. 1.

   Take a sample of size $n$, $x_1,x_2,\mathrm{\dots },x_n$ from the distribution $X$, $f(\cdot )$.
   We assume that the observations are independent.

2. 2.

   Then $𝔼\left[X\right]\approx \frac{1}{n}{\sum }_{i=1}^n{x}_i$.

For example, $X$ could represent the height (in cm) of a 5 year old child. An estimate of $𝔼\left[X\right]$ can then be obtained by taking a class of 20, 5 year old children and measuring their heights.

This argument can be generalised. Suppose that we wish to calculate/evaluate, $\theta =𝔼\left[\varphi (X)\right]$, for some function $\varphi (\cdot )$

$$\theta =\int \varphi (x)f(x)𝑑x.$$

Then if we take a sample $x_1,x_2,\mathrm{\dots },x_n$ from the distribution $X$, we have that

$$\widehat{\theta }=\frac{1}{n}\sum _{i=1}^n\varphi (x_i)$$

is an unbiased estimate of $\theta$. This process is exactly the same as for the expectation, in that, we use the sample mean to estimate the theoretical mean. This approach is very easy to use, even for multi-dimensional distributions.

Properties of Monte Carlo integration

1. 1.

   Unbiased.

2. 2.

   If variance of $\varphi (X)$ exists, then

   $$\text{Var}\left[\widehat{\theta }\right]=\frac{1}{n}\text{Var}\left[\varphi (X)\right].$$

3. 3.

   Central limit theorem

   $$\frac{\sqrt{n}(\widehat{\theta }-\theta )}{\sqrt{\text{Var}\left[\varphi (X)\right]}}\stackrel{D}{⟶}N(0,1)$$

   That is,

   $$\widehat{\theta }\approx N(\theta ,\text{Var}\left[\varphi (X)\right]/n).$$

All the above hold provided $\{X_i\}$ are independent. However, these results can be extended to the case where the $\{X_i\}$ are dependent which will be the case when using MCMC.

Drawbacks of Monte Carlo integration

1. 1.

   The variance of $\varphi (X)$ is often large.
   This problem can be circumvented to some extent by variance reduction methods such as importance sampling.

2. 2.

   It can be difficult to simulate from $X$ (in an efficient manner).
   This is a harder problem to deal with and this is where MCMC will be useful. Simulating from univariate distributions (rejection sampling, inversion of the cdf) is usually fairly straightforward but for multivariate distributions the problem becomes a lot harder as methods such as rejection sampling perform much worse as the dimension of the problem grows. This means a lot of work is required to get even a small sample from the distribution $X$.

There is a lot of deep theory related to the use of Markov chains in Monte Carlo methods and we shall barely scratch the surface.

What is a Markov Chain?
A Markov chain $\{X_n,n\ge 0\}$ is a stochastic process which satisfies the ‘memoryless’ property:

$$P(X_n\in A|X_{n-1}=x_{n-1},X_{n-2}=x_{n-2},\mathrm{\dots },X_0=x_0)=P(X_n\in A|X_{n-1}=x_{n-1}),$$

where $X_n$ denotes the state of the process after $n$ steps. Basically, the future is independent of the past, given the present state of the process. A good introduction to the theory of Markov Chains is given in Grimmett and Stirzaker (1992), Section 6. We restrict attention in our discussion to discrete, time-homogeneous Markov chains, that is, for all $A$, $i$ and $n$,

$$P(X_{n+1}\in A|X_n=i)=P(X_1\in A|X_0=i).$$

A very simple example of a Markov chain is the game Snakes and Ladders. Your future progress only depends upon your current position and not how you have arrived at your current position.

We shall illustrate the basic ideas, assuming that the state space for $X$ is countable, although in most examples of MCMC, the state space (set of possible values) of the parameters will be ${ℝ}^d$ which is uncountable. For all $i,j$, let

$$P_{ij}=P(X_1=j|X_0=i).$$

We shall assume that the Markov chain is time homogeneous, that is, for all $n\ge 1$,

$$P(X_n=j|X_{n-1}=i)=P_{ij},$$

the probability of moving from state $i$ to state $j$ in one step. Suppose that the state space is finite. Let $P$ denote the matrix with elements $\{P_{ij}\}$. Then the ${(i,j)}^{th}$ element of $P^n$ gives $P(X_n=j|X_0=i)$.

Snakes and Ladders Example

Suppose that we are currently at square $37$ after $n$ turns and that there is a ladder at square 40 leading to square 62 and a snake at square 42 leading to square 27. Then rolling a standard 6 sided die, we have

	$P(X_{(n+1)}=38|X_n=37)$	$=$	${\displaystyle \frac{1}{6}}$
	$P(X_{(n+1)}=39|X_n=37)$	$=$	${\displaystyle \frac{1}{6}}$
	$P(X_{(n+1)}=62|X_n=37)$	$=$	${\displaystyle \frac{1}{6}}$
	$P(X_{(n+1)}=41|X_n=37)$	$=$	${\displaystyle \frac{1}{6}}$
	$P(X_{(n+1)}=27|X_n=37)$	$=$	${\displaystyle \frac{1}{6}}$
	$P(X_{(n+1)}=43|X_n=37)$	$=$	${\displaystyle \frac{1}{6}}$

depending upon which number we roll from 1 to 6.

A Markov chain is
irreducible; if for all $i$ and $j$, there exists some $n$ such that $P(X_n=j|X_0=i)>0$;
(we can get from any state to any other state)
aperiodic; if the greatest common divisor of $\{n;P(X_n=j|X_0=i)>0\}=1$;
(does not exhibit periodic behaviour)
recurrent; if $P(\text{Markov chain returns to }i|\text{it started at }i)=1$ for all $i$;
(we eventually return to the starting state)
and positive recurrent; if it is irreducible, aperiodic and recurrent, and there exists a unique collection of probabilities $\{{\pi }_j\}$ such that for all $i$,

$\underset{n\to \mathrm{\infty }}{lim}P(X_n=j|X_0=i)={\pi }_j$

(2.1)

and

$$\sum _i{\pi }_i{P}_{ij}={\pi }_j.$$

The distribution $𝝅$ is called the stationary distribution of the Markov chain. Note that (2.1) says, that regardless of the value of $X_0$, for large $n$, $P(X_n=j|X_0=i)\approx {\pi }_j$. That is, for large $n$, $X_n$ is (approximately) distributed to according to the stationary distribution. How large $n$ needs to be for the approximation to be close depends upon the starting value and the rate of convergence. The rate of convergence can be calculated using probability theory, and depends upon how

$$\underset{x\in S,A\in S}{sup}|P(X_n\in A|X_0=x)-\pi (A)|$$

behaves as $n\to \mathrm{\infty }$, where $S$ is the state space of the parameters (often ${ℝ}^d$ or a subset thereof). In practice with MCMC convergence is often ‘detected’ by eye without any formal procedures.

A key point is that if, for all $i$, $P(X_0=i)={\pi }_i$, then for all $n\ge 1$ and $i$, $P(X_n=i)={\pi }_i$. Therefore if $X_0$ is distributed according to $𝝅$, then for all $n\ge 1$, $X_n$ is distributed according to $𝝅$ (the stationary distribution). Also $𝝅=𝝅P$ in matrix notation.

Example. Suppose that we have four states labeled 1, 2, 3, 4.

Transition matrix

Then $p_{ij}$ is the probability of moving from state $i$ to state $j$. That is,

$$P(X_1=j|X_0=i)=p_{ij}.$$

Note that all the rows must sum to 1, i.e. ${\sum }_{j=1}^4{p}_{ij}=1$.

We consider different transition matrices describing the transitions from one time point to the next.

Matrix 1.

This is an example of a reducible Markov chain since we can’t get from state 2 to state 3. Also the Markov chain is not recurrent since if you start in state 3, the probability of being in state 3 at any later time does not exceed 0.5. (If you go to state 1 or 2 at the first step you will never return to state 3.) The probability of ever returning to state 3 is 0.3.

Matrix 2.

This is periodic, it returns to a given state after an even number of steps. (It alternates between odd states, 1 and 3, and even states, 2 and 4.)

Matrix 3.

This Markov chain is positive recurrent with $𝝅=(0.2,0.5,0.1,0.2)$ such that

$$𝝅=𝝅P.$$

For example,

	${\displaystyle \sum _{j=1}^4}{\pi }_j{P}_{j2}$	$=$	$(0.2\times 0.5)+(0.5\times 0.6)+(1\times 0.1)+(0\times 0.2)$
		$=$	$0.5={\pi }_2$

The aim of MCMC is to construct a Markov chain whose stationary distribution is the posterior distribution that we are interested in. Therefore a sample from the Markov chain (provided the Markov chain is started in stationarity) will form a sample from the posterior distribution of interest. It is important to check that the Markov chain we construct is irreducible and aperiodic, in particular the irreducibility of the Markov chain.

Detailed balance

A discrete Markov chain with transition matrix $P$ is said to satisfy detailed balance if there is a distribution $𝝂$ such that for all $i,j$,

$${\nu }_i{P}_{ij}={\nu }_j{P}_{ji}.$$

Lemma If a Markov chain with transition matrix $P$ satisfies detailed balance with distribution $𝝂$ then $𝝂(=𝝅)$ is the stationary distribution of the chain.
Proof: Suppose that $𝝂$ satisfies detailed balance. Then for any $i$,

	${\displaystyle \sum _j}{\nu }_j{P}_{ji}$	$=$	${\displaystyle \sum _j}{\nu }_i{P}_{ij}$
		$=$	${\nu }_i{\displaystyle \sum _j}{P}_{ij}={\nu }_i,$

since for all $i$, ${\sum }_j{P}_{ij}=1$. (All rows sum to 1.) Hence $𝝂$ is the stationary distribution of the Markov chain. $\mathrm{\square }$

Checking detailed balance is often an easy way of confirming a probability distribution $𝝅$ is the stationary distribution of a Markov chain.

Continuous state space Markov chains

Usually, but not always, in Bayesian statistics we want to obtain samples from a continuous distribution (or a distribution which is a mixture of continuous and discrete variables). Most of the above discussion for discrete state space Markov chains holds over to continuous state space Markov chains in very natural ways with minor modifications. In this short Section, we briefly outline the key differences and in particular the influence this will have on the MCMC discussion which follows.

Firstly, we replace the transition matrix $P$ by a transition kernel $K$. Specifically, if we let $𝕊$ denote the set of values the Markov chain can take ($𝕊$ will often be ${ℝ}^d$ or some subset of it), for all $x,y\in 𝕊$, we have a density $k(x,y)$ such that

$$\int k(x,y)𝑑y=1$$

and for a set $A\subseteq 𝕊$, the probability of moving from state $x$ to a state $y\in A$ is

$${\int }_{A}k(x,y)𝑑y.$$

For a continuous distribution the probability of ever returning to a state $x$ is 0. (For any continuous distribution the probability of observing any single value is 0.) Therefore for recurrence of a continuous state space Markov chain, we talk in terms of returning to a set $A$. We won’t go into the details but for any $x\in 𝕊$, we could define a set of values, $C_x$ say, close to $x$, for example within a certain distance of $x$ and look at the time taken by the Markov chain to return to $x$.

Secondly, we can still talk in terms of a stationary distribution $\pi$, but $\pi$ will now be a probability density function. The stationary distribution $\pi$ satisfies

$$\int \pi (y)k(y,x)𝑑y=\pi (x),$$

where in moving from discrete distributions to continuous distributions we have replaced the summation by integration and the transition matrix by the transition kernel. Similarly, detailed balance is given by

$$\pi (x)k(x,y)=\pi (y)k(y,x).$$

In classical (or frequentist) statistics given a parametric model, we assume that the parameters of the model are fixed, typically unknown, constants. We can then use data to make inference about the parameters, often through the likelihood function. That is, find the maximum likelihood estimate (MLE) of the parameters. We can also find standard errors and confidence intervals for the MLE. See for example the Week 6 notes on the EM algorithm.

Bayesian statistics takes its name from Thomas Bayes. Again we can assume a parametric model for the data. However, rather than assuming that the parameters are fixed constant, we now assume that they are random variables incorporating uncertainty about the parameters. Moreover, we can make use of any prior knowledge we (or others) may have about the parameter values. This is an important feature of Bayesian statistics, the presence of a prior distribution for the parameters. These priors can be either vague or informative, reflecting how much prior knowledge we may have concerning the parameters. The posterior distribution of the parameters is what we are interested in. The posterior distribution for the parameters is dependent upon the parametric model chosen, the data and the prior distribution. From the posterior distribution, we can obtain the modal value for the parameters (these will be very close to the MLE values, however, depending upon the choice of prior distribution these will rarely agree exactly), the mean of the parameters and the variance (covariance) of the parameters. Thus we can obtain a great deal of information from the posterior distribution and there is no need to construct confidence intervals to measure uncertainty in the parameter values as these are obtained directly from the posterior distribution. Note that as the amount of data increases, the effects of the prior will diminish and in the limit as we obtain infinite data, the choice of prior will be inconsequential.

We now consider Bayesian statistics in action and give some motivation of why MCMC methods have been so successful. Consider a parametric model with parameter(s) $\theta$. Let $\pi (\theta )$ denote the pdf (probability density function) of the prior distribution of $\theta$. Let $𝐱=(x_1,x_2,\mathrm{\dots },x_n)$ denote the observed data and let $\pi (𝐱,\theta )$ denote the joint pdf of the data and the parameters. Then by Bayes Theorem,

	$\pi (𝐱,\theta )$	$=$	$\pi (𝐱|\theta )\pi (\theta )$
		$=$	$\pi (\theta |𝐱)\pi (𝐱).$

We are interested in the posterior distribution of $\theta$, that is, the distribution of $\theta$ given the data $𝐱$, ie. $\pi (\theta |𝐱)$. Therefore

$$\pi (\theta |𝐱)=\pi (𝐱|\theta )\pi (\theta )/\pi (𝐱).$$

Since the likelihood $L(\theta ;𝐱)=\pi (𝐱|\theta )$, we have that

$\pi (\theta |𝐱)\propto \text{likelihood}\times \text{prior}.$

(2.2)

(Note that $\pi (𝐱)$ does not depend upon $\theta$.) The likelihood and prior are often relatively easy to obtain, and in such circumstances we can find the posterior distribution up to a constant of proportionality. That is,

$$\pi (\theta |𝐱)=K\times \text{likelihood}\times \text{prior},$$

where $K$ needs to be computed. The computing of $K(=1/\pi (𝐱))$ is often very difficult, if at all analytically possible. This has been a major stumbling block for Bayesian statistics. Fortunately, MCMC techniques can be applied directly to (2.2) to obtain samples from the posterior distribution, thus circumventing the need to compute $K$. In MATH553 you will see alternatives to MCMC for getting around the problem of computing $K$.

Gaussian Example.

This simple example produces a nice analytical answer.

Suppose that we have independent and identically distributed data according to a random variable $X\sim N(\mu ,1)$, where $\mu$ is an unknown parameter. We want to estimate $\mu$.

Prior distribution: Suppose that we know nothing about $\mu$ and take $\pi (\mu )\propto 1$, ie. ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\pi (\mu )𝑑\mu =\mathrm{\infty }$. This is an improper prior since it does not integrate to 1 over the possible values of $\mu$. However, often (but not always) improper prior distributions give rise to proper posterior distributions, so this is not a problem.

Suppose that $x_1,x_2,\mathrm{\dots },x_n$ are independent realisations of $X$.

Then

$$f(x|\mu )=\frac{1}{\sqrt{2\pi }}\exp (-{(x-\mu )}^2/2),$$

and

	$L(\mu ;𝐱)$	$=$	${\displaystyle \prod _{i=1}^n}f(x_i|\mu )$
		$=$	${\displaystyle \prod _{i=1}^n}{\displaystyle \frac{1}{\sqrt{2\pi }}}\exp (-{(x_i-\mu )}^2/2)$
		$=$	${\displaystyle \frac{1}{{(\sqrt{2\pi })}^n}}\exp \left(-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{(x_i-\mu )}^2\right).$

Therefore

	$\pi (\mu |𝐱)$	$\propto$	${\displaystyle \frac{1}{{(\sqrt{2\pi })}^n}}\exp \left(-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{(x_i-\mu )}^2\right)\times 1$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2}}\left({\displaystyle \sum _{i=1}^n}{x}_i^2-2\mu {\displaystyle \sum _{i=1}^n}{x}_i+n{\mu }^2\right)\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2}}\left(n{\mu }^2-2\mu n\overline{x}\right)\right),$

where $\overline{x}=\frac{1}{n}{\sum }_{i=1}^n{x}_i$.

A very useful result. Let $Y$ have probability density function $g(y)$. Then if

$$g(y)\propto \exp \left(-\frac{1}{2}\left(A{y}^2-2By\right)\right),$$

then $Y\sim N(B/A,1/A)$.

Thus

$$\mu |𝐱\sim N(\overline{x},\frac{1}{n}).$$

Therefore the mean of $\mu$ is the sample mean (which is the MLE). Also we have that the variance is $1/n(=\text{Var}\left[X\right]/n)$ which agrees the variance of the MLE. In other words, as the sample size grows, we obtain tighter and tighter estimates of $\mu$.

Conjugate priors

The choice of prior plays an important part in Bayesian statistics. There are two key features, how informative the priors are and conjugacy. We illustrate this idea with an extension of the Gaussian example above to the case where both the mean $\mu$ and variance ${\sigma }^2$ are unknown. The choice of priors, and in particular, conjugacy will be covered in detailed in MATH553.

For $1\le i\le n$, let $X_i|\mu ,{\sigma }^2\sim N(\mu ,{\sigma }^2)$. Suppose also that we have the following prior distributions for $\mu$ and ${\sigma }^2$:

$$\pi (\mu )\sim N({\lambda }_0,1/{\omega }_0)$$

$$\pi ({\sigma }^{-2})\sim \mathrm{Gamma}({\alpha }_0,{\gamma }_0),$$

where $\mu$ and $\sigma$ are assumed to be apriori independent and ${\lambda }_0,{\omega }_0,{\alpha }_0$ and ${\gamma }_0$ are assumed to be known hyperparameters. Hyperparameters are simply the parameters of the prior distributions of the model parameters. Letting $\tau ={\sigma }^{-2}$, ($\tau$ is known as the precision and is the inverse of the variance) we have that, $f(x_i|\tau ,\mu )=\sqrt{\tau /(2\pi )}\exp (-\tau {(x_i-\mu )}^2/2)$, so

$$\pi (\mu ,\tau |𝐱)\propto \prod _{i=1}^n\{{\tau }^{1/2}\exp (-\tau {(x_i-\mu )}^2/2)\}\times \exp (-{\omega }_0{(\mu -{\lambda }_0)}^2/2){\tau }^{{\alpha }_0-1}{e}^{-{\gamma }_0\tau }.$$

Therefore, although the form of the prior distribution is highly tractable, the posterior distribution is a complex two dimensional distribution. This is often the case for multidimensional parameter sets.

However, if we consider the conditional distributions of each of the parameters in turn, these turn out to be rather nice. We exploit the following two observations which are useful for deriving conditional distributions:

1. 1.

   Suppose that the random variable $Y$ has a pdf of the form,

   $$f(y)\propto \prod _{j=1}^m\exp \left(-\frac{{\varphi }_j}{2}{(y-{\delta }_j)}^2\right),$$

   then

   $$Y\sim N(\frac{{\sum }_{j=1}^m{\varphi }_j{\delta }_j}{{\sum }_{j=1}^m{\varphi }_j},\frac{1}{{\sum }_{j=1}^m{\varphi }_j}).$$

   That is, if the pdf of $Y$ is composed as the product of Normal densities then $Y$ has a Normal density with more weight given to those components with larger precisions (smaller variances).

2. 2.

   Suppose that the random variable $W$ has a pdf of the form,

   $$g(w)\propto w^{A-1}\exp (-Bw),$$

   then

   $$W\sim \mathrm{Gamma}(A,B).$$

Firstly, to consider the conditional distribution of $\mu |\tau ,𝐱$ we only need to focus upon terms involving $\mu$. Therefore

	$\pi (\mu |\tau ,𝐱)$	$\propto$	${\displaystyle \prod _{i=1}^n}\left\{\exp \left(-{\displaystyle \frac{\tau }{2}}{(x_i-\mu )}^2\right)\right\}\exp \left(-{\displaystyle \frac{{\omega }_0}{2}}{(\mu -{\lambda }_0)}^2\right)$
		$\sim$	$N({\displaystyle \frac{\tau {\sum }_{i=1}^n{x}_i+{\lambda }_0{\omega }_0}{n\tau +{\omega }_0}},{\displaystyle \frac{1}{n\tau +{\omega }_0}}).$

Similarly, for of $\tau |\mu ,𝐱$ we only need to focus upon terms involving $\tau$, giving

	$\pi (\tau |\mu ,𝐱)$	$\propto$	${\displaystyle \prod _{i=1}^n}\{{\tau }^{1/2}\exp (-\tau {(x_i-\mu )}^2/2)\}\times {\tau }^{{\alpha }_0-1}{e}^{-{\gamma }_0\tau }$
		$=$	${\tau }^{{\alpha }_0+n/2-1}\exp (-\tau ({\displaystyle \frac{1}{2}}i=1^n{(x_i-\mu )}^2+{\gamma }_0))$
		$\sim$	$\mathrm{Gamma}({\alpha }_0+{\displaystyle \frac{n}{2}},{\gamma }_0+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{(x_i-\mu )}^2).$

Thus the conditional distributions of each of the parameters has a nice simple form. This is called conditional conjugacy.

Therefore we will proceed by introducing the Gibbs Sampler which is applicable, where the conditional distributions of each parameter given the other parameters is known but the full joint distribution of the parameters is not known explicitly.

Suppose that we wish to obtain a sample from the multivariate posterior distribution $\pi (𝜽|𝐱)$, where $𝜽=({\theta }_1,{\theta }_2,\mathrm{\dots },{\theta }_d)$ denotes the parameters of the model and $𝐱=(x_1,x_2,\mathrm{\dots },x_N)$ denotes the observed data. The Gibbs sampler does this by successively and repeatedly simulating from the conditional distributions of each component given the other components. This procedure is particularly useful where we have conditional conjugacy, so that the resulting conditional distributions are from standard distributions.

Gibbs Sampler algorithm

1. 1.

   Initialise with ${𝜽}^{(0)}=({\theta }_1^{(0)},\mathrm{\dots },{\theta }_d^{(0)})$.

2. 2.

   For $i=1,2,\mathrm{\dots },n$,

   1. (a)

      Simulate ${\theta }_1^{(i)}$ from the conditional ${\theta }_1|({\theta }_2^{(i-1)},\mathrm{\dots },{\theta }_d^{(i-1)}),𝐱$

   2. (b)

      Simulate ${\theta }_2^{(i)}$ from the conditional ${\theta }_2|({\theta }_1^{(i)},{\theta }_3^{(i-1)},\mathrm{\dots },{\theta }_d^{(i-1)}),𝐱$

   3. (c)

      $\mathrm{\dots }$

   4. (d)

      Simulate ${\theta }_d^{(i)}$ from the conditional ${\theta }_d|({\theta }_1^{(i)},\mathrm{\dots },{\theta }_{d-1}^{(i)}),𝐱$

3. 3.

   Discard the first $k$ iterations and estimate summary statistics of the posterior distribution using ${𝜽}^{(k+1)},{𝜽}^{(k+2)},\mathrm{\dots },{𝜽}^{(n)}$.

How does it work?

Suppose that $𝜽$ comes from $\pi (𝜽|𝐱)$. Then for $i=1,2,\mathrm{\dots },d$, updating the $i^{th}$ component involves drawing a new value of ${\theta }_i$, ${\theta }_i^{\prime }$ from ${\theta }_i|{𝜽}_{-i},𝐱$. Thus ${𝜽}^{\prime }$, the updated set of parameters with ${\theta }_i$ replaced by ${\theta }_i^{\prime }$ also comes from the posterior distribution of $𝜽$.

Thus provided ${𝜽}^{(0)}$ is drawn from the posterior distribution of $𝜽$; ${𝜽}^{(1)},{𝜽}^{(2)},\mathrm{\dots },{𝜽}^{(n)}$ are samples from the posterior distribution. Note that we have a dependent sample.

We have to specify initial values for $𝜽$. If ${𝜽}^{(0)}$ was drawn from the stationary distribution of the Markov chain (distribution we are interested in), then $({\theta }_1^{(0)},\mathrm{\dots },{\theta }_d^{(0)}),\mathrm{\dots },({\theta }_1^{(n)},\mathrm{\dots },{\theta }_d^{(n)})$ would be realisations from the stationary distribution, $\pi (𝜽|𝐱)$. However, we don’t know what the (stationary) distribution is. (If we knew the distribution of $𝜽=({\theta }_1,{\theta }_2,\mathrm{\dots },{\theta }_d)$, there would be no need for MCMC to simulate from it!) Thus typically ${𝜽}^{(0)}$ will not be drawn from the stationary distribution. As $k$ increases ${𝜽}^{(k)}$ increasingly forgets, the initial value ${𝜽}^{(0)}$, and consequently, for large $k$, ${𝜽}^{(k)}$ is approximately from the stationary distribution. Hence, $({\theta }_1^{(k+1)},\mathrm{\dots },{\theta }_d^{(k+1)}),\mathrm{\dots },({\theta }_1^{(n)},\mathrm{\dots },{\theta }_d^{(n)})$ are approximately from the stationary distribution and thus are considered a sample thereof. In some circumstances it is possible to use perfect simulation in which case ${𝜽}^{(0)}$ is drawn from the stationary distribution even though the stationary distribution is unknown. Perfect simulation is not covered in this course and is only viable in some special cases.

We recap the normal example given in the Section 2.4. Suppose that $x_1,x_2,\mathrm{\dots },x_n$ are independent and identically distributed according to $X\sim N(\mu ,1/\tau )$, where $\mu$ and $\tau$ are unknown parameters to be estimated. The following prior distributions are assigned to $\mu$ and $\tau$:-

	$\pi (\mu )$	$\sim$	$N({\lambda }_0,1/{\omega }_0)$
	$\pi (\tau )$	$\sim$	$\mathrm{Gamma}({\alpha }_0,{\gamma }_0).$

We have the following Gibbs sampler algorithm for obtaining a sample of size $N$ from the joint distribution $\theta =(\mu ,\tau )$.

1. 1.

   Initialise with $\theta =({\mu }^{(0)},{\tau }^{(0)})$.
   Any values of ${\mu }_0$ and ${\tau }_0>0$ will be OK. However a reasonable start would be the sample mean and inverse of the sample variance. That is, ${\mu }^{(0)}=\frac{1}{n}{\sum }_{i=1}^n{x}_i$ and ${\tau }^{(0)}=1/\left(\frac{1}{n-1}{\sum }_{i=1}^n{(x_i-{\mu }^{(0)})}^2\right)$.

2. 2.

   For $k=1,2,\mathrm{\dots },N$:-

   1. (a)

      Simulate ${\mu }^{(k)}|{\tau }^{(k-1)}=\tau ,𝐱\sim N(\frac{\tau {\sum }_{i=1}^n{x}_i+{\lambda }_0{\omega }_0}{n\tau +{\omega }_0},\frac{1}{n\tau +{\omega }_0})$

   2. (b)

      Simulate ${\tau }^{(k)}|{\mu }^{(k)}=\mu ,𝐱\sim \mathrm{Gamma}({\alpha }_0+\frac{n}{2},{\gamma }_0+\frac{1}{2}{\sum }_{i=1}^n{(x_i-\mu )}^2)$.

The Gibbs sampler was run to produce a sample of size 1100 from the joint distribution of $\mu$ and $\tau$. We discard the first 100 samples as burn-in using the remaining 1000 samples to estimate $\mu$ and $\sigma (=1/\sqrt{\tau })$.

Unnumbered Figure: Link

Unnumbered Figure: Link

Verification

We will show that $\pi (\theta )$ is indeed the stationary distribution of the Gibbs sampler in the case where $d=2$. The proof for $d>2$ is similar but more long-winded.

Let ${\pi }_1({\theta }_1|{\theta }_2)=\pi ({\theta }_1,{\theta }_2)/\pi ({\theta }_2)$ and let ${\pi }_2({\theta }_2|{\theta }_1)=\pi ({\theta }_1,{\theta }_2)/\pi ({\theta }_1)$, the conditional distributions of ${\theta }_1$ and ${\theta }_2$, respectively. Suppose that ${\theta }^0=({\theta }_1^0,{\theta }_2^0)$ is drawn from $\pi (\cdot )$. Then the pdf of ${\theta }^1$ given ${\theta }^0$ is

$$f({\theta }^1|{\theta }^0)={\pi }_1({\theta }_1^1|{\theta }_2^0){\pi }_2({\theta }_2^1|{\theta }_1^1).$$

We want to show that the pdf of ${\theta }^1$, $f(\cdot )$ is equal to $\pi (\cdot )$.

For any set $A$,

	$P({\theta }^1\in A)$	$=$	${\displaystyle \int 1_{\{{\theta }^1\in A\}}f({\theta }^1)𝑑{\theta }^1}$
		$=$	${\displaystyle \int \int 1_{\{{\theta }^1\in A\}}f({\theta }^1|{\theta }^0)\pi ({\theta }^0)𝑑{\theta }^1𝑑{\theta }^0}$
		$=$	${\displaystyle \int \int \int \int 1_{\{({\theta }_1^1,{\theta }_2^1)\in A\}}{\pi }_1({\theta }_1^1|{\theta }_2^0){\pi }_2({\theta }_2^1|{\theta }_1^1)\pi ({\theta }_1^0,{\theta }_2^0)𝑑{\theta }_1^1𝑑{\theta }_2^1𝑑{\theta }_1^0𝑑{\theta }_2^0}$
		$=$	${\displaystyle \int \int \int 1_{\{({\theta }_1^1,{\theta }_2^1)\in A\}}{\pi }_1({\theta }_1^1|{\theta }_2^0){\pi }_2({\theta }_2^1|{\theta }_1^1)\left[\int \pi ({\theta }_1^0,{\theta }_2^0)𝑑{\theta }_1^0\right]𝑑{\theta }_1^1𝑑{\theta }_2^1𝑑{\theta }_2^0}$
		$=$	${\displaystyle \int \int \int 1_{\{({\theta }_1^1,{\theta }_2^1)\in A\}}{\pi }_1({\theta }_1^1|{\theta }_2^0){\pi }_2({\theta }_2^1|{\theta }_1^1)\pi ({\theta }_2^0)𝑑{\theta }_1^1𝑑{\theta }_2^1𝑑{\theta }_2^0}$
		$=$	${\displaystyle \int \int 1_{\{({\theta }_1^1,{\theta }_2^1)\in A\}}{\pi }_2({\theta }_2^1|{\theta }_1^1)\left[\int {\pi }_1({\theta }_1^1|{\theta }_2^0)\pi ({\theta }_2^0)𝑑{\theta }_2^0\right]𝑑{\theta }_1^1𝑑{\theta }_2^1}$
		$=$	${\displaystyle \int \int 1_{\{({\theta }_1^1,{\theta }_2^1)\in A\}}{\pi }_2({\theta }_2^1|{\theta }_1^1){\pi }_1({\theta }_1^1)𝑑{\theta }_1^1𝑑{\theta }_2^1}$
		$=$	${\displaystyle \int \int 1_{\{({\theta }_1^1,{\theta }_2^1)\in A\}}\pi ({\theta }_1^1,{\theta }_2^1)𝑑{\theta }_1^1𝑑{\theta }_2^1}$
		$=$	${\displaystyle \int 1_{\{{\theta }^1\in A\}}\pi ({\theta }^1)𝑑{\theta }^1}$

as required.

Observations

We make a couple of minor observations about the Gibbs sampler and its output which are very useful in practice.

Firstly, suppose that we interested in the marginal distribution of one or more of the parameters, for example, $\pi ({\theta }_j|𝐱)$. Then ${\theta }_j^{(1)},{\theta }_j^{(2)},\mathrm{\dots },{\theta }_j^{(n)}$ represents a sample from $\pi ({\theta }_j|𝐱)$. That is, we automatically get samples from the marginal distributions of parameters from the Gibbs sampler.

Secondly, whilst it will often be the case that univariate conditional distributions are used to construct the Gibbs sampler, it is possible to update more than one parameter together in a block by utilising multi-dimensional conditional distributions. That is, updating a subset $p$ of the $d$ parameters conditional upon the data and the other $d-p$ parameters. This is particularly useful in regression and similar examples which give rise to multivariate Gaussian distributions as conditional distributions.

We shall use this example to further illustrate the Gibbs sampler in action. The data to be analysed are obtained from Jarrett (1979) and is concerned with the number of coal mining disasters per year, over the period 1851-1962.

A plot of the data is given in the figure below, and suggests that there has been a reduction in the number of disasters per year over the period. The following model was suggested by Carlin et al. (1992):

	$Y_i$	$\sim$	$Po(\theta );i=1,\mathrm{\dots },k;$
	$Y_i$	$\sim$	$Po(\lambda );i=k+1,\mathrm{\dots },n(=112).$

That is, the number of disasters per year is Poisson distributed. Moreover, the first $k$ years have a common mean $\theta$ and the remaining $n-k$ years have a common mean $\lambda$. This is a standard change point problem, in that, there is a behaviour change in the model (here a change in mean) at an unknown point in time.

Unnumbered Figure: Link

We complete the specification of the model by assuming the following prior structure. Let $\theta \sim \mathrm{Gamma}(a_1,b_1)$ and $\lambda \sim \mathrm{Gamma}(a_2,b_2)$, $k$ discrete uniform over $\{1,2,\mathrm{\dots },112\}$, each independent of one another (ie. the parameters are a priori independent) with known hyperparameters $a_1,a_2,b_1$ and $b_2$. More complicated hierarchical models have previously been used to analyse these data, see Carlin et al. (1992).

Therefore

	$L(\theta ,\lambda ,k;𝐲)$	$=$	${\displaystyle \prod _{i=1}^k}{\displaystyle \frac{{\theta }^{y_i}}{y_i!}}{e}^{-\theta }\times {\displaystyle \prod _{i=k+1}^n}{\displaystyle \frac{{\lambda }^{y_i}}{y_i!}}{e}^{-\lambda }$		(2.3)
		$\propto$	${\theta }^{{\sum }_{i=1}^k{y}_i}{e}^{-k\theta }{\lambda }^{{\sum }_{i=k+1}^n{y}_i}{e}^{-(n-k)\lambda }$
		$=$	${\left({\displaystyle \frac{\theta }{\lambda }}\right)}^K{e}^{k(\lambda -\theta )}{\lambda }^{{\sum }_{i=1}^n{y}_i}{e}^{-n\lambda },$		(2.4)

where $K={\sum }_{i=1}^k{y}_i$, and

	$\pi (\theta ,\lambda ,k|𝐲)$	$\propto$	$L(\theta ,\lambda ,k;𝐲)\times \pi (\theta )\times \pi (\lambda )\times \pi (k)$		(2.5)
		$\propto$	$L(\theta ,\lambda ,k;𝐲)\times {\theta }^{a_1-1}{e}^{-b_1\theta }\times {\lambda }^{a_2-1}{e}^{-b_2\lambda }\times {\displaystyle \frac{1}{112}}.$		(2.6)

The conditional distributions of $\theta$, $\lambda$ and $k$, can be obtained from (2.11-2.5) by focussing upon only those terms involving the parameter of interest. Note that $k$ takes discrete values in the range $\{1,2,\mathrm{\dots },112\}$.

The conditional distributions are given as follows:

1. 1.

   Using (2.3) and (2.5), we have that

   	$\pi (\theta |\lambda ,k,𝐲)$	$\propto$	${\theta }^{{\sum }_{i=1}^k{y}_i}{e}^{-k\theta }\times {\theta }^{a_1-1}{e}^{-b_1\theta }$		(2.7)
   		$=$	${\theta }^{{\sum }_{i=1}^k{y}_i+a_1-1}{e}^{-(k+b_1)\theta }.$

   Noting that this is the kernel of a Gamma distribution, see the observation in Section 2.4,

   $\theta |𝐲,\lambda ,k$

   $\sim$

   $\mathrm{Gamma}(a_1+{\displaystyle \sum _{i=1}^k}{y}_i,k+b_1).$

2. 2.

   Using (2.3) and (2.5), we have that

   	$\pi (\lambda |\theta ,k,𝐲)$	$\propto$	${\lambda }^{{\sum }_{i=k+1}^n{y}_i}{e}^{-(n-k)\lambda }\times {\lambda }^{a_2-1}{e}^{-b_2\lambda }$		(2.8)
   		$=$	${\lambda }^{{\sum }_{i=1}^k{y}_i+a_1-1}{e}^{-(k+b_1)\theta }.$

   Noting that this is the kernel of a Gamma distribution, see the observation in Section 2.4,

   $\theta |𝐲,\lambda ,k$

   $\sim$

   $\mathrm{Gamma}(a_2+{\displaystyle \sum _{i=k+1}^n}{y}_i,(n-k)+b_2).$

3. 3.

   Using (2.4) and (2.5), we have that

   $p(k|𝐲,\theta ,\lambda )\propto \exp \{k(\lambda -\theta )\}{(\theta /\lambda )}^K=Q(𝐲;k,\theta ,\lambda ),\text{ say}.$

   (2.9)

   Now $k$ does not have a nice conditional distribution, by which a mean an easily recognisable, well known distribution. However, $k$ is a discrete distribution and there are only 112 possible values. Therefore we can compute the normalising constant directly to give

   $p(k|𝐲,\theta ,\lambda )={\displaystyle \frac{Q(𝐲;k,\theta ,\lambda )}{{\sum }_{l=1}^{112}Q(𝐲;l,\theta ,\lambda )}}.$

   (2.10)

   Sampling from the conditional distribution is straightforward in R using the sample command.

The Gibbs sampler was applied to the coal mining data set and run for 1100 iterations to obtain samples from $(\theta ,\lambda ,k)$. We set $a_1=a_2=b_1=b_2=1$ for the priors. The results of the Gibbs sampler are shown in the figures below.

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

Convergence of the algorithm seems to be rapid (after compensating for a poor choice of starting values, namely, $k=100$). Therefore I deleted the first 100 values and based all the analysis on the remaining 1000 iterations. In other words, I took the first 100 iterations as burn-in.

A histogram of the posterior distribution of year $(k+1850)$ is given below. The posterior mode is at $k=41$, and there is very strong evidence to suggest that the change point occurs about the posterior mode which corresponds to the year 1891.

Unnumbered Figure: Link

It is clear from the trace plots above and kernel density estimates for the posterior distributions that $\lambda$ is (almost certainly) less than $\theta$. This corresponds to less disasters, and hence, a safer working environment after the change point. Note that any summary statistic or distributional quantity related to the joint distribution of $\lambda$ and $\theta$ such as the distribution $\theta -\lambda$ can be estimated using the output from the Gibbs sampler. This is simply done by looking at each pair of realisations ${(\theta ,\lambda )}_i$ giving ${(\theta -\lambda )}_i$ as the $i^{th}$ realisation from the posterior distribution of the difference.

Unnumbered Figure: Link

Unnumbered Figure: Link

Unnumbered Figure: Link

We can easily sample from the predictive distribution. Let $z$ denote a future observation. Future observations are iid according to $\mathrm{Po}(\lambda )$ and does not depend upon $𝐲$, $k$ or $\theta$ given $\lambda$. Therefore the predictive distribution of $z$ given $𝐲$, $\pi (z|𝐲)$ satisfies

$$\pi (z|𝐲)=\int \pi (z|\lambda )\pi (\lambda |𝐲)𝑑\lambda .$$

Samples from $\pi (\lambda |𝐲)$ are given by the Gibbs sampler in the form of ${\lambda }^{(1)},{\lambda }^{(2)},\mathrm{\dots },{\lambda }^{(N)}$. Therefore we can estimate $\pi (z|𝐲)$ using the Monte-Carlo estimate

	$\widehat{\pi (z|𝐲)}$	$=$	${\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\pi (z|{\lambda }^{(i)})$		(2.11)
		$=$	${\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{{({\lambda }^{(i)})}^z}{z!}}\exp (-{\lambda }^{(i)}),$

or simply simulate $z_i|{\lambda }^{(i)}\sim \mathrm{Po}({\lambda }^{(i)})$ and estimate $\pi (z|𝐲)$ by

$\widehat{p(z|𝐲)}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}{1}_{\{z_i=z\}}.$

(2.12)

The estimate in (2.11) has smaller variance, ie. is better. However, (2.12) is easier to compute and only requires the ability to simulate from the distribution.

For the coal mining data set, we can use the estimator to predict the number of disasters in future years (Poisson with rate $\lambda$). A graph of the predictive distribution alongside the corresponding estimates based on the posterior mean ($\pi (z|\overline{\lambda })$) are given for a sample of 1000 realisations in the figure below. Since there is limited variation in the posterior distribution of $\lambda$, the two distributions are virtually identical. However, in general the predictive distribution will demonstrate greater variability owing to the uncertainty in $\lambda$.

Unnumbered Figure: Link

Predictive and estimative (based on $\pi (z|\overline{\lambda })$) distributions of number of coal mining disasters per year (solid line is predictive)

R code for applying the Gibbs sampler to generate samples from a bivariate normal distribution and the outline of a Gibbs sampler algorithm for the coal mining data.

1. 1.

   Bivariate Normal
   The program biv (in biv normR.R) generates values, using the Gibbs sampler, from the bivariate normal distribution:-

   Can you modify the code to produce samples from the bivariate normal distribution distribution:-

   given that

   	$Y_1|Y_2=y_2$	$\sim$	$N({\mu }_1+{\sigma }_1\rho (y_2-{\mu }_2)/{\sigma }_2,(1-{\rho }^2){\sigma }_1^2)$
   	$Y_2|Y_1=y_1$	$\sim$	$N({\mu }_2+{\sigma }_2\rho (y_1-{\mu }_1)/{\sigma }_1,(1-{\rho }^2){\sigma }_2^2).$

2. 2.

   Coal Mining Example
   Complete the Gibbs sampling algorithm for the coal mining data (in coal outlineR.R).

   Remember:-

   	$\theta |𝐲,\lambda ,k$	$\sim$	$Gamma(a_1+{\displaystyle \sum _{i=1}^k}{y}_i,k+b_1)$
   	$\lambda |𝐲,\theta ,k$	$\sim$	$Gamma(a_2+{\displaystyle \sum _{i=k+1}^n}{y}_i,n-k+b_2)$
   	$p(k|𝐲,\theta ,\lambda )$	$=$	$Q(𝐲;k,\theta ,\lambda )/{\displaystyle \sum _{j=1}^n}Q(𝐲;j,\theta ,\lambda ),$

   where

   $$Q(𝐲;k,\theta ,\lambda )=\exp \{k(\lambda -\theta )\}{(\theta /\lambda )}^K.$$

   Hint: To sample $k$ use the sample command in R.

   Apply the algorithm to the coal mining data.