4 Week 9: The Metropolis-Hastings Algorithm

4.1 Introduction to the Metropolis-Hastings algorithm

The Gibbs sampler is a special case of the general Metropolis-Hastings algorithm. All MCMC algorithms are based upon the Metropolis-Hastings algorithm, but as with the Gibbs sampler other special cases are considered in their own right.

Suppose that we wish to simulate from a (multivariate) distribution $f(\theta)$ , which we term the target distribution. We introduce an arbitrary but specified transition probability $q(\theta,\phi)$ from which simulation is straightforward. The transition probability $q(\theta,\phi)$ is the probability density of moving from $\theta$ to $\phi$ .

If $\theta_{0}$ is drawn from $f(\cdot)$ , the Metropolis-Hastings will generate (dependent) samples $\theta_{0},\theta_{1},\ldots$ , from $f(\cdot)$ .

The Metropolis-Hastings algorithm is:

Metropolis-Hastings algorithm

1.

Given the current position $\theta_{n}=\vartheta$ , generate a ‘candidate value’, $\phi^{\ast}$ from $q(\vartheta,\phi)$ .
Note that $q(\cdot,\cdot)$ is often referred to as the proposal distribution, in that, we are proposing a candidate value $\phi^{\ast}$ using it.
2.

Calculate the ‘acceptance probability’, $\alpha(\vartheta,\phi^{\ast})$ , given by

$\alpha(\vartheta,\phi)=\min\left\{\frac{f(\phi)q(\phi,\vartheta)}{f(\vartheta)% q(\vartheta,\phi)},1\right\}.$
3.

With probability $\alpha(\vartheta,\phi^{\ast})$ accept the candidate value and set $\theta_{n+1}=\phi^{\ast}$ ; otherwise reject the candidate value and set $\theta_{n+1}=\vartheta$ .
4.

Repeat until a sample of the desired size is obtained.

It is not difficult to show the $f(\theta)$ is the stationary distribution of this chain. This is done below. Under mild conditions, including making sure that the proposal distribution does not lead to a reducible Markov chain, the chain will converge to this stationary distribution.

To show that the stationary distribution of the Metropolis-Hastings algorithm is indeed $\pi(\cdot)$ , it is necessary to show that, for all $\theta$

\int f(\phi)g(\phi,\vartheta)\,d\phi=f(\vartheta)

where $g(\phi,\vartheta)$ is the probability density, for $\theta_{1}=\vartheta$ given that $\theta_{0}=\phi$ . For $\phi\neq\vartheta$ , the density $g(\phi,\vartheta)$ consists of two parts, first we need to propose the move from $\phi$ to $\vartheta$ , ie. $q(\phi,\vartheta)$ and secondly, we need to accept the move, $\alpha(\phi,\vartheta)$ . Thus $g(\phi,\vartheta)=q(\phi,\vartheta)\alpha(\phi,\vartheta)$ . However, the probability $g(\vartheta,\vartheta)$ is slightly more involved. (Note that typically there is a non-zero probability of staying at $\vartheta$ .) Firstly, we have the probability of proposing that we stay where we are $q(\vartheta,\vartheta)$ which is a move we will always accept! Secondly, we could propose a candidate value $\varphi$ and reject such a move, meaning that we would stay where we are. Thus

g(\vartheta,\vartheta)=q(\vartheta,\vartheta)\alpha(\vartheta,\vartheta)+\int q% (\vartheta,\varphi)(1-\alpha(\vartheta,\varphi))\,d\varphi.

In general it is not a good idea to have $q(\vartheta,\vartheta)>0$ , proposing to stay where you are just slows down the process. Therefore without loss of generality, we take $q(\vartheta,\vartheta)=0$ .

Therefore bringing all of the above arguments together, we have that

$\displaystyle\int f(\phi)g(\phi,\vartheta)\,d\phi$	$\displaystyle=$	$\displaystyle\int f(\phi)q(\phi,\vartheta)\alpha(\phi,\vartheta)\,d\phi+\int f% (\vartheta)q(\vartheta,\varphi)(1-\alpha(\vartheta,\varphi))\,d\varphi$
	$\displaystyle=$	$\displaystyle\int f(\phi)q(\phi,\vartheta)\min\left\{1,\frac{f(\vartheta)q(% \vartheta,\phi)}{f(\phi)q(\phi,\vartheta)}\right\}\,d\phi+\int f(\vartheta)q(% \vartheta,\varphi)\left(1-\min\left\{1,\frac{f(\varphi)q(\varphi,\vartheta)}{f% (\vartheta)q(\vartheta,\varphi)}\right\}\right)\,d\varphi$
	$\displaystyle=$	$\displaystyle\int\min\{f(\phi)q(\phi,\vartheta),f(\vartheta)q(\vartheta,\phi)% \}\,d\phi+\int f(\vartheta)q(\vartheta,\varphi)\,d\varphi-\int\min\{f(\varphi)% q(\varphi,\vartheta),f(\vartheta)q(\vartheta,\varphi)\}\,d\varphi$
	$\displaystyle=$	$\displaystyle\int\min\{f(\phi)q(\phi,\vartheta),f(\vartheta)q(\vartheta,\phi)% \}\,d\phi+f(\vartheta)\int q(\vartheta,\varphi)\,d\varphi-\int\min\{f(\phi)q(% \phi,\vartheta),f(\vartheta)q(\vartheta,\phi)\}\,d\phi$
	$\displaystyle=$	$\displaystyle f(\vartheta),$

as required. This is since $q(\vartheta,\cdot)$ is a probability density, and hence, integrates to 1 and the other two terms cancel with each other.

There exists a simpler proof based upon detailed balance. $f(\cdot)$ is the stationary distribution of a Markov chain with transition kernel $g(\phi,\theta)$ if for all $\phi,\theta$ ,

f(\phi)g(\phi,\theta)=f(\theta)g(\theta,\phi).

Specifically for the Metropolis-Hastings algorithm, $g(\phi,\theta)=q(\phi,\theta)\alpha(\phi,\theta)$ , and therefore,

$\displaystyle f(\phi)g(\phi,\theta)$	$\displaystyle=$	$\displaystyle f(\phi)q(\phi,\theta)\min\left\{1,\frac{f(\theta)q(\theta,\phi)}% {f(\phi)q(\phi,\theta)}\right\}$
	$\displaystyle=$	$\displaystyle\min\left\{f(\phi)q(\phi,\theta),f(\theta)q(\theta,\phi)\right\}$
	$\displaystyle=$	$\displaystyle f(\theta)q(\theta,\phi)\min\left\{\frac{f(\phi)q(\phi,\theta)}{f% (\theta)q(\theta,\phi)},1\right\}$
	$\displaystyle=$	$\displaystyle f(\theta)g(\theta,\phi),$

as required.

There are (infinitely) many possible choices for the proposal distribution $q$ . Two of the most common choices are outlined below.

1.

Independence sampler
For all $\vartheta$ , let $q(\vartheta,\phi)=q(\phi)$ for some probability density $q(\cdot)$ , that is, the proposal density is independent of the current value $\vartheta$ . This idea is similar to rejection sampling, in that we sample from one probability density $q(\cdot)$ and then accept or reject the samples as coming from the probability density $f(\cdot)$ . This method works well when we can find a simple (well understood, easy to simulate from) probability density $q(\cdot)$ which closely approximates $f(\cdot)$ . Note that in this case the acceptance probability is given by

$\alpha(\vartheta,\phi)=\min\left\{\frac{w(\phi)}{w(\vartheta)},1\right\},$

where $w(\varphi)=f(\varphi)/q(\varphi)$ denotes the relative weights of the two densities at $\varphi$ .

