\chapter{Non-conjugate multi-parameter problems} In this chapter we will study some multi-parameter problems in which the prior distribution does not have to be conjugate. Inferences are made by using techniques which simulate realisations from the posterior distribution. These methods are generally referred to as \textit{Markov Chain Monte Carlo} techniques, and often abbreviated to MCMC. There are many different MCMC techniques, but we only have time to look briefly at two of the most fundamental. The first is the \textit{Gibbs sampler}, which was at the forefront of the recent MCMC revolution, and the second is generally known as \textit{Metropolis-Hastings} sampling. In fact, MCMC schemes based on the combination of these two fundamental techniques are still at the forefront of MCMC research. \section{Why is inference not straightforward in non-conjugate problems?} \label{sec:infproblem} \begin{example}~\\ Consider again the problem in section~\ref{sec:nochtheory} in which \label{ex:semiconjnormal} we have a random sample from a normal distribution where both the mean~$\mu$ and the precision~$\tau$ are unknown, that is, $X_i|\mu,\tau\sim N(\mu,1/\tau)$, $i=1,2,\ldots,n$ (independent). In this section, we showed that a $NGa$ prior for $(\mu,\tau)^T$ was conjugate, that is, if we used a $NGa(b,c,g,h)$ prior distribution for $(\mu,\tau)^T$ then the posterior was a $NGa(B,C,G,H)$ distribution. But what if a $NGa(b,c,g,h)$ prior distribution does not adequately represent our prior beliefs? Suppose instead that our prior beliefs are represented by independent priors for the parameters, with \[ \mu\sim N\left(b,\frac{1}{c}\right) \quad\quad\text{and}\quad\quad \tau\sim Ga(g,h) \] for known values $b$, $c$, $g$ and $h$. What is the posterior distribution for $(\mu,\tau)^T$? \end{example} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \textcolor{blue}{ %\begin{soln} Previously we have seen that the likelihood function is \begin{align*} f(\xvec|\mu,\tau)&=\left(\frac{\tau}{2\pi}\right)^{n/2} \exp\left[-\frac{n\tau}{2} \left\{s^2+(\bar x-\mu)^2\right\}\right]. %\end{align*} %\begin{align} \intertext{Using Bayes Theorem, the posterior density is} \pi(\mu,\tau|\xvec)&\propto\pi(\mu,\tau)\,f(\xvec|\mu,\tau) \\ &\propto\pi(\mu)\,\pi(\tau)\,f(\xvec|\mu,\tau) \intertext{and so, for $\mu\in\R$, $\tau>0$} \pi(\mu,\tau|\xvec)&\propto \left(\frac{c}{2\pi}\right)^{1/2}\exp\left\{-\frac{c}{2}(\mu-b)^2\right\} \times \frac{h^g\tau^{g-1}e^{-h\tau}}{\Gamma(g)}\\ &\hspace*{3cm}\times \tau^{\frac{n}{2}} \exp\left[-\frac{n\tau}{2} \left\{s^2+(\bar x-\mu)^2\right\}\right] \\ &\propto\tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}. \end{align*} This is not a $NGa(B,C,G,H)$ density as it cannot be written in the form \[ k\,\tau^{G-\frac{1}{2}}\exp\left\{-\frac{\tau}{2} \left[C(\mu-B)^2+2H\right]\right\} \] for any choice of $B$, $C$, $G$ and $H$. In fact $\pi(\mu,\tau|\xvec)$ is not the density of any standard distribution. Can we plot this density? Before we do, we need to know the constant of proportionality, $k$. As the posterior density must integrate to one \begin{align*} \int_{-\infty}^\infty&\int_0^\infty k\tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau\,d\mu=1\\ k^{-1}&=\int_{-\infty}^\infty\int_0^\infty \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau\,d\mu \end{align*} Therefore the posterior density is, for $\mu\in\R$, $\tau>0$ \[ \pi(\mu,\tau|\xvec)=\bigfrac{ \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}} {\int_{-\infty}^\infty\int_0^\infty \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau\,d\mu}. \] %\end{soln} %%Hack for vspace \vspace{5cm} }}} What is the posterior mean of $\mu$ and of $\tau$? What are their marginal distributions? How can we calculate the moments $E(\mu^{m_1}\tau^{m_2}|\xvec)$ of this posterior distribution? Now \begin{align*} \pi(\mu|\xvec) &=\int_0^\infty\pi(\mu,\tau|\xvec)\,d\tau\\ &=\bigfrac{\int_0^\infty \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau} {\int_{-\infty}^\infty\int_0^\infty \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau\,d\mu} \intertext{and} \pi(\tau|\xvec) &=\int_{-\infty}^\infty\pi(\mu,\tau|\xvec)\,d\mu\\ &=\bigfrac{\int_{-\infty}^\infty \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\mu} {\int_{-\infty}^\infty\int_0^\infty \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau\,d\mu}. \end{align*} In general, the moments are \begin{align*} E(\mu^{m_1}\tau^{m_2}|\xvec) &=\int_{-\infty}^\infty\int_0^\infty\mu^{m_1}\tau^{m_2}\pi(\mu,\tau|\xvec)\,d\tau\,d\mu\\ &=\bigfrac{\int_{-\infty}^\infty\int_0^\infty \mu^{m_1}\tau^{m_2}\times\tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau\,d\mu} {\int_{-\infty}^\infty\int_0^\infty \tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\, d\tau\,d\mu}. \end{align*} These integrals cannot be determined analytically, though it is possible to use numerical integration methods or approximations (for large~$n$). However, in general, the accuracy of the numerical approximation to the integral deteriorates as the dimension of the integral increases. \subsubsection*{Comment} The above shows how not using a conjugate prior distribution can cause many basic problems such as plotting the posterior density or determining posterior moments. But having to use conjugate priors is far too restrictive for many real data analyses: (i)~our prior beliefs may not be captured using a conjugate prior; (ii) most models for complex data do not have conjugate priors. It was for these reasons that until relatively recently (say mid-1990s), practical Bayesian inference for real complex problems was either not feasible or only undertaken by the dedicated few prepared to develop bespoke computer code to numerically evaluate all the integrals etc. \newpage \section{Simulation-based inference} One way to get around the problem of having to work out integrals (like those in the previous section) is to base inferences on simulated realisations from the posterior distribution. This is the fundamental idea behind MCMC methods. If we could simulate from the posterior distribution then we could use a very large sample of realisations to determine posterior means, standard deviations, correlations, joint densities, marginal densities etc. As an example, imagine you wanted to know about the standard normal distribution -- its shape, its mean, its standard deviation -- but didn't know any mathematics so that you couldn't derive say the distribution's zero mean and unit variance. However you've been given a ``black box'' which can simulate realisations from this distribution. Here we'll use the \texttt{R} function \texttt{rnorm()} as the black box simulator. If you decide to generate 1K realisations the output might look something like the top row of Figure~\ref{fig:nrs1}. \begin{figure}[!b] \begin{center} \includegraphics[width=0.9\textwidth]{nrs1} \caption{Summaries of the 1K realisations from the black box simulator} \label{fig:nrs1} \end{center} \end{figure} The top left plot shows the trace plot of the output, that is, the realisations from the black box sampler in the order they are produced. The next plot along the top row shows the autocorrelation (ACF) plot. This shows how correlated the realisations are at different lags. We know that the simulator \texttt{rnorm()} produces independent realisations and so the (sample) correlation between say consecutive values $corr(x_i,x_{i+1})$ will be almost zero. This is also the case for correlations at all positive lags. Finally the lag~0 autocorrelation $corr(x_i,x_i)$ must be one (by definition). The sample ACF plot is consistent with all of these ``claims''. Finally the top right plot is a density histogram of the realisations. This too is consistent with the standard normal density (which is also shown). We can also estimate various quantities of the standard normal distribution; for example: \begin{verbatim} 1st Qu. Median Mean 3rd Qu. St.Dev. -0.65240 -0.00130 -0.00192 0.64810 0.96049 \end{verbatim} Here we see that the mean and median are around zero and the standard deviation is around one. The second row of plots in the figure is another collection of 1K realisations from the black box simulator. These look very similar to those on the top row but are slightly different due to the stochasticity (random nature) of the simulator. This output has the following numerical summaries: \begin{verbatim} 1st Qu. Median Mean 3rd Qu. St.Dev. -0.69880 -0.09637 -0.03274 0.67330 0.99599 \end{verbatim} Again these numerical summaries are slightly different but essentially tell the same story. In fact we know from previous modules that there is sample variability in estimates of means from random samples. So if we use the simulator again (twice) to obtain 10K realisations, we will get even more similar looking output; see Figure~\ref{fig:nrs2}. \begin{figure}[t!] \begin{center} \includegraphics[width=0.9\textwidth]{nrs2} \caption{Summaries of the 10K realisations from the black box simulator} \label{fig:nrs2} \end{center} \end{figure} The numerical summaries from these outputs are \newpage \begin{verbatim} 1st Qu. Median Mean 3rd Qu. St.Dev. -0.66820 -0.00048 0.00370 0.68070 0.99593 -0.67130 0.01354 0.01008 0.67920 1.00691 \end{verbatim} Here we have much less sampling variability in our estimates due to the larger sample size. In fact we can estimate any ``population'' quantity to any required accuracy simply by simulating a large enough collection of realisations. These analyses show how we can make inferences, calculate means, variances, densities etc by using realisations from a distribution. In the rest of this chapter, we will look into how we can construct algorithms for simulating from (complex) posterior distributions, from which we can then make inferences. \section{Motivation for MCMC methods} The example in section~\ref{sec:infproblem} showed that using non-conjugate priors can be problematic. MCMC methods address this problem by providing an algorithm which simulates realisations from the posterior distribution. We consider a generic case where we want to simulate realisations of two random variables $X$ and $Y$ with joint density $f(x,y)$. This joint density can be factorised as \[ f(x,y) = f(x)f(y|x) \] and so we can simulate from $f(x,y)$ by first simulating $X=x$ from $f(x)$, and then simulating $Y=y$ from $f(y|x)$. On the other hand, we can write \[ f(x,y) = f(y)f(x|y) \] and so simulate $Y=y$ from $f(y)$ and then $X=x$ from $f(x|y)$. We have already seen that dealing with conditional posterior distributions is straightforward when the prior is semi-conjugate, so let's assume that simulating from $f(y|x)$ and $f(x|y)$ is straightforward. The key problem with using either of the above methods is that, in general, we can't simulate from the marginal distribution, $f(x)$ and $f(y)$. For the moment, suppose we can simulate from the marginal distribution for~$X$, that is, we have an $X=x$ from $f(x)$. We can now simulate a $Y=y$ from $f(y|x)$ to give a pair $(x,y)$ from the bivariate density. Given that this pair is from the bivariate density, the $y$ value must be from the marginal $f(y)$, and so we can simulate an $X=x'$ from $f(x|y)$ to give a new pair $(x',y)$ also from the joint density. But now $x^\prime$ is from the marginal $f(x)$, and so we can simulate a $Y=y'$ from $f(y|X=x')$ to give a new pair $(x',y')$ also from the joint density. And we can keep going. This alternate sampling from conditional distributions defines a bivariate Markov chain, and the above is an intuitive explanation for why $f(x,y)$ is its stationary distribution. Thus being able to simulate easily from conditional distributions is key to this methodology. \section{The Gibbs sampler} Suppose we want to generate realisations from the posterior density $\pi(\thetavec|\xvec)$, where $\thetavece$, and that we can simulate from the full conditional distributions (FCDs) \[ \pi(\theta_i|\theta_1,\ldots,\theta_{i-1},\theta_{i+1},\ldots,\theta_p,\xvec) = \pi(\theta_i|\cdot), \qquad i=1,2,\ldots,p. \] The Gibbs sampler follows the following algorithm: \begin{enumerate} \item Initialise the iteration counter to $j=1$.\\ Initialise the state of the chain to $\thetavec^{(0)}=(\theta_1^{(0)},\ldots,\theta_p^{(0)})^T$. \item Obtain a new value $\thetavec^{(j)}$ from $\thetavec^{(j-1)}$ by successive generation of values \begin{align*} \theta_1^{(j)} &\sim \pi(\theta_1|\theta_2^{(j-1)},\theta_3^{(j-1)}, \theta_4^{(j-1)}, \theta_5^{(j-1)}, \ldots, \theta_p^{(j-1)},\xvec) \\ \theta_2^{(j)} &\sim \pi(\theta_2|\theta_1^{(j)},\theta_3^{(j-1)}, \theta_4^{(j-1)}, \theta_5^{(j-1)}, \ldots, \theta_p^{(j-1)},\xvec) \\ \theta_3^{(j)} &\sim \pi(\theta_3|\theta_1^{(j)},\theta_2^{(j)},\theta_4^{(j-1)}, \theta_5^{(j-1)}, \ldots, \theta_p^{(j-1)},\xvec) \\ \theta_4^{(j)} &\sim \pi(\theta_4|\theta_1^{(j)},\theta_2^{(j)}, \theta_3^{(j)}, \theta_5^{(j-1)}, \ldots, \theta_p^{(j-1)},\xvec) \\ \vdots \quad & \qquad \vdots \qquad\qquad \vdots \\ \theta_p^{(j)} &\sim \pi(\theta_p|\theta_1^{(j)},\theta_2^{(j)},\ldots,\theta_{p-1}^{(j)},\xvec) \end{align*} \item Change counter $j$ to $j+1$, and return to step 2. \end{enumerate} This algorithm defines a homogeneous Markov chain as each simulated value depends only on the previous simulated value and not on any other previous values or the iteration counter $j$. It can be shown that $\pi(\thetavec|\xvec)$ is the stationary distribution of this chain and so if we simulate realisations by using a Gibbs sampler, eventually the Markov chain will converge to the required posterior distribution. \subsection{Processing output from a Gibbs sampler} \subsubsection*{Burn-in period} First we have to determine how many iterations are needed before the Gibbs sampler has reached its stationary distribution. This is known as the \textit{burn-in} period. There are many diagnostic tests available to help determine how long this is but, in general, the most effective method is simply to look at a trace plot of the posterior sample and detect the point after which the realisations look to be from the same distribution. Figure~\ref{fig:burnin} illustrates typical output from a Gibbs sampler. Here we see the output when using three different starting points. Initially there is a transient stage and then the distribution of the output becomes the same for each chain (perhaps after iteration 500). \begin{figure}[t!] \begin{center} \includegraphics[width=0.9\textwidth]{burnin} \caption{Demonstration of burn-in using a Gibbs sampler and three initial values} \label{fig:burnin} \end{center} \end{figure} Once the Gibbs sampler has reached its stationary distribution, all subsequent iterates are realisations from the posterior distribution. Suppose that $\thetavec=(\mu,\tau)^T$ and we have run the Gibbs sampler for $N$ iterations after convergence giving a posterior sample \[ \{(\mu^{(1)},\tau^{(1)}),(\mu^{(2)},\tau^{(2)}),\ldots,(\mu^{(N)},\tau^{(N)})\}. \] We can use this sample to calculate any features of the posterior distribution. For example, we can estimate the marginal posterior densities $\pi(\mu|\xvec)$ and $\pi(\tau|\xvec)$ by using histograms of the $\mu^{(j)}$ and of the $\tau^{(j)}$ respectively. Of course, as $N$ is finite we cannot determine these densities exactly. We can also estimate other features of the posterior distribution such as the posterior means, variances and correlation by using their equivalents in the posterior sample: $\bar\mu$, $\bar\tau$, $s_\mu^2$, $s_\tau^2$ and $r_{\mu\tau}$. Again these estimates will not be exact as $N$ is finite. However, we can make them as accurate as we want by taking a sufficiently large posterior sample, that is, by taking $N$ large enough. \subsubsection*{Dealing with autocorrelation} Another problem with the output from a Gibbs sampler (after convergence) is that it is not a random sample. It should not be surprising that successive realisations from the sampler are autocorrelated, after all the output is a realisation from a Markov chain! To understand the dependence structure, we look at the sample autocorrelation function for each variable. For example, the sample autocorrelation function (ACF) for~$\mu$ at lag~$k$ is \[ r(k)=Corr(\mu^{(j)},\mu^{(j+k)}). \] Looking at a plot of this ACF can give an idea as to how much to \textit{thin} the output before it becomes un-autocorrelated (the sample autocorrelations at lags 1,2,3,\ldots are small). Thinning here means not taking every realisation by say taking say every $m$th realisation. In general, an appropriate level of thinning is determined by the largest lag $m$ at which any of the variables have a non-negligible autocorrelation. If doing this leaves a (thinned) posterior sample which is too small then the original Gibbs sampler should be re-run (after convergence) for a sufficiently large number of iterations until the thinned sample is of the required size. To get a clearer idea of how thinning works, we now look at some output from a moving average $MA(1)$ process. Figure~\ref{fig:thin1} shows this output together with its sample autocorrelation function. Theoretically output from this process should have zero autocorrelations at lags greater than one, and this is what we see (with sample noise) in the figure. If we now thin the output by taking every other value then it's clear that the autocorrelations at non-zero lags should be zero. The lower two plots show the thinned output and its sample autocorrelation function. This ACF plot suggests that the thinned output is not autocorrelated. We now look at some output from an autoregressive $AR(1)$ process. Figure~\ref{fig:thin2} shows this output together with its sample autocorrelation function. As mentioned previously, theoretically output from this process has autocorrelations which decrease geometrically, and this is what we see (with sample noise) in the figure. The smallest lag after which the autocorrelations are negligible is around 7--10. If we now thin the output by taking every 10th value then we should get output which is almost un-autocorrelated. The lower two plots show the thinned output and its sample autocorrelation function and the ACF plot is consistent with the thinned output being un-autocorrelated. If the MCMC output is un-autocorrelated then the accuracy of $\bar\mu$ is roughly $\pm 2s_\mu/\sqrt{N}$. However if the Markov chain followed an autoregressive $AR(1)$ process, its autocorrelations would decrease geometrically, with $r(k)\simeq r(1)^k$, $k\geq 2$. In this case it can be shown that the accuracy of $\bar\mu$ is roughly $\pm 2s_\mu/\sqrt{N\{1-r(1)\}^2}$, that is, because the process has autocorrelation, the amount of information in the data is equivalent to a random sample with size $N_{eff}=N\{1-r(1)\}^2$. This effective random sample size calculation gets more complicated for processes with non-zero higher order autocorrelations and this is why we usually adopt the simplistic method of thinning. It's worth noting that, in general, MCMC output with positive autocorrelations has $N_{eff}N$. \newpage \begin{figure}[!h] \begin{center} \includegraphics[width=0.9\textwidth]{ma1-trace-acf} \includegraphics[width=0.9\textwidth]{ma1thin-trace-acf} \caption{Effect of thinning output from a $MA(1)$ process} \label{fig:thin1} \end{center} \end{figure} \newpage \begin{figure}[!h] \begin{center} \includegraphics[width=0.9\textwidth]{ar1-trace-acf} \includegraphics[width=0.9\textwidth]{ar1thin-trace-acf} \caption{Effect of thinning output from a $AR(1)$ process} \label{fig:thin2} \end{center} \end{figure} \newpage \subsubsection*{Strategy} \begin{enumerate} \item Determine the \textit{burn-in} period, after which the Gibbs sampler has reached its stationary distribution. This may involve thinning the posterior sample as slowly snaking trace plots may be due to high autocorrelations rather than a lack of convergence. \item After this, determine the level of thinning needed to obtain a posterior sample whose autocorrelations are roughly zero. \item Repeat steps 1 and 2 several times using different initial values to make sure that the sample really is from the stationary distribution of the chain, that is, from the posterior distribution. \end{enumerate} \subsection{Bayesian inference using a Gibbs sampler} \begin{example}~\\ Construct a Gibbs sampler for the posterior distribution in Example~\ref{ex:semiconjnormal}. \label{ex:semiconj} \end{example} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \textcolor{blue}{ %\begin{soln} We first need to determine the conditional posterior density for $\mu$ and for~$\tau$.\\ The conditional posterior density for $\mu$ is, for $\mu\in\mathbb{R}$ \begin{align*} \pi(\mu|\tau,\xvec)&=\dfrac{\pi(\mu,\tau|\xvec)}{\pi(\tau|\xvec)}\\ &\propto\pi(\mu,\tau|\xvec)\\ &\propto\tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\\ &\propto\exp\left\{-\frac{c}{2}(\mu-b)^2 -\frac{n\tau}{2}(\bar x-\mu)^2\right\} \intertext{after moving all multiplicative constants not involving $\mu$ into the proportionality sign, and so after completing the square in the exponent} \pi(\mu|\tau,\xvec)&\propto\exp\left\{-\frac{c+n\tau}{2} \left(\mu-\frac{bc+n\tau\bar x}{c+n\tau}\right)^2\right\},\mu\in\mathbb{R} \intertext{that is} \mu|\tau,\xvec&\sim N\left(\frac{bc+n\tau\bar x}{c+n\tau},\frac{1}{c+n\tau}\right). \end{align*} \vspace{2cm} }}} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \textit{Solution to Example 4.2 continued...} \vspace{0.7cm} \textcolor{blue}{ Note that we could have obtained this result by using Example~\ref{ex:normal} with $d=c$.\\ The conditional posterior density for $\tau$ is, for $\tau>0$ \begin{align*} \pi(\tau|\mu,\xvec)&=\dfrac{\pi(\mu,\tau|\xvec)}{\pi(\mu|\xvec)}\\ &\propto\pi(\mu,\tau|\xvec)\\ &\propto\tau^{g+\frac{n}{2}-1}\exp\left\{-\frac{c}{2}(\mu-b)^2-h\tau -\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}\\ &\propto\tau^{g+\frac{n}{2}-1}\exp\left\{-\left[h +\frac{n}{2}\left\{s^2+(\bar x-\mu)^2\right\}\right]\tau\right\}\\ \intertext{after moving all multiplicative constants not involving $\tau$ into the proportionality sign, and so} \tau|\mu,\xvec&\sim Ga\left(g+\frac{n}{2},h +\frac{n}{2}\left\{s^2+(\bar x-\mu)^2\right\}\right). \end{align*} We also need to initialise the algorithm. We might use the prior mean ($\mu^{(0)}=b$, $\tau^{(0)}=g/h$) or the mle ($\mu^{(0)}=\bar x$, $\tau^{(0)}=1/s^2$). Alternatively, if we wanted to compare different runs of the Gibbs sampler, we might simulate $(\mu^{(0)},\tau^{(0)})$ from the prior distribution.\\ Therefore the Gibbs sampler is \begin{enumerate} \item Initialise the iteration counter to $j=1$.\\ Initialise the state of the chain by taking $(\mu^{(0)},\tau^{(0)})$ as $(b,g/h)$ or $(\bar x,1/s^2)$ or as a random draw from the prior distribution. \item Obtain new values $\mu^{(j)}$ and $\tau^{(j)}$ from $\mu^{(j-1)}$ and $\tau^{(j-1)}$ via \begin{align*} \mu^{(j)} &\sim N\left(\frac{bc+n\tau^{(j-1)}\bar x}{c+n\tau^{(j-1)}}, \frac{1}{c+n\tau^{(j-1)}}\right)\\ \tau^{(j)} &\sim Ga\left(g+\frac{n}{2},h +\frac{n}{2}\left\{s^2+\left(\bar x-\mu^{(j)}\right)^2\right\}\right) \end{align*} \item Change counter $j$ to $j+1$, and return to step 2. \end{enumerate} \vspace{6cm} }}} %\end{soln} \subsubsection*{Comments} Notice that, since $\mu$ and $\tau$ independent \textit{a priori}, $\mu|\tau\sim N(b,1/c)$. Therefore, given $\tau$, the normal prior for $\mu$ is conjugate. Similarly, $\tau|\mu\sim Ga(g,h)$ and so, given $\mu$, the gamma prior for $\tau$ is conjugate. Therefore, both conditional priors (for $\mu|\tau$ and $\tau|\mu$) are conjugate. Such priors are called \textit{semi-conjugate}. \subsubsection*{Producing and analysing output from this Gibbs sampler} The \texttt{R} function \texttt{gibbsNormal} in the library \texttt{nclbayes} implements this Gibbs sampling algorithm. The library also contains the functions \texttt{mcmcProcess} which can be used to remove the burn--in and thin the output, and \texttt{mcmcAnalysis} which analyses the MCMC output. Let us consider the case in which the data have size $n=100$, mean $\bar x=15$ and standard deviation $s=4.5$ and the prior distribution has $\mu\sim N(10,1/100)$ and $\tau\sim Ga(3,12)$, independently. The following code produces output from this Gibbs sampler, initialising at $(10,0.25)$ and then analyses the output. %\hline \begin{verbatim} library(nclbayes) posterior=gibbsNormal(N=1000,initial=c(10,0.25), priorparam=c(10,1/100,3,12),n=100,xbar=15,s=4.5) posterior2=mcmcProcess(input=posterior,burnin=10,thin=1) op=par(mfrow=c(2,2)) plot(posterior,col=c(1:length(posterior)),main="All realisations") plot(posterior,type="l",main="All realisations") plot(posterior2,col=c(1:length(posterior2)),main="After deleting first 10") plot(posterior2,type="l",main="After deleting first 10") par(op) mcmcAnalysis(posterior,rows=2,show=F) mcmcAnalysis(posterior2,rows=2,show=F) \end{verbatim} %\hline The first block of code runs the function \texttt{gibbsNormal}, with initial values \texttt{initial} taken as the prior means, to obtain the output. The next block then uses the function \texttt{mcmcProcess} to post-process the Gibbs output by deleting an initial $\texttt{burnin}=10$ values and then not thinning by taking $\texttt{thin}=1$. The next block of code produces the plots in Figure~\ref{fig:gibbsNormal1}. After this the code analyses the Gibbs sampler output and produces the plots in Figure~\ref{fig:gibbsNormal2}. In Figure~\ref{fig:gibbsNormal1}, the top left plot shows the values produced by the Gibbs sampler: notice the initial value $(\mu^{(0)},\tau^{(0)})=(b=10,g/h=0.25)$ appears at the top left part of the plot and the other values towards the bottom right part (in different colours). The top right plot is another representation of this output but here each consecutive pair $(\mu^{(j)},\tau^{(j)})$ are joined by a line. This clearly shows that all pairs after the first one remain in the same vicinity (the bottom right part). The lower plots are the equivalent ones to the upper plots but only use the Gibbs sampler output after deleting the first 10 pairs. \begin{figure}[t!] \begin{center} \includegraphics[width=0.9\textwidth]{gibbsNormal1} \caption{Progress of the MCMC scheme} \label{fig:gibbsNormal1} \end{center} \end{figure} In Figure~\ref{fig:gibbsNormal2}, the top two rows of plots summarise the Gibbs sampler output using all realisations and the bottom two rows of plots are the equivalent plots but only use the Gibbs sampler output after deleting the first 10 pairs. The first column of plots shows the trace plot of the output, that is, the values for each variable ($\mu$ and $\tau$) as the sampler iterates from its first value to its final value. The top two first column plots clearly show the initial value $(\mu^{(0)},\tau^{(0)})=(b=10,g/h=0.25)$, after which the subsequent values all look to be from the same distribution. In particular, there looks to be no change in the range of values or in the mean value. The bottom two first column plots emphasise these points; here the first 10 values have been deleted as burn-in, though probably we needed only to delete the first 2 or 3 values. Note that the benefit of using \texttt{R} code that runs quickly is that adopting a conservative strategy which deletes too many values as burn-in, does not have significant time implications (though this is not a sensible strategy if the code is very slow and takes months to run!). \begin{figure}[p!] \begin{center} \includegraphics[width=0.9\textwidth]{gibbsNormal2} \bigskip\hrule\bigskip\bigskip \includegraphics[width=0.9\textwidth]{gibbsNormal3} \caption{Trace plots, autocorrelation plots and histograms of the Gibbs sampler output. Upper plots: all realisations. Lower plots: after deleting the first 10 iterations} \label{fig:gibbsNormal2} \end{center} \end{figure} The second column shows the autocorrelation function for each variable. Note that the spike at lag~0 is due to $r(0)=Corr(\mu^{(j)},\mu^{(j)})=1$. The plots show that the autocorrelations at all other lags are negligible, and so no thinning is needed. The final column shows histograms of the Gibbs sampler output. If using a burn-in of 10 iterations is okay (and it is here!) then the subsequent output can be taken as our posterior sample and therefore the lower two histograms will be good estimates of the marginal densities: good because the output is (almost) uncorrelated and the sample size is quite large. \newpage If we use the command: \begin{verbatim} mcmcAnalysis(posterior2,rows=2) \end{verbatim} that is, don't use the \texttt{show=F} option, then the function will produce the plots and also various useful numerical summaries. In this run of the Gibbs sampler it gave: \begin{verbatim} N = 990 iterations mu tau Min. :13.66 Min. :0.03025 1st Qu.:14.68 1st Qu.:0.04663 Median :15.01 Median :0.05068 Mean :14.99 Mean :0.05110 3rd Qu.:15.29 3rd Qu.:0.05557 Max. :16.41 Max. :0.07515 Standard deviations: mu tau 0.448562708 0.006743413 \end{verbatim} We can calculate other features of the joint posterior distribution such as its correlation \[ Corr(\mu,\tau|\xvec)=-0.002706 \] using the command \texttt{cor(posterior2)}, and summarise the posterior sample with a plot of its values and its marginal distributions; see Figure~\ref{fig:gibbsNormal3}. We can also determine $100(1-\alpha)\%$ equi-tailed confidence intervals as follows. Suppose we have $N$ realisations from our Gibbs sampler. If we sort the values into increasing order then the confidence interval will have end points which are $N\alpha/2$th and $N(1-\alpha/2)$th values. The \texttt{nclbayes} package has a function \texttt{mcmcCi} to do this. In this case we would use \texttt{mcmcCi(posterior2,level=0.95)} and, for this output, obtain the 95\% confidence intervals as \begin{align*} \mu&:\quad (14.071,~15.803)\\ \tau&:\quad (0.03818,~0.06499). \end{align*} \begin{figure}[t!] \begin{center} \includegraphics[width=0.5\textwidth]{gibbsNormal3a}\hspace*{-1mm} \includegraphics[width=0.5\textwidth]{gibbsNormal3b} \caption{Plot of the bivariate posterior sample and their marginal distributions. Left plot: $(\mu,\tau)^T$; right plot: $(\mu,\sigma)^T$.} \label{fig:gibbsNormal3} \end{center} \end{figure} Now we have a sample from the posterior distribution, we can determine the posterior distribution for any function of the parameters. For example, if we want the posterior distribution for $\sigma=1/\sqrt\tau$ then we can easily obtain realisations of $\sigma$ as $\sigma^{(j)}=1/\sqrt{\tau^{(j)}}$, from which we can produce a plot of its values and its marginal distributions (see Figure~\ref{fig:gibbsNormal3}) and also obtain its numerical summaries: \begin{verbatim} > sigma=1/sqrt(posterior2[,2]) > summary(sigma) Min. 1st Qu. Median Mean 3rd Qu. Max. 3.648 4.242 4.442 4.453 4.631 5.750 > sd(sigma) [1] 0.2977992 > quantile(sigma,probs=c(0.025,0.975)) 2.5% 97.5% 3.922480 5.109632 \end{verbatim} \subsubsection*{Summary} We can use the (converged and thinned) MCMC output to do the following. \begin{itemize} \item Obtain the posterior distribution for any (joint) functions of the parameters, such as $\sigma=1/\sqrt\tau$ or $(\theta_1=\mu-\tau,\theta_2=e^{\mu+\tau/2})^T$ \item Look at bivariate posterior distributions via scatter plots \item Look at univariate marginal posterior distributions via histograms or boxplots \item Obtain numerical summaries such as the mean, standard deviation and confidence intervals for single variables and correlations between variables. \end{itemize} \newpage \newpage \begin{example}~\\ Suppose we have a random sample of size~$n$ from a gamma \label{ex:gammars} $Ga(\alpha,\lambda)$ distribution in which both the index~$\alpha>0$ and scale parameter~$\lambda>0$ are unknown, that is, $X_i|\alpha,\lambda\sim Ga(\alpha,\lambda)$, $i=1,2,\ldots,n$ (independent). We shall assume independent prior distributions for these parameters, with $\alpha\sim Ga(a,b)$ and $\lambda\sim Ga(c,d)$ for known values $a$, $b$, $c$ and $d$. Determine the posterior density for $(\alpha,\lambda)^T$ and hence the posterior conditional densities for $\alpha|\lambda$ and $\lambda|\alpha$. \end{example} %\newnotes %\begin{soln} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \textcolor{blue}{ The likelihood function is \begin{align*} f(\xvec|\alpha,\lambda)&=\prod_{i=1}^n \dfrac{\lambda^\alpha x_i^{\alpha-1}e^{-\lambda x_i}}{\Gamma(\alpha)}\\ &=\dfrac{\lambda^{n\alpha}\left(\prod x_i\right)^{\alpha-1}e^{-n\bar x\lambda}}{\Gamma(\alpha)^n}\\ &=\dfrac{\lambda^{n\alpha}\bar x_g^{n(\alpha-1)}e^{-n\bar x\lambda}}{\Gamma(\alpha)^n} \intertext{where $\bar x_g=\sqrt[n]{\prod x_i}$ is the geometric mean of the data. Using Bayes Theorem, the posterior density is} \pi(\alpha,\lambda|\xvec)&\propto\pi(\alpha,\lambda)\,f(\xvec|\alpha,\lambda) \intertext{and so, for $\alpha>0$, $\lambda>0$} \pi(\alpha,\lambda|\xvec)&\propto\alpha^{a-1}e^{-b\alpha}\times \lambda^{c-1}e^{-d\lambda}\times \dfrac{\lambda^{n\alpha}\bar x_g^{n(\alpha-1)}e^{-n\bar x\lambda}}{\Gamma(\alpha)^n}\\ &\propto\dfrac{\alpha^{a-1}\bar x_g^{n\alpha} \lambda^{c+n\alpha-1}}{\Gamma(\alpha)^n}\exp\left\{-b\alpha-(d+n\bar x)\lambda\right\}. \intertext{This is not the density of a standard distribution. The conditional posterior density for $\alpha$ is, for $\alpha>0$} \pi(\alpha|\lambda,\xvec)&=\dfrac{\pi(\alpha,\lambda|\xvec)}{\pi(\lambda|\xvec)}\\ &\propto\pi(\alpha,\lambda|\xvec)\\ &\propto\dfrac{\alpha^{a-1}\bar x_g^{n\alpha} \lambda^{c+n\alpha-1}}{\Gamma(\alpha)^n}\exp\left\{-b\alpha-(d+n\bar x)\lambda\right\}\\ &\propto\dfrac{\alpha^{a-1}\bar x_g^{n\alpha} \lambda^{n\alpha}e^{-b\alpha}}{\Gamma(\alpha)^n}\\ &\propto\dfrac{\alpha^{a-1}e^{-(b-n\log \bar x_g-n\log\lambda)\alpha}}{\Gamma(\alpha)^n}. \end{align*} }}} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \textit{Solution to Example 4.2...} \vspace{0.7cm} \textcolor{blue}{ This looks like it might be a gamma density but (i) we would require $b-n\log \bar x_g-n\log\lambda>0$ which can't be guaranteed (ii) the divisor term is wrong: it would have to be $\Gamma(a)$ for this to be a gamma density.\\ The conditional posterior density for $\lambda$ is, for $\lambda>0$ \begin{align*} \pi(\lambda|\alpha,\xvec)&=\dfrac{\pi(\alpha,\lambda|\xvec)}{\pi(\alpha|\xvec)}\\ &\propto\pi(\alpha,\lambda|\xvec)\\ &\propto\dfrac{\alpha^{a-1}\bar x_g^{n\alpha} \lambda^{c+n\alpha-1}}{\Gamma(\alpha)^n}\exp\left\{-b\alpha-(d+n\bar x)\lambda\right\}\\ &\propto\lambda^{c+n\alpha-1}\,e^{-(d+n\bar x)\lambda}, \end{align*} that is, $\lambda|\alpha,\xvec\sim Ga(c+n\alpha,\,d+n\bar x)$. %\end{soln} \vspace{11cm} }}} %%Hack for vspace %\vspace{-3cm} The full conditional distribution (FCD) for $\lambda$ is a standard distribution and so it is straightforward to simulate from this distribution. However, the FCD for $\alpha$ is not a standard distribution and not easy to simulate from. Therefore we cannot use a Gibbs sampler to simulate from the posterior $\pi(\alpha,\lambda|\xvec)$. Fortunately, there are other methods available to help us in these situations. \newpage \section{Metropolis-Hastings sampling} The Gibbs sampler is a very powerful tool but is only useful if the full conditional distributions (FCDs) are standard distributions (which are easy to simulate from). Fortunately there is a class of methods which can be used when the FCDs are non-standard. These methods are known as Metropolis-Hastings schemes. Suppose we want to simulate realisations from the posterior density $\pi(\thetavec|\xvec)$ and all of the FCDs are non-standard. Suppose further that we have a \emph{proposal distribution} with density $q(\thetavec^*|\thetavec)$, which is easy to simulate from. This distribution gives us a way of proposing new values~$\thetavec^*$ from the current value~$\thetavec$. Consider the following algorithm: \enlargethispage*{1cm} \begin{enumerate} \item Initialise the iteration counter to $j=1$, and initialise the chain to $\thetavec^{(0)}$. \item Generate a \emph{proposed} value $\thetavec^*$ using the proposal distribution $q(\thetavec^*|\thetavec^{(j-1)})$. \item Evaluate the \emph{acceptance probability} $\alpha(\thetavec^{(j-1)},\thetavec^*)$ of the proposed move, where \[ \alpha(\thetavec,\thetavec^*) = \min\left\{1, \frac{\pi(\thetavec^*|\xvec)\,q(\thetavec|\thetavec^*)}{\pi(\thetavec|\xvec)\,q(\thetavec^*|\thetavec)} \right\}. \] \item Set $\thetavec^{(j)}=\thetavec^*$ with probability $\alpha(\thetavec^{(j-1)},\thetavec^*)$, and set $\thetavec^{(j)}=\thetavec^{(j-1)}$ otherwise. \item Change the counter from $j$ to $j+1$ and return to step 2. \end{enumerate} In other words, at each stage, a new value is generated from the proposal distribution. This is either accepted, in which case the chain moves, or rejected, in which case the chain stays where it is. Whether or not the move is accepted or rejected depends on an acceptance probability which itself depends on the relationship between the density of interest and the proposal distribution. Note that the posterior density $\pi(\cdot|\xvec)$ only enters into the acceptance probability as a ratio, and so the method can be used when it is known up to a scaling constant, that is, \begin{align*} \alpha(\thetavec,\thetavec^*) &= \min\left\{1, \frac{\pi(\thetavec^*)\,f(\xvec|\thetavec^*)\,q(\thetavec|\thetavec^*)}{\pi(\thetavec)\,f(\xvec|\thetavec)\,q(\thetavec^*|\thetavec)} \right\}, \intertext{since} \pi(\thetavec|\xvec)&=\dfrac{\pi(\thetavec)\,f(\xvec|\thetavec)}{f(\xvec)}. \end{align*} It can be shown that the above algorithm defines a Markov chain with $\pi(\thetavec|\xvec)$ as its stationary distribution. Notice that the above description holds for all possible proposal distributions (subject to them generating realisations from the full parameter space). But are some choices better than others? We now discuss some commonly used proposal distributions. \subsection{Symmetric chains (Metropolis method)} The simplest case is the Metropolis sampler and uses a symmetric proposal distribution, that is, one with $q(\thetavec^*|\thetavec)=q(\thetavec|\thetavec^*)$, $\forall~\thetavec,\thetavec^*$. In this case the acceptance probability simplifies to \[ \alpha(\thetavec,\thetavec^*) = \min\left\{ 1, \frac{\pi(\thetavec^*|\xvec)}{\pi(\thetavec|\xvec)} \right\}, \] and hence does not involve the proposal density at all. Consequently proposed moves which will take the chain to a region of higher posterior density are always accepted, while moves which take the chain to a region of lower posterior density are accepted with probability proportional to the ratio of the two densities --- moves which will take the chain to a region of very low density will be accepted with very low probability. Note that any proposal of the form $q(\thetavec^*|\thetavec)=f(|\thetavec^*-\thetavec|)$ is symmetric, where $f(\cdot)$ is some zero mean density function, as $|\thetavec^*-\thetavec|=|\thetavec-\thetavec^*|$. In this case, the proposal value is a symmetric displacement from the current value. This motivates the following. \subsubsection{Random walk proposals} Consider the random walk proposal in which the proposed value $\thetavec^*$ depends on the current value~$\thetavec$ via \[ \thetavec^*=\thetavec+\vec{w}, \] where $\vec{w}$ is a random $p\times 1$ vector from the zero mean density $f(\cdot)$ which is \textit{symmetric} about its mean, and is independent of the state of the chain. We can generate our \emph{proposal} value by first simulating an \emph{innovation} $\vec{w}$, and then set the \emph{proposal} value to $\thetavec^*=\thetavec+\vec{w}$. Clearly $q(\thetavec^*|\thetavec)=f(\thetavec^*-\thetavec)=f(\vec{w})$. Also $\thetavec=\thetavec^*-\vec{w}$ and so $q(\thetavec|\thetavec^*)=f(\thetavec-\thetavec^*)=f(-\vec{w})$. However, as $f(\cdot)$ is a zero mean symmetric density, we have that $f(\vec{w})=f(-\vec{w})$ and so $q(\thetavec^*|\thetavec)=q(\thetavec|\thetavec^*)$. But, what distribution should we use for~$f(\cdot)$? A distribution which is simple and easy to simulate from would be good, with obvious choices of the uniform or normal distributions, though the normal distribution is generally better, but is a bit more expensive to simulate. What variance should we use for the distribution we choose? This choice will affect the acceptance probability, and hence the overall proportion of accepted moves. If the variance of the innovation is too low, then most proposed values will be accepted, but the chain will move very slowly around the space --- the chain is said to be too ``cold''. On the other hand, if the variance of the innovation is too large, very few proposed values will be accepted, but when they are, they will often correspond to quite large moves --- the chain is said to be too ``hot''. Theoretically it has been shown that the optimal acceptance rate is around 0.234 --- this is an asymptotic result (for large samples of data) --- but experience suggests that an acceptance rate of around 20--30\% is okay. Thus, the variance of the innovation distribution should be ``tuned'' to get an acceptance rate of around this level. \subsubsection{Normal random walk proposals} A symmetric normal random walk proposal takes the form $\theta^*|\theta\sim N(\theta,k^2)$ for some innovation size~$k>0$. This is a symmetric proposal because $\theta^*=\theta+w$, where $w\sim N(0,k^2)$ has a density which is symmetric about zero. Also the proposal ratio is %\begin{ggap} \[ \frac{q(\theta|\theta^*)}{q(\theta^*|\theta)} =\frac{\dfrac{1}{\sqrt{2\pi k^2}}\exp\left\{-\dfrac{(\theta-\theta^*)^2}{2k^2}\right\}} {\dfrac{1}{\sqrt{2\pi k^2}}\exp\left\{-\dfrac{(\theta^*-\theta)^2}{2k^2}\right\}} =1. \] %\end{ggap} \subsubsection{Uniform random walk proposals} A symmetric uniform random walk proposal takes the form $\theta^*|\theta\sim U(\theta-a,\theta+a)$ for some innovation size~$a>0$. This is a symmetric proposal because $\theta^*=\theta+w$, where $w\sim U(-a,a)$ has a density which is symmetric about zero. Also %\begin{ggap} \[ \frac{q(\theta|\theta^*)}{q(\theta^*|\theta)} =\frac{1/(2a)}{1/(2a)} =1. \] %\end{ggap} \section{Hybrid methods} We have now seen how we can use the Gibbs sampler to sample from multivariate distributions provided that we can simulate from the full conditional distributions. We have also seen how we can use Metropolis-Hastings methods to sample from awkward FCDs. If we wish, we can combine these in order to form hybrid Markov chains whose stationary distribution is a distribution of interest. \subsection{Componentwise transitions} Given a posterior distribution with full conditional distributions that are awkward to sample from directly, we can define a Metropolis-Hastings scheme for each full conditional distribution, and apply them to each component in turn for each iteration. This is like the Gibbs sampler, but each component update is a Metropolis-Hastings update, rather than a direct simulation from the full conditional distribution. Each of these steps will require its own proposal distribution. The algorithm is as follows: \begin{enumerate} \item Initialise the iteration counter to $j=1$.\\ Initialise the state of the chain to $\thetavec^{(0)}=(\theta_1^{(0)},\ldots,\theta_p^{(0)})^T$. \item Let $\thetavec_{-i}^{(j)}=\bigl(\theta_1^{(j)},\ldots,\theta_{i-1}^{(j)},\theta_{i+1}^{(j-1)},\ldots,\theta_p^{(j-1)}\bigr)^T$, $i=1,2,\ldots,p$.\\ Obtain a new value $\thetavec^{(j)}$ from $\thetavec^{(j-1)}$ by successive generation of values \begin{itemize} \item $\theta_1^{(j)}\sim \pi(\theta_1|\thetavec_{-1}^{(j)},\xvec)$ using a Metropolis--Hastings step with proposal distribution $q_1(\theta_1|\theta_1^{(j-1)},\thetavec_{-1}^{(j)})$ \item $\theta_2^{(j)}\sim \pi(\theta_2|\thetavec_{-2}^{(j)},\xvec)$ using a Metropolis--Hastings step with proposal distribution $q_2(\theta_2|\theta_2^{(j-1)},\thetavec_{-2}^{(j)})$ \item[] $\vdots$ \item $\theta_p^{(j)}\sim \pi(\theta_p|\thetavec_{-p}^{(j)},\xvec)$ using a Metropolis--Hastings step with proposal distribution $q_p(\theta_p|\theta_p^{(j-1)},\thetavec_{-p}^{(j)})$ \end{itemize} \item Change counter $j$ to $j+1$, and return to step 2. \end{enumerate} This is in fact the original form of the Metropolis algorithm. Note that the distributions $\pi(\theta_i|\thetavec_{-i}^{(j)},\xvec)$ are just the FCDs. \subsection{Metropolis within Gibbs} Given a posterior distribution with full conditional distributions, some of which may be simulated from directly, and others which have Metropolis-Hastings updating schemes, the Metropolis within Gibbs algorithm goes through each in turn, and simulates directly from the full conditional, or carries out a Metropolis-Hastings update as necessary. This algorithm is, in fact, just the above algorithm but uses the full conditional distributions as the proposal distributions when they are easy to simulate from. To see this, suppose that we can simulate from the FCD $\pi(\theta_i|\thetavec_{-i}^{(j)},\xvec)$ and use this as the proposal distribution, that is, take $\theta_i^*\sim\pi(\theta_i|\thetavec_{-i}^{(j)},\xvec)$. Then the acceptance probability for this step is \begin{align*} \alpha(\theta_i,\theta_i^*) &=\min\left\{1, \frac{\pi(\theta_i^*|\thetavec_{-i}^{(j)},\xvec)}{\pi(\theta_i|\thetavec_{-i}^{(j)},\xvec)} \frac{q(\theta_i|\theta_i^*,\thetavec_{-i}^{(j)})}{q(\theta_i^*|\theta_i,\thetavec_{-i}^{(j)})}\right\}\\ &=\min\left\{1, \frac{\pi(\theta_i^*|\thetavec_{-i}^{(j)},\xvec)}{\pi(\theta_i|\thetavec_{-i}^{(j)},\xvec)} \frac{\pi(\theta_i|\thetavec_{-i}^{(j)},\xvec)}{\pi(\theta_i^*|\thetavec_{-i}^{(j)},\xvec)}\right\}\\ &=\min(1,1)\\ &=1, \intertext{that is, we always accept the proposal from the FCD.} \end{align*} \newpage \begin{example}~\\ Construct an MCMC scheme for the problem in Example~\ref{ex:gammars} where we had a random sample of size~$n$ from a gamma $Ga(\alpha,\lambda)$ distribution and independent gamma $Ga(a,b)$ and $Ga(c,d)$ prior distributions for $\alpha$ and $\lambda$ respectively. Recall that the FCDs were \begin{align*} \pi(\alpha|\lambda,\xvec) &\propto\dfrac{\alpha^{a-1}e^{(-b+n\log \bar x_g+n\log\lambda)\alpha}} {\Gamma(\alpha)^n},\quad\alpha>0 \intertext{and} \pi(\lambda|\alpha,\xvec) &\propto\lambda^{c+n\alpha-1}\,e^{-(d+n\bar x)\lambda},\quad\lambda>0. \end{align*} \end{example} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \textcolor{blue}{We have that $\lambda|\alpha,\xvec\sim Ga(c+n\alpha,\,d+n\bar x)$ but the distribution of $\alpha|\lambda,\xvec$ is non-standard. Therefore we will use a Metropolis within Gibbs algorithm which uses a Gibbs step for $\lambda$ and a M--H step for $\alpha$. In the M--H step, we will use a normal random walk proposal distribution, with $\alpha^*|\alpha\sim N(\alpha,\Sigma_\alpha)$, in which the proposal $\alpha^*$ has acceptance probability $\min(1,A)$, where \begin{align*} A&=\frac{\pi(\alpha^*|\lambda,\xvec)}{\pi(\alpha|\lambda,\xvec)}\\ &=\frac{{\alpha^*}^{a-1}e^{(-b+n\log \bar x_g+n\log\lambda)\alpha^*}} {\Gamma(\alpha^*)^n}\times \frac{\Gamma(\alpha)^n}{\alpha^{a-1}e^{(-b+n\log \bar x_g+n\log\lambda)\alpha}},\quad\alpha^*>0\\ &=\left(\frac{\alpha^*}{\alpha}\right)^{a-1}\, \left\{\frac{\Gamma(\alpha)}{\Gamma(\alpha^*)}\right\}^n\, e^{(-b+n\log \bar x_g+n\log\lambda)(\alpha^*-\alpha)},\quad\alpha^*>0, \end{align*} and zero otherwise. \vspace{8cm} }}} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \textit{Solution to Example 4.4...} \vspace{0.7cm} \textcolor{blue}{ We need to find a value for $\alpha$ or $\lambda$ to initialise the MCMC algorithm. We could use a value simulated from the prior or use the prior mean. Alternatively we could use the mle but this is rather complicated to determine. However, the moment estimates are rather more straightforward, and equating first and second population and sample moments gives \begin{align*} \text{mean}:\quad\frac{\tilde\alpha}{\tilde\lambda}&=\bar x & \text{variance}:&\quad\frac{\tilde\alpha}{\tilde\lambda^2}=s^2 \intertext{and so} \tilde\lambda&=\frac{\bar x}{s^2} & \tilde\alpha&=\frac{\bar x^2}{s^2}. \end{align*} Therefore the MCMC algorithm is \begin{enumerate} \item Initialise the iteration counter to $j=1$, and initialise the chain to $\alpha^{(0)}=(\bar x/s)^2$. \item Obtain a new value \[ \lambda^{(j)}\sim Ga(c+n\alpha^{(j-1)},\,d+n\bar x) \] \item Generate a proposed value $\alpha^*\sim N(\alpha^{(j-1)},\Sigma_\alpha)$ \item Evaluate the acceptance probability $\min(1,A)$ at $\alpha^*=\alpha^*$, $\alpha=\alpha^{(j-1)}$ and $\lambda=\lambda^{(j)}$ \item Set $\alpha^{(j)}=\alpha^*$ with probability $\min(1,A)$, and set $\alpha^{(j)}=\alpha^{(j-1)}$ otherwise. \item Change the counter from $j$ to $j+1$ and return to step 2. \end{enumerate} \vspace{1.15cm} %\end{soln} }}} The \texttt{R} function \texttt{mwgGamma} in the library \texttt{nclbayes} implements this Metropolis within Gibbs algorithm. The following code produces posterior output from an analysis of a dataset with $n=50$, $\bar x=0.62$, $\bar x_g=0.46$ and $s=0.4$, with prior beliefs represented by $a=2$, $b=1$, $c=3$ and $d=1$, and uses a normal random walk proposal with variance $\Sigma_\alpha=0.9^2$ as this gives a reasonable acceptance probability of~0.237. The initial value is taken as the moment estimate $\tilde\alpha=(\bar x/s)^2$. \begin{verbatim} posterior=mwgGamma(N=20000,initial=(0.46/0.4)^2,innov=0.9, priorparam=c(2,1,3,1),n=50,xbar=0.62,xgbar=0.46,show=TRUE) mcmcAnalysis(posterior,rows=2,show=F) posterior2=mcmcProcess(input=posterior,burnin=10,thin=20) mcmcAnalysis(posterior2,rows=2) \end{verbatim} The upper plots in Figure~\ref{fig:gammars} show all the output of this MCMC scheme and the lower plots show the output after deleting the first 10~iterations as burn--in and then thinning by only taking every 20th iterate to reduce the autocorrelations. \begin{figure}[p!] \begin{center} \includegraphics[width=0.9\textwidth]{metropwithingibbs} \bigskip\hrule\bigskip\bigskip \includegraphics[width=0.9\textwidth]{metropwithingibbs2} \caption{Trace plots, autocorrelation plots and histograms of the Metropolis with Gibbs output. Upper plots: all realisations. Lower plots: with $\text{burn-in}=10$, $\text{thin}=20$.} \label{fig:gammars} \end{center} \end{figure} \newpage \subsection*{Comments} \begin{enumerate} \item If you're unsure whether the proposal distribution is symmetric then it's quite straightforward to examine the proposal ratio. In this last example, we have proposal $\alpha^*|\alpha\sim N(\alpha,\Sigma_\alpha)$ and so \[ \frac{q(\alpha|\alpha^*)}{q(\alpha^*|\alpha)} =\frac{\dfrac{1}{\sqrt{2\pi\Sigma_\alpha}}\exp\left\{-\dfrac{(\alpha-\alpha^*)^2}{2\Sigma_\alpha}\right\}} {\dfrac{1}{\sqrt{2\pi\Sigma_\alpha}}\exp\left\{-\dfrac{(\alpha^*-\alpha)^2}{2\Sigma_\alpha}\right\}} =1. \] \item A normal random walk proposal $\alpha^*$ is not accepted if it is negative as, in this case, $A=0$. This can be wasteful. An alternative is to use a proposal distribution which only generates positive proposal values, such as $\alpha^*|\alpha\sim LN(\log\alpha,\Sigma_\alpha)$. Using this skewed proposal distribution, we have \[ \frac{q(\alpha|\alpha^*)}{q(\alpha^*|\alpha)} =\frac{\dfrac{1}{\alpha\sqrt{2\pi\Sigma_\alpha}}\exp\left\{-\dfrac{(\log\alpha-\log\alpha^*)^2}{2\Sigma_\alpha}\right\}} {\dfrac{1}{\alpha^*\sqrt{2\pi\Sigma_\alpha}}\exp\left\{-\dfrac{(\log\alpha^*-\log\alpha)^2}{2\Sigma_\alpha}\right\}} =\frac{\alpha^*}{\alpha}. \] Therefore the acceptance probability for a proposed value $\alpha^*$ is $\min(1,B)$ where \begin{align*} B&=\frac{\pi(\alpha^*|\lambda,\xvec)}{\pi(\alpha|\lambda,\xvec)}\times \frac{q(\alpha|\alpha^*)}{q(\alpha^*|\alpha)} \\ &=\frac{\alpha^*\,\pi(\alpha^*|\lambda,\xvec)}{\alpha\,\pi(\alpha|\lambda,\xvec)}. \end{align*} The acceptance probability is still quite straightforward to calculate, and with this proposal distribution we never reject proposal values that are inconsistent with the parameter space (here $\alpha>0$). Incidentally, $\log X\sim N(\mu,\sigma^2)$ if $X\sim LN(\mu,\sigma^2)$ and so using a log-normal proposal is the same as using a normal random walk on the log scale, that is, a normal random walk for $\log\alpha$. Also log-normal proposals are easy to simulate because if $Y\sim N(\mu,\sigma^2)$ then $e^Y\sim LN(\mu,\sigma^2)$. \item Dealing with a constraint such as $\alpha>0$ in optimisation methods or here in MCMC methods, can be solved by re-parameterising the model. Here, for example, we could work with $A=\log\alpha$ and obtain realisations from the posterior distribution for~$A$. Once we have these realisations we can easily obtain realisations from the posterior distribution for~$\alpha=e^A$. Working in $A$ rather than $\alpha$ means we have to simulate realisations from the conditional posterior \[ \pi_A(A|\lambda,\xvec)=\pi_\alpha(e^A|\lambda,\xvec)\times \left|\frac{d}{da}\, e^A\right| =e^A\,\pi_\alpha(e^A|\lambda,\xvec) \] using \eqref{eq:trans1}. If we also use a normal random walk for proposing new values for~$A$ (as it's unconstrained) then a proposal~$A^*$ is accepted with probability $\min(1,C)$ where \begin{align*} C&=\frac{\pi(A^*|\lambda,\xvec)}{\pi(A|\lambda,\xvec)}\times \frac{q(A|A^*)}{q(A^*|A)} \\ &=\frac{\pi(A^*|\lambda,\xvec)}{\pi(A|\lambda,\xvec)}\qquad\text{since the proposal distribution is symmetric about zero}\\ &=\frac{e^{A^*}\,\pi_\alpha(e^{A^*}|\lambda,\xvec)}{e^A\,\pi_\alpha(e^A|\lambda,\xvec)}\\ &=\frac{\alpha^*\,\pi_\alpha(\alpha^*|\lambda,\xvec)}{\alpha\,\pi_\alpha(\alpha|\lambda,\xvec)}. \end{align*} Notice that the acceptance probabilities $B$ and $C$ are the same, that is, there is no (algorithmic) difference between using a log-normal random walk for a positive parameter or working on the log-scale and using a symmetric normal random walk. \end{enumerate} \section{Summary} \begin{enumerate}[(i)] \item Bayesian inference can be complicated when not using a conjugate prior distribution. \item One solution is to use Markov chain Monte Carlo (MCMC) methods. \item These work by producing realisations from the posterior distribution by constructing a Markov chain which has the posterior distribution as its stationary distribution. \item The MCMC methods we have studied are the Gibbs sampler, Metropolis within Gibbs algorithm and the Metropolis--Hastings algorithm. \item When obtaining output from these algorithms, we need to assess whether there needs to be a burn-in and whether the output needs to be thinned (by looking at traceplots and autocorrelation plots) using \texttt{mcmcAnalysis} and \texttt{mcmcProcess}. \item The (converged and thinned) MCMC output are realisations from the posterior distribution. It can be used to \begin{itemize} \item obtain the posterior distribution for any (joint) functions of the parameters (such as $\sigma=1/\sqrt\tau$ or $(\theta_1=\mu-\tau,\theta_2=e^{\mu+\tau/2})^T$); \item look at bivariate posterior distributions via scatter plots; \item look at univariate marginal posterior distributions via histograms or boxplots; \item obtain numerical summaries such as the mean, standard deviation and confidence intervals for single variables and correlations between variables. \end{itemize} \item Equi-tailed posterior confidence intervals can be determined from the MCMC output using \texttt{mcmcCi}. \end{enumerate} \section{Learning objectives} By the end of this chapter, you should be able to: \begin{itemize} \item explain why not using a conjugate prior generally causes problems in determining the posterior distribution \item describe the Gibbs sampler, explain why it is a Markov chain and give an outline as to why its stationary distribution is the posterior distribution \item describe the issues of processing MCMC output (burn-in, autocorrelation, thinning etc.) and interpret numerical/graphical output \item derive the full conditional densities for any posterior distribution and name these distributions if they are ``standard'' distributions given in the notes or on the exam paper \item describe a Metropolis-Hastings algorithm in general terms and when using either symmetric or non-symmetric random walk proposals or independence proposals \item describe the hybrid methods \textit{componentwise transitions} and \textit{Metropolis within Gibbs} \item provide a detailed description of \textbf{any} of the MCMC algorithms as they apply to generating realisations from \textbf{any} posterior distribution \end{itemize}