Example. Mixture of normal distributions as the target (stationary) distribution, $X$ , with a normal distribution for the proposal distribution, $Y$ .

$\displaystyle f(x)$ $\displaystyle=$ $\displaystyle\frac{1}{2}\left\{\frac{1}{\sqrt{2\pi}}\exp(-(x+1)^{2}/2)+\frac{1% }{2\sqrt{2\pi}}\exp(-x^{2}/8)\right\}\hskip 14.226378pt\mbox{(black)}$

$\displaystyle q(y)$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2\pi}\sqrt{5}}\exp(-(y+0.5)^{2}/10)\hskip 14.22637% 8pt\mbox{(red)}$

That is,

$\displaystyle X$ $\displaystyle\sim$ $\displaystyle\left\{\begin{array}[]{ll}N(-1,1)&w.p.~{}0.5\\ N(0,4)&w.p.~{}0.5\end{array}\right.$

$\displaystyle Y$ $\displaystyle\sim$ $\displaystyle N(-0.5,5)$

Unnumbered Figure: Link
2.

Random-walk Metropolis
Suppose that we are currently at $\vartheta$ . In stationarity it is likely that $\vartheta$ has high posterior probability. (Note that in stationarity the chances of being at $\vartheta$ is proportional to $f(\vartheta)$ .) Therefore it might seem sensible to propose a candidate value, $\phi$ from a distribution which is centred about the current value $\vartheta$ . For example, $\phi|\vartheta\sim N(\vartheta,\sigma^{2})$ , that is, we propose $\phi$ from a normal distribution with mean $\vartheta$ . If the proposal distribution is symmetric about $\vartheta$ as in the case of the normal distribution, then we have the nice property that

$q(\vartheta,\phi)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2\sigma^{2}}(% \phi-\vartheta)^{2}\right)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2% \sigma^{2}}(\vartheta-\phi)^{2}\right)=q(\phi,\vartheta).$

Hence in this case, the acceptance probability $\alpha(\vartheta,\phi)$ , simplifies to

$\alpha(\vartheta,\phi)=\min\left\{\frac{f(\phi)}{f(\vartheta)},1\right\}.$

Example. Mixture of normal distributions is the target with a normal distribution for the proposal.

$\displaystyle f(x)$ $\displaystyle=$ $\displaystyle\frac{1}{2}\left\{\frac{1}{\sqrt{2\pi}}\exp(-(x+1)^{2}/2)+\frac{1% }{2\sqrt{2\pi}}\exp(-x^{2}/8)\right\}\hskip 14.226378pt\mbox{(black)}$

$\displaystyle q(y|x)$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2\pi}}\exp(-(y-x)^{2}/2)\hskip 14.226378pt\mbox{(% red)}$

That is,

$\displaystyle X$ $\displaystyle\sim$ $\displaystyle\left\{\begin{array}[]{ll}N(-1,1)&w.p.~{}0.5\\ N(0,4)&w.p.~{}0.5\end{array}\right.$

$\displaystyle Y$ $\displaystyle\sim$ $\displaystyle N(x,1)$

For $x=1$ (marked with a vertical green line),

Unnumbered Figure: Link

(Note that the proposal density has been divided by 2 in the above figure.)

Why choose $\alpha(\vartheta,\phi)=\min\left\{\frac{f(\phi)q(\phi,\vartheta)}{f(\vartheta)% q(\vartheta,\phi)},1\right\}$ ?
The answer is given by Peskun (1973) and is as follows. $\alpha(\vartheta,\phi)$ is the largest value that the acceptance probability can take for the given choice of $q(\cdot,\cdot)$ which ensures that the stationary distribution of the Markov chain is given by $f(\cdot)$ . Choosing $\alpha(\cdot,\cdot)$ means the Markov chain moves around the sample space as fast as possible (given $q(\cdot,\cdot)$ ) and therefore converges as quickly as possible.

Accepting a move.
The simplest way to accept a proposed move computationally is to generate a random number $u$ from $U(0,1)$ and accept the move if $u\leq\alpha(\vartheta,\phi)$ (reject otherwise). This leads us to a couple of observations which can assist in implementing MCMC.

1.

Let $R(\vartheta,\phi)=\frac{f(\phi)q(\phi,\vartheta)}{f(\vartheta)q(\vartheta,\phi)}$ , so that $\alpha(\vartheta,\phi)=\min\{1,R(\vartheta,\phi)\}$ . It suffices to accept the move if $u\leq R(\vartheta,\phi)$ (reject otherwise), since if $R(\vartheta,\phi)\neq\alpha(\vartheta,\phi)$ , $R(\vartheta,\phi)>1$ and the proposed move is accepted anyway.
2.

It suffices to accept the move if $\log(u)\leq\log(R(\vartheta,\phi))$ (reject otherwise) as $\log$ is a monotonic function. This can help with numerical stability, although this unlikely to be the case for problems we consider in this course.

4.2 The Metropolis-Hastings algorithm in a Bayesian context

Thus far we have described the Metropolis-Hastings algorithm for obtaining samples from a probability distribution with probability density function $f(\cdot)$ . Typically, $f(\cdot)$ is the posterior distribution for a Bayesian statistical problem. That is, $f(\theta)=\pi(\theta|\mathbf{x})$ . The key to the success of the Metropolis-Hastings algorithm (MCMC) is that we only need to know $\pi(\theta|\mathbf{x})$ up to a constant of proportionality. To see this note that

\pi(\theta|\mathbf{x})=\pi(\mathbf{x}|\theta)\pi(\theta)/\pi(\mathbf{x}).

Then the acceptance probability of the Metropolis-Hastings algorithm can be written as

$\displaystyle\alpha(\vartheta,\phi)$	$\displaystyle=$	$\displaystyle\min\left\{\frac{\pi(\phi\|\mathbf{x})q(\phi,\vartheta)}{\pi(% \vartheta\|\mathbf{x})q(\vartheta,\phi)},1\right\}$
	$\displaystyle=$	$\displaystyle\min\left\{\frac{\{\pi(\mathbf{x}\|\phi)\pi(\phi)/\pi(\mathbf{x})% \}q(\phi,\vartheta)}{\{\pi(\mathbf{x}\|\vartheta)\pi(\vartheta)/\pi(\mathbf{x})% \}q(\vartheta,\phi)},1\right\}$
	$\displaystyle=$	$\displaystyle\min\left\{\frac{\pi(\mathbf{x}\|\phi)\pi(\phi)q(\phi,\vartheta)}{% \pi(\mathbf{x}\|\vartheta)\pi(\vartheta)q(\vartheta,\phi)},1\right\}.$

That is, it is not necessary to compute the marginal likelihood, the constant of proportionality of the posterior distribution.

The general Metropolis-Hastings algorithm has the major advantage over the Gibbs sampler in that we do not need to know any of the conditional distributions. All that is required is to be able to simulate from $q$ which we can choose arbitrarily subject to a few basic constraints. Moreover, and this can be of crucial importance, we only need to know $\pi(\cdot|\mathbf{x})$ up to proportionality. This is because we only need to consider the ratio of the posterior density at different values, and hence, the proportionality constants cancel. For example, with the independence sampler we only need to know $w$ up to proportionality. Whilst, we have great freedom in choosing $q$ , it is important to make good choices concerning $q$ . A poor choice of $q$ can seriously affect the efficiency and convergence of the Metropolis-Hastings algorithm.

It is often desirable to create hybrid MCMC chains which are hybrid of different types of MCMC chains. A particularly useful application is when using the Gibbs sampler. It is often impossible to obtain the conditional distributions of all the parameters. In this case, it is natural to use a one-dimensional Metropolis-Hastings step to simulate from the one-dimensional conditionals. For example, use the independence sampler or random walk Metropolis for one of the parameters which does not have a nice conditional distribution. Then the accepted value will be from the correct one-dimensional conditional distribution but we have not needed to simulate directly from the conditional distribution. This type of algorithm is sometimes referred to as Metropolis-within-Gibbs or componentwise MCMC.

4.2.1 Componentwise MCMC

Suppose that we want to obtain samples from $\pi(\mbox{\boldmath$\theta$}|\mathbf{x})=\pi(\mathbf{x}|\mbox{\boldmath$\theta% $})\pi(\mbox{\boldmath$\theta$})/\pi(\mathbf{x})$ (a posterior distribution), where $\mbox{\boldmath$\theta$}=(\theta_{1},\ldots,\theta_{d})$ are the parameters of interest and $\mathbf{x}$ denotes data. Componentwise MCMC follows the same procedure as the Gibbs sampler by updating subsets of the parameters, often one parameter, sequentially keeping the other parameters fixed.

That is, within each iteration, for $k=1,2,\ldots,d$ .

1.

Propose $\theta_{k}^{\prime}$ according to some proposal $q_{k}(\mbox{\boldmath$\theta$},\theta_{k}^{\prime})$ , where $\theta$ denotes the current value of the Markov chain.
Note that the proposal distribution $q_{k}(\cdot,\cdot)$ could depend on the whole of $\theta$ but often it will only depend on $\theta_{k}$ .
2.

Set $\mbox{\boldmath$\theta$}^{\prime}=(\theta_{1},\ldots,\theta_{k-1},\theta_{k}^{% \prime},\theta_{k+1},\ldots\theta_{d})$ . (Only the $k^{th}$ element differs from $\theta$ .)

Compute the acceptance probability $\alpha_{k}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{\prime})$ for the proposed move from $\theta$ to $\mbox{\boldmath$\theta$}^{\prime}$ .

	$\displaystyle\alpha_{k}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{% \prime})$	$\displaystyle=$	$\displaystyle\min\left\{1,\frac{\pi(\mbox{\boldmath$\theta$}^{\prime}\|\mathbf{% x})}{\pi(\mbox{\boldmath$\theta$}\|\mathbf{x})}\frac{q_{k}(\mbox{\boldmath$% \theta$},\theta_{k}^{\prime})}{q_{k}(\mbox{\boldmath$\theta$}^{\prime},\theta_% {k})}\right\}$
		$\displaystyle=$	$\displaystyle\min\left\{1,\frac{\pi(\mathbf{x}\|\mbox{\boldmath$\theta$}^{% \prime})\pi(\mbox{\boldmath$\theta$}^{\prime})}{\pi(\mathbf{x}\|\mbox{\boldmath% $\theta$})\pi(\mbox{\boldmath$\theta$})}\frac{q_{k}(\mbox{\boldmath$\theta$},% \theta_{k}^{\prime})}{q_{k}(\mbox{\boldmath$\theta$}^{\prime},\theta_{k})}% \right\}.$

4.

With probability $\alpha_{k}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{\prime})$ , accept the proposed move set $\mbox{\boldmath$\theta$}=\mbox{\boldmath$\theta$}^{\prime}$ . Otherwise keep $\theta$ unchanged (reject the proposed move).

The above algorithm can be adapted to update a subset of the parameters at once. In other words, componentwise MCMC is not restricted to single parameter updates.

4.3 Independence sampler

4.3.1 Introduction

The ideal proposal distribution for an independence sampler is $q(\phi)=f(\phi)$ .
Why? The acceptance rate is always 1 and we get iid samples from $f(\cdot)$ .

Why is this of no direct use?
If we can sample from $f(\cdot)$ directly there is no need for MCMC.

How might it, nonetheless, help us to choose a sensible proposal distribution?
We can look for $q(\cdot)$ which is a good approximation for $f(\cdot)$ .

What might a good acceptance rate for the independence sampler be?
As close to 1 as possible whilst exploring the whole of the posterior distribution.

4.3.2 Genetics example (revisited)

To illustrate the independence sampler, and to illustrate the effect of the choice of $q$ on efficiency, we return to the genetics example, analysed using the Gibbs sampler in Week 8. Note that using the Metropolis-Hastings algorithm, we can use

(f(\theta)=)\pi(\theta|\mathbf{x})\propto(2+\theta)^{x_{1}}(1-\theta)^{x_{2}+x% _{3}}\theta^{x_{4}},

directly.

Remember that there are 197 animals divided into four categories. The number of animals in each category are

\mathbf{x}=(x_{1},x_{2},x_{3},x_{4})=(125,18,20,34),

with category probabilities

(p_{1},p_{2},p_{3},p_{4})=\left(\frac{1}{2}+\frac{\theta}{4},\frac{1}{4}(1-% \theta),\frac{1}{4}(1-\theta),\frac{\theta}{4}\right).

For the proposal $q(\theta)$ we shall use $\theta\sim{\rm Beta}(a,b)$ and consider efficiency for different choices of $a$ and $b$ . This gives

w(\theta)=C\frac{(2+\theta)^{x_{1}}(1-\theta)^{x_{2}+x_{3}}\theta^{x_{4}}}{% \theta^{a-1}(1-\theta)^{b-1}},

where $C$ is a constant independent of $\theta$ which is not important for the MCMC.

We shall consider three choices for $(a,b)$ , the uniform $(1,1)$ , $(6,4)$ and $(1,4)$ . We can see from the figures below that $(6,4)$ works best. This is because the mean of the proposal $0.6$ is close to the mean of the posterior distribution $0.624$ , and the shape of the density $q(\cdot)$ is the best approximation of $\pi(\cdot|\mathbf{x})$ , see the plot comparing the density of $q(\cdot)$ (solid line) with $\pi(\cdot|\mathbf{x})$ (dashed line) below. Even better results are obtained by using say $(31,19)$ . However, there are a number of reasons for preferring a proposal distribution $q(\cdot)$ which is over-dispersed relative to the posterior distribution $\pi(\cdot|\mathbf{x})$ . Roughly speaking the proposal distribution has a larger variance than the posterior distribution. Note that at the other extreme, $(a,b)=(1,4)$ is very inefficient since the proposal distribution does not reflect the posterior distribution $\pi(\cdot|\mathbf{x})$ , ie. very few of the proposed moves are from high density regions of the posterior distribution. Finally, note that the naive uniform proposal performs adequately.

sigene1.pdf

Unnumbered Figure: Link

In the final figure the solid line denotes ${\rm Beta}(6,4)$ distribution and the dashed line denotes $\pi$ .

4.3.3 Why should $q$ be over-dispersed relative to $f$ ?

The closer $q$ approximates $f$ , the better. Thus we look to have $w(\theta)\approx 1$ for all $\theta$ . However, the theoretical (and practical) performance of the MCMC algorithm relies upon $\sup_{\theta}w(\theta)$ . In particular, there are convergence problems if $\sup_{\theta}w(\theta)=\infty$ . In fact the algorithm is reducible if $\sup_{\theta}w(\theta)=\infty$ . There are no such problems caused by having $w(\theta)\approx 0$ , although the efficiency of the algorithm can be affected. The main cause of $\sup_{\theta}w(\theta)=\infty$ is choosing $q$ to have too small a variance (under-dispersed relative to $f$ ).

Suppose that $\sup_{\theta}w(\theta)=K$ . Then the probability of accepting a move away from $\theta^{\ast}$ , where $w(\theta^{\ast})=K$ is

$\displaystyle\int_{-\infty}^{\infty}\left(1\wedge\frac{w(\phi)}{w(\theta^{\ast% })}\right)q(\phi)\,d\phi$	$\displaystyle=$	$\displaystyle\int_{-\infty}^{\infty}\frac{w(\phi)}{K}q(\phi)\,d\phi$
	$\displaystyle=$	$\displaystyle\int_{-\infty}^{\infty}\frac{f(\phi)}{K}\,d\phi$
	$\displaystyle=$	$\displaystyle\frac{1}{K}.$

Therefore if the independence sampler is started at $\theta^{\ast}$ , it takes a geometrically distributed number of iterations to move from $\theta^{\ast}$ with success probability $1/K$ . As $K\rightarrow\infty$ , the probability of accepting a move away from $\theta^{\ast}$ converges to 0 which is clearly problematic. It is worth noting that for any point $\vartheta$ , the number of iterations the independence sampler spends at the point $\vartheta$ before moving is geometrically distributed with success probability (of moving)

\int_{-\infty}^{\infty}\left(1\wedge\frac{w(\phi)}{w(\vartheta)}\right)q(\phi)% \,d\phi.

A measure of the efficiency of an independent sampler is the average acceptance rate. That is,

			$\displaystyle\int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty}\left(1\wedge% \frac{w(\phi)}{w(\vartheta)}\right)q(\phi)\,d\phi\right]f(\vartheta)\,d\vartheta$
		$\displaystyle=$	$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\left(f(\phi)q(% \vartheta)\wedge f(\vartheta)q(\phi)\right)\,d\phi\,d\vartheta.$

The acceptance rate is high, when $f$ and $q$ are in close agreement. Furthermore, if $q(\cdot)=f(\cdot)$ , then the acceptance rate is

\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\vartheta)f(\phi)\,d\vartheta% \,d\phi=1.

4.4 Independence sampler within Gibbs: Censored data

We consider a censored data example. The approach has similarities to that taken in Section 1.6 for analysing censored data using the EM algorithm, in that again we use data augmentation.

Suppose that $x_{1},x_{2},\ldots,x_{n}$ are independent and identically distributed according to $f(x|\theta)$ . However, assume that there is right censoring at the value $a$ . That is, for all observations greater than $a$ , the only information we have is that the observation exceeds $a$ . Therefore, we have a sample $y_{1},y_{2},\ldots,y_{n}$ , where $y_{i}=\min\{x_{i},a\}$ . The common notation is to denote censored values with an asterisk superscript. Usually for convenience we rearrange the sample as

y_{1},\ldots,y_{m},y_{m+1}^{\ast},\ldots,y_{n}^{\ast},

where $y_{j}^{\ast}=a$ . Then with $\mathbf{y}=(y_{1},y_{2},\ldots,y_{n})$ denoting the observed data, we have that

L(\theta|\mathbf{y})=\prod_{i=1}^{m}f(y_{i}|\theta)\{1-F_{\theta}(a)\}^{n-m},

where $F_{\theta}$ is the cdf corresponding to $f$ . Typically it is difficult to work with $L(\theta|\mathbf{y})$ because of the presence of the non-standard term $\{1-F_{\theta}(a)\}^{n-m}$ . Let $\mathbf{z}=(z_{m+1},z_{m+2},\ldots,z_{n})$ represent the true but unobserved values of $y_{m+1}^{\ast},\ldots,y_{n}^{\ast}$ . Then the pair $(\mathbf{y},\mathbf{z})$ corresponds to complete information and

	$\displaystyle L(\theta;\mathbf{y},\mathbf{z})$	$\displaystyle=$	$\displaystyle\pi(\mathbf{y},\mathbf{z}\|\theta)$
		$\displaystyle=$	$\displaystyle\prod_{i=1}^{m}f(y_{i}\|\theta)\prod_{i=m+1}^{n}f(z_{i}\|\theta),$

which is generally straightforward to work with. (Note the obvious similarities with the EM algorithm in Week 6.)

Suppose that $X\sim{\rm Gamma}(2,\delta)$ , where $\delta$ is an unknown with a ${\rm Gamma}(\alpha,\beta)$ prior on $\delta$ . Suppose that we have the following 20 observations with the data censored at $a=2.0$ :

	$\displaystyle 0.7596,\;1.5408,\;2.0000^{\ast},\;0.7261,\;0.2157,\;1.1026,\;2.0% 000^{\ast}$
	$\displaystyle 2.0000^{\ast},\;1.3579,\;0.9520,\;1.2688,\;2.0000^{\ast},\;2.000% 0^{\ast},\;2.0000^{\ast},$
	$\displaystyle 1.5395,\;1.8802,\;0.8471,\;0.4374,\;1.5395,\;1.2576,$

where $\ast$ denotes censored observations.

Note that cumulative distribution function of $X$ is

F_{X}(x)=1-e^{-\delta x}(1+\delta x)\;\;\;(x>0),

whilst the pdf of $X$ is $f(x)=f(x)=\delta^{2}xe^{-\delta x}$ $(x>0)$ .

Therefore if we reorder the data into ascending order, we have that

L(\delta;\mathbf{y})=\prod_{i=1}^{14}\delta^{2}y_{i}e^{-\delta y_{i}}\left\{e^% {-2\delta}(1+2\delta)\right\}^{6},

and it is difficult to write down the posterior distribution of $\delta$ . We could proceed as above and use an independence sampler to obtain samples from $\delta$ , possibly using a Gamma proposal distribution. In other words, take an approach similar to that used for the genetics example above.

However letting $\mathbf{x}=(\mathbf{y},\mathbf{z})$ ,

	$\displaystyle L(\delta;\mathbf{x})$	$\displaystyle=$	$\displaystyle\prod_{i=1}^{20}\delta^{2}x_{i}e^{-\delta x_{i}}$
		$\displaystyle\propto$	$\displaystyle\delta^{40}\exp\left(-\delta\sum_{i=1}^{20}x_{i}\right).$

Hence

	$\displaystyle\pi(\delta\|\mathbf{x})$	$\displaystyle\propto$	$\displaystyle\delta^{40}\exp\left(-\delta\sum_{i=1}^{20}x_{i}\right)\times% \delta^{\alpha-1}\exp(-\delta\beta)$
		$\displaystyle\propto$	$\displaystyle\delta^{40+\alpha-1}\exp\left(-\delta\left(\sum_{i=1}^{20}x_{i}+% \beta\right)\right),$

giving $\delta|\mathbf{y},\mathbf{z}\sim{\rm Gamma}(40+\alpha,\sum_{i=1}^{20}x_{i}+\beta)$ .

For $i=1,2,\ldots,m$ , $y_{i}$ is observed, so no updating is required. For $i=m+1,m+2,\ldots,n$ , $z_{i}|y_{i}=a^{\ast},\delta\sim X^{a}$ , where $X^{a}$ is distributed according to ${\rm Gamma}(2,\delta)$ conditioned to be greater than $a$ . Simulating values from $X^{a}$ directly is not possible and the following independence sampler is proposed. Propose $z_{i}^{\prime}$ from $q(z)=24/z^{4}$ $(z>2)$ and $q(z)=0$ otherwise. Simulation from this distribution can be done using inversion of the cdf. That is, simulate $u_{i}\sim(0,1)$ and set $z_{i}^{\prime}=(8/u_{i})^{1/3}$ . Thus, if $f_{a}(z)=f(z)/P(X>a)$ $(z>a)$ is the pdf of $X^{a}$ , for $z>a=2$ ,

$\displaystyle w(z)$	$\displaystyle=$	$\displaystyle\frac{f_{a}(z)}{q(z)}$
	$\displaystyle=$	$\displaystyle\frac{z\exp(-\delta z)/P(X>2\|\delta)}{24/z^{4}}$
	$\displaystyle\propto$	$\displaystyle z^{5}\exp(-\delta z)=h(z),\mbox{ say}.$

Note that in computing $w(z)$ , we need only compute $w(z)$ up to proportionality in $z$ as $\delta$ remains fixed in the update of $z_{i}$ . That is, it suffices to use $h(z)$ when updating a censored observation.

Independence sampler within Gibbs

1.

Choose initial values for $\delta$ and $\mathbf{z}=(z_{15},z_{16},\ldots,z_{20})$ .
2.
For $k=1,2,\ldots,n$ ;
- •
  
  Update $\delta|\mathbf{y},\mathbf{z}\sim{\rm Gamma}(40+\alpha,\sum_{i=1}^{20}x_{i}+\beta)$ ;
- •
  
  For $i=15,16,\ldots,20$ , draw $u_{i}\sim U(0,1)$ and set $z_{i}^{\prime}=(8/u_{i})^{1/3}$ .
  Accept the proposed value $z_{i}^{\prime}$ with probability $\min\{1,h(z_{i}^{\prime})/h(z_{i})\}$ , where $h(z)=z^{5}\exp(-\delta z)$ .
- •
  
  Store $\delta$ .

The algorithm works well, see the traceplot of $\delta$ below. The output is based upon a sample of 1000 iterations following a burn-in of 100 iterations. The estimated posterior mean and standard deviation of $\delta$ are $1.165$ and $0.2095$ , respectively.

Unnumbered Figure: Link

4.5 Diagnostics and rules of thumb

4.5.1 Motivation

1.

How do I decide when my chain has converged?
2.

When my chain has converged how do I decide if it is mixing quickly enough (or equivalently, if the chain is long enough to get accurate answers)?
3.

How do I choose my algorithm so that I do not need ”too many” iterations? (e.g. Would a $N(0,1)$ proposal have been better for the Independence Sampler? Would a $N(0,1/2)$ jump proposal have been better for the RWM?)
4.

What do we mean by ”better” in the previous question?!

Below we begin to answer the above questions. To assist us with answering the questions we consider four scenarios for obtaining samples from $X\sim N(0,1)$ .

•

Random walk Metropolis with $x^{(0)}=2$ and proposal distribution $Y|X=x\sim N(x,1)$ .
•

Independence Sampler with $x^{(0)}=2$ and proposal distribution $Y|X=x\sim N(0,2^{2})$ .
•

Random walk Metropolis with $x^{(0)}=100$ and proposal distribution $Y|X=x\sim N(x,1)$ .
•

Random walk Metropolis with $x^{(0)}=2$ and proposal distribution $Y|X=x\sim N(x,(0.1)^{2})$ .

Unnumbered Figure: Link

4.5.2 The Gelman Rubin Convergence Diagnostic

Let $\sigma^{2}=\mbox{Var}\left[{X_{i}}\right]$ then, if the $X_{i}$ are independent,

\mbox{Var}\left[{\overline{X}}\right]=\mbox{Var}\left[{\frac{1}{n}\sum_{i=1}^{% n}X_{i}}\right]=\frac{1}{n^{2}}\sum_{i=1}^{n}\mbox{Var}\left[{X_{i}}\right]=% \frac{\sigma^{2}}{n}.

Recap ANOVA for m groups of size n - comparison of two variance estimates; between groups (from the group means) and within groups (about the group means).

Suppose that you have $\{X_{ij}\}$ $(i=1,\ldots,m;j=1,\ldots,n)$ .

Let $\bar{X}_{i}=\frac{1}{n}\sum_{j=1}^{n}X_{ij}$ with $\mbox{Var}\left[{\bar{X}_{i}}\right]=\sigma^{2}/n$ and $\bar{X}=\frac{1}{m}\sum_{i=1}^{m}\bar{X}_{i}$ .

Note that

$\displaystyle\mathbb{E}\left[{\frac{1}{m-1}\sum_{i=1}^{m}(\bar{X}_{i}-\bar{X})% ^{2}}\right]$	$\displaystyle=$	$\displaystyle\frac{1}{m-1}\sum_{i=1}^{m}\mathbb{E}\left[{\bar{X}_{i}^{2}-2\bar% {X}_{i}\bar{X}+\bar{X}^{2}}\right]$
	$\displaystyle=$	$\displaystyle\frac{1}{m-1}\sum_{i=1}^{m}\left(\mathbb{E}\left[{\bar{X}_{i}^{2}% }\right]-\mathbb{E}\left[{\bar{X}^{2}}\right]\right)$
	$\displaystyle=$	$\displaystyle\frac{m}{m-1}\left\{\mathbb{E}\left[{\bar{X}_{1}^{2}}\right]-% \left(\frac{1}{m}\mathbb{E}\left[{\bar{X}_{1}^{2}}\right]+\frac{m-1}{m}\mathbb% {E}\left[{\bar{X}_{1}}\right]\mathbb{E}\left[{\bar{X}_{2}}\right]\right)\right\}$
	$\displaystyle=$	$\displaystyle\frac{m}{m-1}\frac{m-1}{m}\left(\mathbb{E}\left[{\bar{X}_{1}^{2}}% \right]-\mathbb{E}\left[{\bar{X}_{1}}\right]^{2}\right)$
	$\displaystyle=$	$\displaystyle\mbox{Var}\left[{\bar{X}_{1}}\right]=\frac{\sigma^{2}}{n}$

Hence $B=\frac{n}{m-1}\sum_{i=1}^{m}(\bar{x}_{i}-\bar{x})^{2}$ is an unbiased estimate of $\sigma^{2}$ . (This estimate is based upon between group variability.)

Also

$\displaystyle\mathbb{E}\left[{\frac{1}{m}\sum_{i=1}^{m}\frac{1}{n-1}\sum_{j=1}% ^{n}(X_{ij}-\bar{X}_{i})^{2}}\right]$	$\displaystyle=$	$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}\left[{\frac{1}{n-1}\sum_{j=1}% ^{n}(X_{ij}-\bar{X}_{i})^{2}}\right]$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left[{\frac{1}{n-1}\sum_{j=1}^{n}(X_{1j}-\bar{X}_{1})^% {2}}\right]$
	$\displaystyle=$	$\displaystyle\sigma^{2}$

Hence $W=\frac{1}{m}\sum_{i=1}^{m}\frac{1}{n-1}\sum_{j=1}^{n}(x_{ij}-\bar{x}_{i})^{2}$ is also unbiased estimate of $\sigma^{2}$ . (This estimate is based upon within group variability.)

Finally, $V=\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(x_{ij}-\bar{x})^{2}$ is also an unbiased estimate of $\sigma^{2}$ . (This estimate is based upon total variability.)
Let

$\displaystyle V$	$\displaystyle=$	$\displaystyle\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(x_{ij}-\bar{x})^{2}$
	$\displaystyle=$	$\displaystyle\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(x_{ij}-\bar{x}_{i}+% \bar{x}_{i}-\bar{x})^{2}$
	$\displaystyle=$	$\displaystyle\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(x_{ij}-\bar{x}_{i})^{2% }+\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(\bar{x}_{i}-\bar{x})^{2}$
	$\displaystyle=$	$\displaystyle\frac{m(n-1)}{mn-1}W+\frac{m-1}{(mn-1)}B$
	$\displaystyle\approx$	$\displaystyle\left(1-\frac{1}{n}\right)W+\frac{1}{n}B.$

Thus $B$ , $W$ and $V$ are estimates of the variance based upon the between group, within group and total variation, respectively.

We will run $m$ (e.g. 5) separate chains of length $n$ (e.g. 1000) from different starting points which are more spread out than (”overdispersed with respect to”) the posterior and compare two estimates of the variance of the stationary distribution.

The GR Convergence Diagnostic Procedure

1.

Run $m$ chains for $2n$ iterations and discard the first $n$ iterations of each chain as burn-in. (Keep the last $n$ iterations from each chain.)
Let $\bar{\theta}_{i}$ denote the mean parameter value of chain $i$ and let $\bar{\theta}_{\cdot}$ denote the overall mean. Let $s_{i}^{2}$ denote the variance of chain $i$ . Then

$B=\frac{n\sum_{i=1}^{m}(\bar{\theta}_{i}-\bar{\theta}_{\cdot})^{2}}{m-1}$

and

$W=\frac{1}{m}\sum_{i=1}^{m}s_{i}^{2}.$
2.

Calculate the between chain variance estimate, $B$ , and the within chain variance estimate, $W$ .
3.

Calculate $V=(1-1/n)W+(1/n)B$ (nearly the total variance of all $m$ chains together).
4.

If the chains have all converged to the same distribution then $W,~{}B,$ and $V$ should all be estimates of the variance of the stationary distribution, $\sigma^{2}$ . i.e. They should all be similar.
5.

Evaluate $\hat{R}:=\sqrt{V/W}$ .
6.

If the chains have not yet converged then $V$ should be $>\sigma^{2}$ (because the chains are still too spread out) and $W$ should be $<\sigma^{2}$ (since each chain has not yet traversed the state-space).
7.

So if the chains have not converged then $\hat{R}>1$ otherwise $\hat{R}\approx 1$ .
8.

Gelman (1996) suggests accepting convergence if $\hat{R}<1.2$ .

The figure shows a plot of the GR statistics using 2 chains: the fast mixing RWM started at $x^{\left({0}\right)}=2$ , and the independence sampler started at $x^{\left({0}\right)}=2$ .

Unnumbered Figure: Link

The figure shows a plot of the GR statistics using 3 chains: the fast mixing RWM started at $x^{\left({0}\right)}=2$ , the independence sampler started at $x^{(0)}=2$ , and the RWM started at $x^{(0)}=100$ .

What does the Gelman-Rubin statistics tell us?

It is very helpful for checking whether or not the distribution is multi-model. In that, if we start at a wide range of different values we will converge to different modes if they exist. Then we need to ensure the MCMC can efficiently move between modes. It also gives an indication of a safe burn-in period. For example, if for $n=1000$ $\hat{R}\approx 1$ , this tell us that a burn-in of 1000 iterations should be sufficient to ensure convergence.

4.6 Independent versus dependent samples

If you had $1000$ independent samples $x_{1},\dots,x_{1000}$ from a random variable $X$ with distribution $f$ , then you could estimate

\hat{\mu}=\frac{1}{1000}\sum_{i=1}^{1000}x_{i}\approx\mu=\mathbb{E}\left[{X}% \right],

and

\hat{\sigma^{2}}=\frac{1}{999}\sum_{i=1}^{1000}(x^{2}_{i}-\hat{\mu})^{2}% \approx\sigma^{2}=\mbox{Var}\left[{X}\right].

But you could also estimate the error in each of these approximations

\mbox{Var}\left[{\hat{\mu}}\right]=\frac{\sigma^{2}}{n}\approx\frac{\widehat{% \sigma}^{2}}{n},~{}~{}~{}~{}\mbox{Var}\left[{\hat{\sigma^{2}}}\right]=\frac{2% \sigma^{4}}{n-1}\approx\frac{2\hat{\sigma}^{4}}{n-1}.

Suppose that we have a chain of length $1000$ , $x^{\left({1}\right)},\dots,x^{\left({1000}\right)}$ . Each element of the chain depends on the previous element. Indeed, because the previous element depends on the one before that, there is dependency between $x^{\left({i}\right)}$ and $x^{\left({i-2}\right)}$ etc. Because of this (positive) correlation there is generally less information in the 1000 elements of the chain than there would have been in 1000 independent draws.

If the chain is stationary then $\rho_{1}:=\mbox{Corr}\left[{X^{\left(i\right)},X^{\left({i-1}\right)}}\right]$ should not depend on $i$ . $\rho_{1}$ is called the lag-1 autocorrelation. Similarly define $\rho_{k}:=\mbox{Corr}\left[{X^{\left(i\right)},X^{\left({i-k}\right)}}\right]$ as the lag-k autocorrelation. Note that the lag-0 autocorrelation is 1, in that, any random variable is perfectly correlated with itself!

The plots below show the autocorrelations up to lag-50 for four chains sampling from $N(0,1)$ : RWM with N(0,1), Indep. Sampler with N(0,4), RWM with N(0,1) (start at 100), RWM with N(0,0.01). These correspond to the traceplots given at the start of the Section.

Unnumbered Figure: Link

The RWM with proposal $N(0,1)$ performs quite well with $\rho_{k}\approx 0$ for $k>20$ . $\rho_{1}=0.8$ and $\rho_{k}\approx\rho_{1}^{k}$ . The Independence Sampler works very well with $\rho_{1}=0.25$ and all other $\rho_{k}\approx 0$ . This means that there is a small cost for this algorithm over an independent sample. For $x^{\left({0}\right)}=100$ , the effect of the starting value has a significant effect on the acf plot. The RWM with proposal $N(0,0.01)$ mixes slowly and consequently $\rho_{k}$ decreases only slowly as $k$ increases.

4.6.1 Integrated Autocorrelation Time and Effective Sample Size

The integrated autocorrelation time for an infinite Markov chain is defined as

I_{ACT}:=1+2\sum_{i=1}^{\infty}\rho_{i}.

In practice we only have a finite chain and so calculate

\hat{I}_{ACT}:=1+2\sum_{i=1}^{K-1}\rho_{i},

where $K$ is the first time that $\rho_{K}$ is below some threshold (e.g. $0.05$ ). For the examples above $\hat{I}_{ACT}$ is $11.3,~{}1.6,~{}159.2,~{}278.2$ . The larger the integrated ACT the higher the correlation and so the less information there is in the chain. In practice we only estimate the $I_{ACT}$ after burn-in (since we only use the part of the chain after burn-in). The $I_{ACT}$ estimate for the final 800 iterations of chain started at $x^{\left({0}\right)}=100$ is $14.6$ .

For $n$ independent samples

\mbox{Var}\left[{\hat{\mu}}\right]=\frac{\sigma^{2}}{n}\approx\frac{\hat{% \sigma}^{2}}{n}.

Define the effective sample size $n_{eff}:=n/I_{ACT}$ .

Lemma For a stationary series of correlated variables with $\sum_{i=1}^{\infty}i\left|{\rho_{i}}\right|<\infty$ ,

\mbox{Var}\left[{\hat{\mu}}\right]\approx\frac{\hat{\sigma}^{2}}{n_{eff}}.

Proof:

$\displaystyle\mbox{Var}\left[{\hat{\mu}}\right]$	$\displaystyle=$	$\displaystyle\mbox{Var}\left[{\frac{1}{n}\sum_{i=1}^{n}X^{(i)}}\right]$
	$\displaystyle=$	$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\mbox{Cov}\left[{X^{(i% )},X^{(j)}}\right]$
	$\displaystyle=$	$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\mbox{Var}\left[{X^{(i)}}\right]+% \frac{2}{n^{2}}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}\mbox{Cov}\left[{X^{(i)},X^{(j)% }}\right]$
	$\displaystyle=$	$\displaystyle\frac{1}{n^{2}}\left\{n\mbox{Var}\left[{X^{(1)}}\right]+2\sum_{i=% 1}^{n-1}(n-i)\mbox{Cov}\left[{X^{(0)},X^{(i)}}\right]\right\}$
	$\displaystyle\approx$	$\displaystyle\frac{1}{n}\left\{\mbox{Var}\left[{X^{(1)}}\right]+2\sum_{i=1}^{n% -1}\mbox{Cov}\left[{X^{(0)},X^{(i)}}\right]\right\}$
	$\displaystyle=$	$\displaystyle\frac{1}{n}\sigma^{2}\left\{1+2\sum_{i=1}^{n-1}\rho_{i}\right\}$
	$\displaystyle\approx$	$\displaystyle\frac{1}{n}\hat{\sigma}^{2}\left\{1+2\sum_{i=1}^{\infty}\rho_{i}\right\}$
	$\displaystyle=$	$\displaystyle\frac{1}{n}\hat{\sigma}^{2}I_{ACT}=\frac{\hat{\sigma}^{2}}{n_{eff% }}.$

The effective sample sizes for our four chains are therefore: $1000/11.3=88.5$ , $1000/1.6=625$ , $800/14.6=54.8$ (excluding burn-in), and $1000/278.2=3.6$ . The $1000$ values in the final chain contain less information than would a sample of $4$ independent values!

4.6.2 Answers

1.

How do I decide when my chain has converged? Both by eye and (if feasible) using the Gelman-Rubin statistic.
2.

When my chain has converged how do I decide if it is mixing quickly enough (or equivalently, if the chain is long enough to get accurate answers)? By eye and look at the effective sample size and the standard error of your estimate.
3.

How do I choose my algorithm so that I do not need ”too many” iterations? (e.g. Would a $N(0,1)$ proposal have been better for the Independence Sampler? Would a $N(0,1/2)$ jump proposal have been better for the RWM?) To be discussed!
4.

What do we mean by ”better” in the previous question?! Larger effective sample size for the same length of chain.

4.7 R Practice Questions: Independence Sampler

1.

Write R code to implement the independence sampler for the genetics example. (Section 4.3.2).
2.

Write R code to implement the independence sampler within Gibbs for the censored data. (Section 4.4).
3.
Suppose that, there exists $\beta>0$ such that, $\mathbf{X}$ has pdf,

$\pi_{\beta}(\mathbf{x})\propto\left\{\begin{array}[]{ll}\exp(-\beta(x_{1}-x_{2% })^{2})&\mbox{if }0\leq x_{1},x_{2}\leq 10\\ 0&\mbox{otherwise}.\end{array}\right.$

For fixed values of $\beta$ . Write and implement the following independence sampler programs for obtaining samples from $\mathbf{X}$ .
1. (a)
  
  $q(\mathbf{x})=\frac{1}{100}$ for $0\leq x_{1},x_{2}\leq 10$ , a uniform proposal on $[0,10]$ .
2. (b)
  
  $x_{1}\sim U(0,10)$ and $x_{2}|x_{1}\sim N(x_{1},\sigma^{2})$ giving
  
  $q(\mathbf{x})=\frac{1}{10}\times\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}% {2\sigma^{2}}(x_{2}-x_{1})^{2}\right),\hskip 14.226378pt0\leq x_{1}\leq 10,x_{% 2}\in\mathbb{R}.$
How do the two independent samplers compare with $\beta=2$ ?

Theoretical Question: For $\beta=2$ , comment upon the choice of $\sigma$ .

4.8 Lab Session: Convergence diagnostics

The files background.R, RWM1d.R and Indep1d.R have the code for applying the random walk Metropolis algorithm and the independence sampler to five 1-d distributions. The files background.R gives the background functions whilst RWM1d.R and Indep1d.R run the random walk and independence sampler algorithms.

One observation with the R code is that it uses the $\log$ acceptance probability.

For the 1D functions you have a choice of target (stationary) distributions

•

Gaussian: $\pi(x)\propto e^{-\frac{1}{2}x^{2}}$
•

Laplace: $\pi(x)\propto e^{-\left|{x}\right|}$
•

t5: $\pi(x)\propto 1/(1+x^{2}/5)^{3}$
•

Cauchy: $\pi(x)\propto 1/(1+x^{2})$
•

Bimodal: $\pi(x)\propto 0.2e^{-\frac{1}{2}(x+8)^{2}}+0.8\times\frac{1}{2}e^{-\frac{1}{8}% (x-6)^{2}}$

Of course in real life it is easy to obtain an independent sample from each of these distributions and we would never run an MCMC chain for this purpose. However these distributions have different properties and it can be instructive to see how the behaviour of the chain varies in these simple examples. These rules of thumb can then be carried over to more complex scenarios.

You may choose your proposal distribution from Gaussian, Laplace, t5, or Cauchy; you also choose the scale parameter $\lambda$ for the proposal. For example if you chose a Gaussian proposal with scale parameter $\lambda=3$ when calling RWM1d then the proposed jump is $N(0,9)$ , i.e. $x_{prop}\sim N(x_{curr},9)$ . The same proposal for the independence sampler would mean $x_{prop}\sim N(0,9)$ .

Try running RWM1d or Indep1d using the example function call at the bottom of the file - you will see that it produces some basic output.

Important for intuition Note that the Gaussian has lighter tails than the Laplace which has lighter tails than the $t_{5}$ , which has lighter tails than the Cauchy.

Introducing coda Spend no more than 15 minutes on this section!

Run the 1D RWM function for $1000$ iterations, starting at $x_{0}=0$ , using a Gaussian proposal, a Gaussian target, and a scale parameter of $0.5$

a1=RWM1d(1000,0.5,"Gaussian","Gaussian",0)

a1 is a vector of length $1000$ (check!). To use the coda package we first load the library

library(coda)

We then need to turn a1 into an mcmc object.

a1.mc=mcmc(a1)

Try out the following, work out and note down what they do; if necessary use the help facility.
summary(a1.mc)
rejectionRate(a1.mc) $~{}~{}$ How does this relate to acceptance rate?
plot(a1.mc,smooth=F)
traceplot(a1.mc,smooth=F)
densityplot(a1.mc,smooth=F)
effectiveSize(a1.mc)
autocorr.diag(a1.mc)
autocorr.plot(a1.mc,lag.max=50)
cumuplot(a1.mc) (Why might this be useful?)

To use the Gelman-Rubin convergence test we need to have a few more chains started at very different points.
a2=RWM1d(1000,0.5,"Gaussian","Gaussian",4)
a3=RWM1d(1000,0.5,"Gaussian","Gaussian",-4)

Turn these into mcmc objects a2.mc and a3.mc.

Create an mcmc list
a.list $\leftarrow$ mcmc.list(a1=a1.mc,a2=a2.mc,a3=a3.mc)

Now try the following
gelman.diag(a.list)
gelman.plot(a.list)

You should find that the chains take about 200 iterations to converge (check using the trace plots). You would therefore only use the last 800 iterations to calculate any quantities of interest (e.g. mean and variance).

coda has a number of other functions - feel free to investigate in your own time.
https://cran.r-project.org/web/packages/coda/coda.pdf

The main idea of this part of the lab is that you should investigate how the choice of scaling and proposal distribution affects behaviour for different targets. The instructions below comprise a minimal set of tasks - there are many further investigations if you have time!

Throughout this part of the lab it will be better for you to complete one investigation (e.g. Qn1 below) than to partially complete many.

Mixing

First try RWM1d(); each time you should run the algorithm for 10000 iterations and start at $x_{0}=0$ which is in the main mass of the target distribution (so we do not have to worry about convergence). You should try a variety of different scalings of the proposal from $\lambda=0.1$ to $\lambda=10$ . For each chain, look at the trace plot to see how the chain behaves; note down the acceptance rate, and the effective sample size.

Repeat some of the simulations at least 3 times to check for consistency.

For each of the four stages below (and any other possibilities which you wish to try) answer the following:

•

Which scaling gives the best effective sample size?
•

What is that best effective sample size?
•

What is the acceptance rate at this best scaling?
•

How did the chains behave?

Be sure to examine the traceplots as well as calculating the effective sample size.

1.

Use a Gaussian target and Gaussian proposal.
2.

Repeat (1) for a Gaussian target but instead using a Cauchy proposal.
3.

Repeat (1) and (2) using a Cauchy target with a Gaussian proposal (also try a Cauchy proposal if time allows).
4.

Repeat (1) and (2) but using the bimodal target with a Gaussian proposal (also try a Cauchy proposal if time allows).

What do you conclude? Can you explain your findings?

Now look at Indep1d() following exactly the same instructions as above; your findings and conclusions should be different!

Convergence

Note If time is short then the minimal set of investigations for this section comprises of Qn1 and Qn3 below for RWM1d().

For all chains in this section you should start with $x_{0}=50$ . In each case run a chain with $\lambda=\hat{\lambda}$ , the best scaling that you found in the previous section - and with $\lambda=2\hat{\lambda}$ . In each case investigate the number of iterations which it takes for the chain to converge, either using the traceplot alone or by running several chains, looking at the traceplots and using the GR plot.

Investigate RWM1d using the following combinations

1.

Gaussian target and Gaussian proposal.
2.

Gaussian target and Cauchy proposal.
3.

Cauchy target and Gaussian proposal.
4.

Cauchy target and Cauchy proposal.

Feel free to try other combinations if there is time!

Now investigate Indep1d using the same combinations of target and proposal as for RWM1d().

4.9 Practice Questions

1.
Suppose that $f(x)\propto\sqrt{1+x^{2}}\exp(-x^{2}/2)$ $(x\in\mathbb{R})$ . The following 3 independence samplers are suggested for obtaining a sample from $f(x)$ .
- Sampler 1.
  
  Proposal distribution $N(0,1)$ with proposal probability density function $q(x)=\exp(-x^{2}/2)/\sqrt{2\pi}$ .
- Sampler 2.
  
  Proposal distribution $Cauchy(50,1)$ with proposal probability density function $q(x)=\{1/(\pi(1+(x-50)^{2}))$ .
- Sampler 3.
  
  Proposal distribution $t_{5}$ with proposal probability density function $q(x)=\frac{8}{\pi\sqrt{45}}(1+x^{2}/5)^{-3}$ .
1. (a)
  
  Why is Sampler 3 the best choice for obtaining samples from $f(x)$ , and why are the other two independence samplers unsuitable?
2. (b)
  
  A sample of size 1000 was obtained using Sampler 3 from $f(x)$ and the output is summarised in the form of a traceplot and acf plot in Figure 1 with an acceptance rate of 0.883.
  
  Unnumbered Figure: Link
  
  Traceplot and acf plot of sample of size 1000 from $f(x)$ using Sampler 3. Comment upon the performance of Sampler 3 from $f(x)$ .
2.

A doctor is interested in the survival time of patients (in years) following major surgery. It is proposed that the survival time, $S$ , follows a Gamma distribution with shape parameter 2 and scale parameter $\delta$ . That is, $S$ has probability density function,

$f_{S}(s)=\delta^{2}s\exp(-\delta s)\hskip 14.226378pt(s\geq 0).$
Suppose that $N$ patients were monitored and their survival times $s_{1},s_{2},\ldots,s_{N}$ are recorded.
1. (a)
  
  Given that the prior distribution on $\delta$ is ${\rm Gamma}(A,B)$ , derive the posterior distribution $\delta$ given ${\mathbf{s}}=(s_{1},s_{2},\ldots,s_{N})$ .
Due to limited resources patients are only monitored for 5 years following surgery. Therefore if the survival time of the patient is less than or equal to 5 years it is known exactly. However, if the survival time of the patient is more than 5 years all that is known is that the survival time exceeds 5 years (censored data). Censored survival times are denoted by $\ast$ .
Suppose that for $1\leq M<N$ , $s_{1},s_{2},\ldots,s_{M}$ are observed exactly (survival time less than or equal to 5 years) and for the remaining $N-M$ patients their survival times $s_{M+1},s_{M+2},\ldots,s_{N}$ are censored at 5 years.
1. (b)
  
  The survival time, $S$ , given that it exceeds 5 years has probability density function
  
  $f_{S>5}(s|\delta)=\frac{\delta^{2}s\exp(-\delta s)}{(5\delta+1)\exp(-5\delta)}% \hskip 14.226378pt(s\geq 5).$
  
  Describe an independence sampler for simulating from $S^{\ast}$ .
2. (c)
  
  Describe a Gibbs sampler for obtaining samples from the posterior distribution of $\delta$ given the observed data $s_{1},s_{2},\ldots,s_{M},5^{\ast},5^{\ast},\ldots,5^{\ast}$ .
3.

Suppose that $X$ is a non-negative random variable. That is, for $x<0$ , $f(x)=0$ . A reflected random walk proposal is used. That is, if the MCMC is currently at $x$ , the proposed value $y$ , is generated from

$Y\sim|N(x,\sigma^{2})|.$

Let $\phi(z)$ denote the probability density function of a standard Normal distribution.
1. (a)
  
  Show that $q(x,y)=\phi((x-y)/\sigma)+\phi((x+y)/\sigma)$ .
2. (b)
  
  Hence, or otherwise, show that the acceptance probability, $\alpha(x,y)$ , for the proposed move from $x$ to $y$ satisfies
  
  $\alpha(x,y)=\min\left\{1,\frac{f(y)}{f(x)}\right\}.$

	$\displaystyle f(x)$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left\{\frac{1}{\sqrt{2\pi}}\exp(-(x+1)^{2}/2)+\frac{1% }{2\sqrt{2\pi}}\exp(-x^{2}/8)\right\}\hskip 14.226378pt\mbox{(black)}$
	$\displaystyle q(y)$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{2\pi}\sqrt{5}}\exp(-(y+0.5)^{2}/10)\hskip 14.22637% 8pt\mbox{(red)}$

	$\displaystyle X$	$\displaystyle\sim$	$\displaystyle\left\{\begin{array}[]{ll}N(-1,1)&w.p.~{}0.5\\ N(0,4)&w.p.~{}0.5\end{array}\right.$
	$\displaystyle Y$	$\displaystyle\sim$	$\displaystyle N(-0.5,5)$

4 Week 9: The Metropolis-Hastings Algorithm

4.1 Introduction to the Metropolis-Hastings algorithm

4.2 The Metropolis-Hastings algorithm in a Bayesian context

4.2.1 Componentwise MCMC

4.3 Independence sampler

4.3.1 Introduction

4.3.2 Genetics example (revisited)

4.3.3 Why should q be over-dispersed relative to f?

4.4 Independence sampler within Gibbs: Censored data

4.5 Diagnostics and rules of thumb

4.5.1 Motivation

4.5.2 The Gelman Rubin Convergence Diagnostic

4.6 Independent versus dependent samples

4.6.1 Integrated Autocorrelation Time and Effective Sample Size

4.6.2 Answers

4.7 R Practice Questions: Independence Sampler

4.8 Lab Session: Convergence diagnostics

4.9 Practice Questions

4.3.3 Why should $q$ be over-dispersed relative to $f$ ?