\slidecolour
\mysection{Non-conjugate multi-parameter problems}

\frame{

\vspace{1cm}

\huge
\textbf{Chapter 4}

\vspace{1cm}

\Large
\textbf{Non-conjugate multi-parameter problems}

}

\slidecolour
\frame{
\frametitle{Consequences of strike action for MAS3902}
\begin{itemize}
\item Probably won't get all of Chapter 4 covered :-(  \pause
\item Will get enough done to make sure you're OK for the exam \pause
\item No time to start revision early -- will have three revision classes in teaching week 11 (w/b 8th May) \pause
\item Assignment 4:  Due \textcolor{blue}{\bf 4pm, Friday 5th May} \pause
\item Specimen paper already on Canvas: Solutions posted soon
\end{itemize}
}

\slidecolour
\begin{frame}
  \frametitle{Chapter \thesection. Non-conjugate multi-parameter
    problems}
\begin{block}{}
\begin{itemize}
\item Here we will study some multi-parameter problems in which the
  prior distribution does not have to be conjugate
\item Inferences are made by using techniques which simulate
  realisations from the posterior distribution
\item These methods are generally referred to as \alert{Markov Chain
    Monte Carlo} techniques, and often abbreviated to MCMC
\item Two main MCMC techniques:
\begin{enumerate}[1.]
\item the \alert{Gibbs sampler} -- at the forefront of the
  recent MCMC revolution
\item \alert{Metropolis-Hastings} sampling
\end{enumerate}
\item MCMC schemes based on the combination of these two fundamental
  techniques are still at the forefront of MCMC research
\end{itemize}
\end{block}
\end{frame}

\subsection{Why is inference not straightforward in non-conjugate problems?}
\begin{frame}
\frametitle{{\thesubsection} Why is inference not straightforward in non-conjugate problems?}
\vspace*{-3mm}
\begin{myexample}
\begin{itemize}
\item Consider again the problem in section~\ref{sec:nrs}
\label{ex:semiconjnormal}
\item[] $X_i|\mu,\tau\sim N(\mu,1/\tau)$, $i=1,2,\ldots,n$
  (independent)
\item Here we showed that a $NGa$ prior for
  $(\mu,\tau)^T$ was conjugate
\item But what if a $NGa(b,c,g,h)$ prior
  distribution does not adequately represent our prior beliefs?
\item Suppose instead that our prior beliefs are represented by
  independent priors for the parameters
{\small 
\[
\mu\sim N\left(b,\frac{1}{c}\right) 
\quad\quad\text{and}\quad\quad
\tau\sim Ga(g,h)
\]}
What is the posterior distribution for $(\mu,\tau)^T$?
\end{itemize}
\end{myexample}
\vspace*{-2mm}
\begin{soln}{}\end{soln}
\end{frame}

\begin{frame}
\begin{block}{In this case \ldots}
\begin{itemize}
\item What is the posterior mean of $\mu$ and of $\tau$? 
\item What are their marginal distributions?  
\item How can we calculate the moments $E(\mu^{m_1}\tau^{m_2}|\xvec)$
  of this posterior distribution?
\end{itemize}
\end{block}

\pause
\begin{block}{Marginal posterior density for $\mu$}
\small
\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}
\end{align*}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Marginal posterior density for $\tau$}
\small
\begin{align*}
\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*}
\end{block}

\pause
\begin{block}{General moments}
\scriptsize
\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*}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Comments}
\begin{itemize}
\item These integrals cannot be determined analytically
\item It is possible to use numerical integration methods or
  approximations (for large~$n$)
\item However, in general, the accuracy of
  the numerical approximation to the integral deteriorates as the
  dimension of the integral increases
\item How do we analyse models with a large number of parameters?
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Comment}
\begin{itemize}
\item Not using a conjugate prior distribution can cause many basic
  problems such as difficulty in plotting the marginal posterior density or
  determining posterior moments
\item But having to use conjugate priors is far too restrictive for
  many real data analyses: 
\begin{enumerate}[(i)]
\item our prior beliefs may not be captured using a conjugate prior
\item most models for complex data do not have conjugate priors
\end{enumerate}
\item Until relatively recently, 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.
% \item So how does using semi-conjugate priors help?
\end{itemize}
\end{block}
\end{frame}



\frame{
\frametitle{Examples}
\begin{itemize}
\item Poisson distribution for our data, $X_{i}|\lambda \sim Po(\lambda)$:
\begin{eqnarray*}
\lambda \sim Ga(g, h) &\longrightarrow &\pause \lambda|\boldsymbol{x} \sim \pause Ga(g+n\bar{x}, h+n) \pause \\
\lambda \sim \pause Weibull(a, b) &\longrightarrow & \lambda|\boldsymbol{x} \sim \textbf{??} \pause
\end{eqnarray*}
\item Normal distribution for our data, $X_{i}|\mu, \tau \sim N(\mu, 1/\tau)$:
\begin{eqnarray*}
\mu|\tau \sim N\left(b, \frac{1}{c\tau}\right), \tau\sim Ga(g,h) &\textcolor{white}{\longrightarrow} &\textcolor{white}{(\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H)} \pause \\ 
\longrightarrow (\mu, \tau)^{T} \sim NGa(b, c, g, h)		&\longrightarrow& (\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H) \\ \pause
\mu \sim N\left(b, \frac{1}{c}\right); \quad \tau\sim Ga(g,h) &\longrightarrow & (\mu,\tau)^{T}|\boldsymbol{x} \sim \textbf{??}
\end{eqnarray*}
\textcolor{white}{\bf All hell breaks loose.}
\end{itemize}
}

\frame{
\frametitle{Examples}
\begin{itemize}
\item Poisson distribution for our data, $X_{i}|\lambda \sim Po(\lambda)$:
\begin{eqnarray*}
\lambda \sim Ga(g, h) &\longrightarrow & \lambda|\boldsymbol{x} \sim  Ga(g+n\bar{x}, h+n)  \\
\textcolor{blue}{\lambda \sim  Weibull(a, b)} &\textcolor{blue}{\longrightarrow} & \textcolor{blue}{\lambda|\boldsymbol{x} \sim \textbf{??}}
\end{eqnarray*}
\item Normal distribution for our data, $X_{i}|\mu, \tau \sim N(\mu, 1/\tau)$:
\begin{eqnarray*}
\mu|\tau \sim N\left(b, \frac{1}{c\tau}\right), \tau\sim Ga(g,h) &\textcolor{white}{\longrightarrow} &\textcolor{white}{(\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H)}  \\ 
\longrightarrow (\mu, \tau)^{T} \sim NGa(b, c, g, h)		&\longrightarrow& (\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H) \\ 
\textcolor{blue}{\mu \sim N\left(b, \frac{1}{c}\right); \quad \tau\sim Ga(g,h)} &\textcolor{blue}{\longrightarrow} & \textcolor{blue}{(\mu,\tau)^{T}|\boldsymbol{x} \sim \textbf{??}}
\end{eqnarray*}
\textcolor{blue}{\bf All hell breaks loose.}
\end{itemize}
}



\slidecolour
\subsection{Simulation-based inference}
\begin{frame}
\frametitle{{\thesubsection} Simulation-based inference}
\vspace*{-2mm}
\begin{block}{}
\begin{itemize}
\item Can get around the problem of having to work out integrals
\item Base inferences on simulated realisations from the posterior
  distribution
\item This is the fundamental idea behind \textcolor{red}{\bf Markov Chain Monte Carlo} (MCMC) methods 
\item 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.
\end{itemize}
\end{block}

\begin{center}
{\large \textbf{Black box simulator demo in R}}

\textit{(Ignore the upcoming slides if you'd rather try out the R code!)}
\end{center}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item Imagine you wanted to know about the standard normal
  distribution -- its shape, its mean, its standard deviation 
\item But didn't know any mathematics so that you couldn't derive say
  the distribution's zero mean and unit variance 
\item However you've been given a ``black box'' which can simulate
  realisations from this distribution
\item Here we'll use the \texttt{R} function \texttt{rnorm()} as the
  black box simulator
\item If you decide to generate 1K realisations the output might look
  something like 
\end{itemize}
\begin{figure}[t!]
\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}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{In these plots}
\begin{itemize}
\item[Left:] the trace plot of the output -- the realisations
  from the black box sampler in the order they are produced
\item[Mid:] the autocorrelation (ACF) plot -- shows the correlation between $\{x_1,\ldots,x_{N-k}\}$ and $\{x_{k+1},\ldots,x_N\}$ for $k\geq0$. Call $k$ as the lag. It shows how correlated the realisations are at different lags. The lag~0 autocorrelation $corr(x_i,x_i)$ must be one (by definition). Here the (sample) correlation between say consecutive values $corr(x_i,x_{i+1})$ will be almost zero, which shows that the
simulator \texttt{rnorm()} produces independent realisations. This is also the case for
correlations at all positive lags. This sample ACF plot tells that the simulator generates independent sample.
\item[Right:] the density histogram of the realisations. This tells that what the standard normal distribution is like.
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item If you simulate another 1K realisations then the output you get
  will be a different set of realisations but from the same
  distribution
\item Both sets of realisations might look like
\end{itemize}
\begin{figure}[t!]
\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}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item We can also estimate various quantities of the standard normal
  distribution from these simulated datasets of 1K realisations

\bigskip
\begin{tabular}{lcrrcc}
& 1st Qu. &  Median &    Mean~ & 3rd Qu. & St.Dev.\\
\hline
Dataset 1 & -0.65240 & -0.00130 & -0.00192 & 0.64810 & 0.96049\\
Dataset 2 & -0.69880 & -0.09637 & -0.03274 & 0.67330 & 0.99599\\
\end{tabular}
\bigskip
\item The numerical summaries are slightly different but essentially
  tell the same story
\item In each Dataset, the mean and median are around zero and the
  standard deviation is around one.
\item We know from previous modules that there is sample variability
  in estimates of means from random samples
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item Now look at two datasets with 10K simulated realisations:

\bigskip
\begin{tabular}{lcrrcc}
& 1st Qu. &  Median &    Mean~ & 3rd Qu. & St.Dev.\\
\hline
Dataset 3 & -0.66820 & -0.00048 &  0.00370 & 0.68070 & 0.99593\\
Dataset 4 & -0.67130 &  0.01354 &  0.01008 & 0.67920 & 1.00691\\
\end{tabular}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.8\textwidth]{nrs2}
\caption{Summaries of the 10K realisations from the black box
  simulator}
%\label{fig:nrs1}
\end{center}
\end{figure}
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item Now estimates have much less sampling variability due to the
  larger sample size
\item In fact we can estimate any ``population'' quantity to any
  required accuracy simply by simulating a large enough collection of
  realisations
\bigskip
\item 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
\end{itemize}
\end{block}
\end{frame}

\slidecolour
\subsection{Motivation for MCMC methods}
\begin{frame}
\frametitle{{\thesubsection} Motivation for MCMC methods}
\vspace*{-2mm}
\begin{block}{}
\begin{itemize}
\item We consider a generic case where we want to simulate
  realisations of two random variables $X$ and $Y$ with joint density
  $f(x,y)$
\item 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)$
\item 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)$

\item Assume that simulating from $f(y|x)$ and $f(x|y)$ is
  straightforward
\end{itemize}
\end{block}
\end{frame}

\begin{frame}[fragile]
\frametitle{Example from MAS2906}
Suppose $X\sim U(0.1,0.5)$ and 
\begin{equation*}
f_{Y|X}(Y|X=x) = \begin{cases}
xe^{-xy}, &\text{when\ }y\geq 0,\\
0&\text{otherwise}.
\end{cases}
\end{equation*}
Write an R function to generate a $2\times n$ matrix of samples from this distribution, and produce a scatter plot. 

Note that
\begin{equation*}
Y|X=x \sim Exp(x).
\end{equation*}
\begin{verbatim}
 example.sim = function(n) {
     output = matrix(0, 2, n) 
     x = runif(n, 0.1, 0.5)
     y = rexp(n, rate = x)
     output[1, ] = x
     output[2, ] = y
     output
 }
\end{verbatim}
\end{frame}

\begin{frame}[fragile]
\frametitle{Another example from MAS2906}
Suppose $(X,Y)^{T} \sim N_{2}(\boldsymbol{\mu}, \Sigma)$;
\begin{itemize}
\item  $\boldsymbol{\mu} = (\mu_x, \mu_y)^{T}$ 	\\
\item  $\Sigma = \left(\begin{array}{c c} 
					\sigma_{x}^{2}	&\sigma_{xy}^{2} 	\\
								& \\
					\sigma_{xy}^{2}	&\sigma_{y}^{2} \end{array} \right)$
\end{itemize}\pause

\vspace{0.3cm}

What does the theory tell us about $X$, $Y$, $X|Y$ and $Y|X$? \pause
\begin{itemize}
\item [1.] $X \sim N(\mu_x, \sigma_{x}^{2}), Y \sim N(\mu_y, \sigma_{y}^{2})$ \pause
\item [2.] $X|Y=y \sim N\left(\mu_x + \rho \displaystyle \frac{\sigma_x}{\sigma_y}(y-\mu_{y}), \sigma_{x}^{2}(1-\rho^{2})\right)$ (same for $Y|X=x$)
\end{itemize}

\vspace{0.3cm}

{\color{gray}
\textit{Recall that the correlation, $\rho$, is}
\[
\rho = \frac{\sigma_{xy}^{2}}{\sigma_x \sigma_y} \longrightarrow \sigma_{xy}^{2} = \rho \sigma_x \sigma_y.
\]}
\end{frame}

\begin{frame}[fragile]
\frametitle{Another example from MAS2906}
\begin{itemize}
\item Suppose we want to investigate the marginal distributions of $X$ and $Y$ 
\item Suppose further that we did not know much theory about the bivariate Normal; we simply knew the results in (2) telling us about the conditional distributions
\item Could use a \textcolor{blue}{\bf Gibbs sampler}, in place of statistical theory, to investigate the marginals for $X$ and $Y$! \pause
\item Example: 
\[
\left( \begin{array}{c}
		X \\
		Y \end{array}\right) \sim N_{2}\left(\left[\begin{array}{c}
 								2	\\
								3	\end{array} \right]  ,\left[\begin{array}{c c}
												      1	&0.5	\\
													0.5	&1 \end{array} \right] \right).
\] \pause
\item
We see that $\mu_x = 2, \mu_y = 3; \sigma_x=1, \sigma_y=1, \sigma_{xy}^{2} = 0.5 \rightarrow \rho = 0.5$ \pause
\item From result 2, we can write down that
\begin{eqnarray*}
X|Y=y &\sim& N(2+0.5(y-3), 0.75) \quad \text{and}\\
Y|X=x &\sim& N(3+0.5(x-2), 0.75).
\end{eqnarray*}  \pause

\vspace{-0.5cm}

\item What about the marginals for $X$ and $Y$? 
\end{itemize}

\textcolor{gray}{\it Recall that we know from the theory that $X \sim N(2, 1)$ and $Y\sim N(3, 1)$.}
\end{frame}

\begin{frame}[fragile]
\frametitle{Another example from MAS2906}
\begin{center}
{\large \textbf{Bivariate Normal demo in R}}
\end{center}
\end{frame}

%%%%%% AFTER EASTER %%%%%%%%%%%

\frame{
\frametitle{Welcome back!}
\begin{itemize}
\item \textcolor{blue}{\bf This week:} Three sessions -- Monday (2), Tuesday (11), Friday (12).  Might use the Friday session as an extra lecture to make sure we finish early for revision \pause

\vspace{0.15cm}

\item \textcolor{blue}{\bf Next week:} Three sessions -- Monday (2), Tuesday (11), Thursday (2).  Although the Thursday class is a computer session, it will be really useful in helping you understand some of the methods discussed so far in Chapter 4, and should help you with exam prep for this material.  Assignment 4 due in Friday 4pm \pause

\vspace{0.15cm}

\item \textcolor{blue}{\bf Week after:} Three sessions -- Tuesday (11), Friday (12-2).  Revision
\end{itemize} 
}


\begin{frame}
\frametitle{Welcome back!}
What is Chapter 4 all about? \pause

\vspace{0.15cm}

\begin{itemize}
\item How to look at the posterior when we don't use a conjugate prior!

\vspace{0.15cm}

\item Why can't we just use the usual Bayes Rule?

\vspace{0.15cm}

\item What do we do?
\begin{itemize}

\item[--] Use \textcolor{blue}{\bf Simulation} to obtain a (random) sample of realisations from the posterior

\item [--] Simplest method: \textcolor{blue}{\bf Gibbs sampler}, using conditional distributions
\end{itemize}
\end{itemize}
\end{frame}



\frame{
\frametitle{Examples}
\begin{itemize}
\item Poisson distribution for our data, $X_{i}|\lambda \sim Po(\lambda)$:
\begin{eqnarray*}
\lambda \sim Ga(g, h) &\longrightarrow &\pause \lambda|\boldsymbol{x} \sim \pause Ga(g+n\bar{x}, h+n) \pause \\
\lambda \sim \pause Weibull(a, b) &\longrightarrow & \lambda|\boldsymbol{x} \sim \textbf{??} \pause
\end{eqnarray*}
\item Normal distribution for our data, $X_{i}|\mu, \tau \sim N(\mu, 1/\tau)$:
\begin{eqnarray*}
\mu|\tau \sim N\left(b, \frac{1}{c\tau}\right), \tau\sim Ga(g,h) &\textcolor{white}{\longrightarrow} &\textcolor{white}{(\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H)} \pause \\ 
\longrightarrow (\mu, \tau)^{T} \sim NGa(b, c, g, h)		&\longrightarrow& (\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H) \\ \pause
\mu \sim N\left(b, \frac{1}{c}\right); \quad \tau\sim Ga(g,h) &\longrightarrow & (\mu,\tau)^{T}|\boldsymbol{x} \sim \textbf{??}
\end{eqnarray*}
\textcolor{white}{\bf All hell breaks loose.}
\end{itemize}
}

\frame{
\frametitle{Examples}
\begin{itemize}
\item Poisson distribution for our data, $X_{i}|\lambda \sim Po(\lambda)$:
\begin{eqnarray*}
\lambda \sim Ga(g, h) &\longrightarrow & \lambda|\boldsymbol{x} \sim  Ga(g+n\bar{x}, h+n)  \\
\textcolor{blue}{\lambda \sim  Weibull(a, b)} &\textcolor{blue}{\longrightarrow} & \textcolor{blue}{\lambda|\boldsymbol{x} \sim \textbf{??}}
\end{eqnarray*}
\item Normal distribution for our data, $X_{i}|\mu, \tau \sim N(\mu, 1/\tau)$:
\begin{eqnarray*}
\mu|\tau \sim N\left(b, \frac{1}{c\tau}\right), \tau\sim Ga(g,h) &\textcolor{white}{\longrightarrow} &\textcolor{white}{(\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H)}  \\ 
\longrightarrow (\mu, \tau)^{T} \sim NGa(b, c, g, h)		&\longrightarrow& (\mu,\tau)^{T}|\boldsymbol{x} \sim NGa(B, C, G, H) \\ 
\textcolor{blue}{\mu \sim N\left(b, \frac{1}{c}\right); \quad \tau\sim Ga(g,h)} &\textcolor{blue}{\longrightarrow} & \textcolor{blue}{(\mu,\tau)^{T}|\boldsymbol{x} \sim \textbf{??}}
\end{eqnarray*}
\textcolor{blue}{\bf All hell breaks loose.}
\end{itemize}
}
\begin{frame}[fragile]
\frametitle{Example from MAS2906}
\begin{itemize}
\item Suppose we want to investigate the marginal distributions of $X$ and $Y$ 
\item Suppose further that we did not know much theory about the bivariate Normal; we simply knew results telling us about the conditional distributions
\item Could use a \textcolor{blue}{\bf Gibbs sampler}, in place of statistical theory, to investigate the marginals for $X$ and $Y$! \pause
\item Example: 
\[
\left( \begin{array}{c}
		X \\
		Y \end{array}\right) \sim N_{2}\left(\left[\begin{array}{c}
 								2	\\
								3	\end{array} \right]  ,\left[\begin{array}{c c}
												      1	&0.5	\\
													0.5	&1 \end{array} \right] \right).
\] \pause
\item
We see that $\mu_x = 2, \mu_y = 3; \sigma_x=1, \sigma_y=1, \sigma_{xy}^{2} = 0.5 \rightarrow \rho = 0.5$ \pause
\item From the theory in MAS2906, we can write down that
\begin{eqnarray*}
X|Y=y &\sim& N(2+0.5(y-3), 0.75) \quad \text{and}\\
Y|X=x &\sim& N(3+0.5(x-2), 0.75).
\end{eqnarray*}  \pause

\vspace{-0.5cm}

\item What about the marginals for $X$ and $Y$? 
\end{itemize}

\textcolor{gray}{\it Recall that we know from the theory that $X \sim N(2, 1)$ and $Y\sim N(3, 1)$.}
\end{frame}

\begin{frame}[fragile]
\frametitle{Example from MAS2906}
\begin{center}
{\large \textbf{Bivariate Normal demo in R}}
\end{center}
\end{frame}

\frame{
\frametitle{Key ingredients}
\begin{itemize}
\item To simulate bivariate pairs $(X,Y)$, we need to know the \textcolor{blue}{\bf conditional distributions} for $X|Y$ and $Y|X$ 
\item We then provide a \textcolor{blue}{\bf starting value} $X^{(0)}$ or $Y^{(0)}$ and then: \pause
\begin{itemize}
\item [--] Simulate $Y^{(1)}$ using the distribution for $Y|X=X^{(0)}$ \\
\textcolor{gray}{(or simulate $X^{(1)}$ using $X|Y = Y^{(0)}$)}
\item [--]Simulate $X^{(1)}$ using the distribution for $X|Y = Y^{(1)}$ \pause
\item [--]Simulate $Y^{(2)}$ using the distribution for $Y|X = X^{(1)}$
\item [--]Simulate $X^{(2)}$ using the distribution for $X|Y = Y^{(2)}$ \pause
\item [--]$\vdots$
\item [--]Simulate $Y^{(N)}$ using the distribution for $Y|X = X^{(N-1)}$
\item [--]Simulate $X^{(N)}$ using the distribution for $X|Y = X^{(N)}$
\end{itemize} \pause

\item ``Diagnostics'':
\begin{itemize}
\item [--] Trace plots -- to check the \textcolor{blue}{\bf burn-in period}
\item [--] ACF plots -- to check the \textcolor{blue}{\bf thinning interval}
\end{itemize} \pause

\item Now treat the thinned output as a random sample from the target distribution!
\end{itemize}
}




\slidecolour
\subsection{The Gibbs sampler}
\begin{frame}
\frametitle{{\thesubsection} The Gibbs sampler: Formal set-up}
\vspace*{-4mm}
\begin{block}{}
\begin{itemize}
\item 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
\]

\item 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)},\ldots,\theta_p^{(j-1)},\xvec) \\
\theta_2^{(j)} &\sim \pi(\theta_2|\theta_1^{(j)},\theta_3^{(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}
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Comments}
\begin{itemize}
\item 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$
\item 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
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{More comments}
\begin{itemize}
\item This is an iterative scheme
\item It needs a period to get to convergence: the \textit{burn-in}
  period
\item As it is a Markov chain, successive iterates are not independent
\item We can get accurate values for $E(\thetavec|\xvec)$,
  $SD(\thetavec|\xvec)$ or even $\pi(\thetavec|\xvec)$ by running the
  sampler for a long time
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Illustration of burn-in}
\begin{center}
\includegraphics[width=0.8\textwidth]{burnin}
\end{center}

\vspace*{-5mm}
\begin{itemize}
\item Different starting points 
\item But distributions of the realisations after iteration 500 are the ``same''
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Illustration of burn-in in two dimensions}
\begin{center}
\includegraphics[width=0.8\textwidth]{2-dim-Gibbs-trajectories.png}
\end{center}

\vspace*{-5mm}
\begin{itemize}
\item Different starting points 
\item Chains converge to the ``same'' distribution
\end{itemize}
\end{block}
\end{frame}



\begin{frame}
\begin{block}{{\thesubsection}.1 Processing output from a Gibbs sampler: Burn-in}
\begin{itemize}
\item Consider a $p=2$ parameter problem, with $\thetavec=(\mu,\tau)^T$ 
\item Run Gibbs sampler for $N$ iterations after burn-in, giving
\[
\{(\mu^{(1)},\tau^{(1)}),(\mu^{(2)},\tau^{(2)}),\ldots,(\mu^{(N)},\tau^{(N)})\}
\]

\item Can calculate features of the posterior distribution using their
  sample equivalents: $\bar\mu$, $\bar\tau$, $s_\mu^2$, $s_\tau^2$ and
  $r_{\mu\tau}$
\end{itemize}
\end{block}
\pause
\begin{block}{How accurate are these ``estimates''?}
\begin{itemize}
\item Difficult to determine exactly as they are not a \textit{random
    sample}
\item If it was a random sample then (for large $N$) approximate 95\%
  CI for $\mu$ is
\[
\bar\mu\pm \frac{1.96s_\mu}{\sqrt{N}}
\]
\item But we don't have a random sample! What to do?
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{{\thesubsection}.1 Processing output from a Gibbs sampler: autocorrelation}
\begin{itemize}
\item \textbf{Method one:} Assume the output follows some Markov chain model, and use CI of this model to approximate the true CI.
\item E.g. the autoregressive model with order $1$, AR(1), is 
\[
\mu^{(t)}-\mu=\theta (\mu^{(t-1)}-\mu)+\varepsilon_t,\ \mbox{where }\varepsilon_t\sim N(0,1).
\]
The approximate $95\%$ CI of $\mu$ depends on the sample autocorrelation $r(1)=Corr(\mu^{(t)},\mu^{(t+1)})$, and is 
\[
\bar\mu\pm \frac{1.96s_\mu}{\sqrt{N\{1-r(1)\}^2}}.
\]
Effective sample size: $N_{eff}=N\{1-r(1)\}^2$.
\item In general, MCMC output with positive autocorrelations has $N_{eff}<N$.
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{How to deal with autocorrelation}
\begin{itemize}
\item \textbf{Method two (more popular):} \texttt{thin} the output -- don't take every iterate
\item Sample autocorrelation function with lag~$k$ for~$\mu$ is
\[
r(k)=Corr(\mu^{(j)},\mu^{(j+k)}).
\]
\item For example, if $r(k)\simeq 0$, $k\geq 2$ then taking every 2nd
  iterate would give a posterior sample
\[
\{(\mu^{(1)},\tau^{(1)}),(\mu^{(3)},\tau^{(3)}),\ldots,(\mu^{(2j+1)},\tau^{(2j+1)},\ldots)\}
\]
with autocorrelation function $r^*(k)=0$, $k\geq 1$, that is, this
output has very low autocorrelation
\item Generally we look at the ACF plot and choose the level of
  thinning needed to give low autocorrelations
\item We will see how this works in some examples
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{$N(0,1)$ random sample}
\begin{center}
\includegraphics[width=\textwidth]{nrs-trace-acf}
\end{center}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Sample from a $MA(1)$ process}
\begin{align*}
x_t=\varepsilon_t+\theta\varepsilon_{t-1},\ \mbox{where }\varepsilon_t\sim N(0,1),\mbox{and }r(1)=\theta, r(k)=0, k\geq2.
\end{align*}
\begin{center}
\includegraphics[width=0.8\textwidth]{ma1-trace-acf}
\end{center}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Thinned sample from a $MA(1)$ process ($\texttt{thin}=2$)}
\begin{center}
\includegraphics[width=\textwidth]{ma1thin-trace-acf}
\end{center}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Sample from a $AR(1)$ process}
\begin{align*}
x_t=\theta x_{t-1}+\varepsilon_t,\ \mbox{where }\varepsilon_t\sim N(0,1),\mbox{and }r(k)\simeq r(1)^k, k\geq2.
\end{align*}
\begin{center}
\includegraphics[width=0.8\textwidth]{ar1-trace-acf}
\end{center}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Thinned sample from a $AR(1)$ process ($\texttt{thin}=10$)}
\begin{center}
\includegraphics[width=\textwidth]{ar1thin-trace-acf}
\end{center}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{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 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 posterior
  distribution
\end{enumerate}
\end{block}
\end{frame}


\begin{frame}
\frametitle{Recap from before Easter: Example 4.1}
\begin{itemize}
\item Consider again the problem in section~\ref{sec:nrs}
\label{ex:semiconjnormal}
\item[] $X_i|\mu,\tau\sim N(\mu,1/\tau)$, $i=1,2,\ldots,n$
  (independent)
\item Here we showed that a $NGa$ prior for
  $(\mu,\tau)^T$ was conjugate
\item But what if a $NGa(b,c,g,h)$ prior
  distribution does not adequately represent our prior beliefs?
\item Suppose instead that our prior beliefs are represented by
  independent priors for the parameters
{\small 
\[
\mu\sim N\left(b,\frac{1}{c}\right) 
\quad\quad\text{and}\quad\quad
\tau\sim Ga(g,h)
\]}
What is the posterior distribution for $(\mu,\tau)^T$?
\end{itemize}
\end{frame}
\begin{frame}
\frametitle{Recap from before Easter: Example 4.1}
\begin{soln}{
\begin{itemize}
\item No longer conjugate -- lots of horrible integrals to work out if we want to plot the posterior, or work out posterior summaries
\item Until the mid-1990s this made Bayesian inference in the non-conjugate case almost impossible
\end{itemize}
}
\end{soln}
\end{frame}

\begin{frame}
\frametitle{Recap from before Easter: Example 4.2}
\begin{myexample}
  Construct a Gibbs sampler for the posterior distribution in
  Example 4.1:
\bigskip
\begin{itemize}
\item Data: $X_i|\mu,\tau\sim N(\mu,1/\tau)$, $i=1,2,\ldots,n$ (independent)
\item Prior: 
\[
\mu\sim N\left(b,\frac{1}{c}\right) 
\quad\quad\text{and}\quad\quad
\tau\sim Ga(g,h),\qquad \text{independent}
\]
\end{itemize}
\end{myexample}

\begin{soln}{}\end{soln}
\end{frame}


\begin{frame}
\begin{block}{R functions in the nclbayes library}
\begin{itemize}
\item \R{gibbsNormal} -- implements this algorithm
\item \R{mcmcProcess} -- used to remove the \texttt{burn-in} and
  \texttt{thin} the output
\item \R{mcmcAnalysis} -- analyses the MCMC output
\end{itemize}
\end{block}

\pause
\begin{block}{Look at a particular numerical example}
\begin{itemize}
\item Data: $n=100$, $\bar x=15$ and $s=4.5$
\item Prior: $\mu\sim N(10,1/100)$ and $\tau\sim Ga(3,12)$,
  independently
\item Initialise Gibbs sampler at $(10,0.25)$
\end{itemize}
\end{block}


\begin{center}
{\large \textbf{Gibbs sampler demo in R}}

\textit{(Ignore the upcoming slides if you'd rather try out the R code!)}
\end{center}
\end{frame}

\begin{frame}[fragile]
\begin{block}{R commands used}
\color{red}{
\begin{verbatim}
library(nclbayes)
posterior=gibbsNormal(N=2000,initial=c(10,0.25),
    priorparam=c(10,1/100,3,12),n=100,xbar=15,s=4.5)
posterior2=mcmcProcess(input=posterior,burnin=1000,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 1000")
plot(posterior2,type="l",main="After deleting first 1000")
par(op)

mcmcAnalysis(posterior,rows=2,show=F)
mcmcAnalysis(posterior2,rows=2,show=F)
\end{verbatim}
}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Progress of the MCMC scheme}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=\textwidth]{gibbsNormal1}
\label{fig:gibbsNormal1}
\end{center}
\end{figure}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Analysis of MCMC output (all realisations)}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.85\textwidth]{gibbsNormal2}
\caption{Trace plots, autocorrelation plots and histograms of the
  Gibbs sampler output}
\label{fig:gibbsNormal2}
\end{center}
\end{figure}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Analysis of MCMC output (after deleting the first 10 iterations)}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.85\textwidth]{gibbsNormal3}
\caption{Trace plots, autocorrelation plots and histograms of the
  Gibbs sampler output}
\end{center}
\end{figure}
\end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{}
\begin{itemize}
\item Summaries of the output (after deleting the first 10 iterations)
\bigskip
\begin{verbatim}
N = 1000 iterations
       mu             tau         
 Min.   :13.54   Min.   :0.03017  
 1st Qu.:14.69   1st Qu.:0.04595  
 Median :15.01   Median :0.05089  
 Mean   :14.99   Mean   :0.05110  
 3rd Qu.:15.28   3rd Qu.:0.05562  
 Max.   :16.52   Max.   :0.07253  
Standard deviations:
         mu         tau 
0.444299473 0.007124057 
\end{verbatim}
\bigskip
\item Also $Corr(\mu,\tau|\xvec)=-0.02153$
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item Equi-tailed confidence intervals are calculated as follows:
\begin{itemize}
\item MCMC output has $N$ realisations $(\mu^{(j)},\tau^{(j)})$
\item Sort the $\mu^{(j)}$ and the $\tau^{(j)}$ into increasing order
\item CI end-points will be the $N\alpha/2$th and the
  $N(1-\alpha/2)$th values
\item Use \R{mcmcCi} command in the \R{nclbayes} package
\item Gives 95\% confidence intervals
\begin{align*}
\mu&:\quad  (14.123,~15.847)\\ 
\tau&:\quad (0.037949,~0.065571)
\end{align*}
\end{itemize}
\end{itemize}
\end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{Distributions of functions of the parameters}
\begin{itemize}
\item Now we have a sample from the posterior distribution, we can
  determine the posterior distribution for any function of the
  parameters
\item 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\left/\sqrt{\tau^{(j)}}\right.\]

\medskip
\item Summaries of the output for $\sigma$
\bigskip
\begin{verbatim}
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.713   4.240   4.433   4.457   4.665   5.758 
\end{verbatim}
\medskip
\item[] st.dev. 0.3178429
\medskip
\item[] Equi-tailed 95\% confidence interval for $\sigma$ is
$(3.9001,~5.1192)$.
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Joint and marginal distributions of $(\mu,\tau)^T$}
\begin{figure}[t!]
\vspace*{-3mm}
\begin{center}
\includegraphics[width=0.7\textwidth]{gibbsNormal3a}
\end{center}
\label{fig:gibbsNormal3}
\end{figure}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Joint and marginal distributions of $(\mu,\sigma)^T$}
\begin{figure}[t!]
\vspace*{-3mm}
\begin{center}
\includegraphics[width=0.7\textwidth]{gibbsNormal3b}
\end{center}
\end{figure}
\end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{R commands used}
\color{red}{
\begin{verbatim}
library(nclbayes)
posterior=gibbsNormal(N=2000,initial=c(10,0.25),
    priorparam=c(10,1/100,3,12),n=100,xbar=15,s=4.5)
posterior2=mcmcProcess(input=posterior,burnin=1000,thin=1)

mcmcAnalysis(posterior,rows=2,show=F)
mcmcAnalysis(posterior2,rows=2,show=F)

cor(posterior2)
mcmcCi(posterior2,level=0.95)

sigma=1/sqrt(posterior2[,2])
summary(sigma)
sd(sigma)
\end{verbatim}
}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{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}
\end{block}
\end{frame}

\begin{frame}
\begin{myexample}
\begin{itemize}
\item Gibbs sampling can also be used when using a conjugate prior
\item Construct a Gibbs sampler for the problem in
  Example~3.\ref{ex:cavbiv}: Cavendish's data on the earth's density
\item Data: random sample from a normal distribution with unknown
  mean~$\mu$ and precision~$\tau$, that is
\[
X_i|\mu,\tau\sim N(\mu,1/\tau),\quad i=1,2,\ldots,n\quad\text{
  (independent)}
\]
\item Prior: $NGa$ distribution for $(\mu,\tau)^T$
\end{itemize}
\end{myexample}
\begin{soln}{}\end{soln}
\end{frame}

\begin{frame}
\begin{block}{R function in the nclbayes library}
\begin{itemize}
\item \R{gibbsNormal2} -- implements this algorithm
\end{itemize}
\end{block}

\pause
\begin{block}{Reanalysis of Cavendish's data}
\begin{itemize}
\item Data: $n=23$, $\bar x=5.4848$, $s=0.1882$
\item Prior: $NGa(b=5.41,c=0.25,g=2.5,h=0.1)$ 
\item Have seen in Example~3.\ref{ex:cavbiv} that
\begin{itemize}
\item posterior is $NGa(B=5.4840,C=23.25,G=14,H=0.5080)$
\item marginals are
\vspace*{-2mm}
\begin{align*}
\mu|\xvec&\sim t_{2G=28}(B=5.4840,H/(GC)=0.001561)\\
\noalign{\smallskip}
\tau|\xvec&\sim Ga(G=14,H=0.5080)
\end{align*}
\end{itemize}
\medskip
\item We now analyse this problem using the Gibbs sampler
\item Can check that the Gibbs sampler is producing realisations from
  the correct distribution by comparing its output with the above
  theoretical marginal distributions
\end{itemize}
\end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{Main R commands used}
\color{red}{
\begin{verbatim}
library(nclbayes)
posterior=gibbsNormal2(N=2000,initial=c(5.41,25),
    priorparam=c(5.41,0.25,2.5,0.1),n=23,xbar=5.4848,s=0.1882)
posterior2=mcmcProcess(input=posterior,burnin=1000,thin=1)

mcmcAnalysis(posterior,rows=2,show=F)
mcmcAnalysis(posterior2,rows=2,show=F)
\end{verbatim}
}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Analysis of MCMC output (after deleting the first 1000 iterations)}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.85\textwidth]{gibbsNormal4}
\caption{Trace plots, autocorrelation plots and histograms of the
  Gibbs sampler output}
\end{center}
\label{fig:gibbsNormal4}
\end{figure}
\end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{}
\begin{itemize}
\item Summaries of the output (after deleting the first 1000 iterations)
\medskip
\begin{verbatim}
N = 1000 iterations
       mu             tau       
 Min.   :5.376   Min.   :11.01  
 1st Qu.:5.456   1st Qu.:22.09  
 Median :5.483   Median :27.06  
 Mean   :5.484   Mean   :27.66  
 3rd Qu.:5.512   3rd Qu.:32.60  
 Max.   :5.605   Max.   :56.16  
Standard deviations:
        mu        tau 
0.04090114 7.47836536 
\end{verbatim}
\item Check:
\vspace*{-4mm}
\begin{align*}
E(\mu|\xvec)&=B=5.4840,& SD(\mu|\xvec)&=\sqrt{\frac{H}{(G-1)C}}=0.04100\\
E(\tau|\xvec)&=\frac{G}{H}=27.559,& SD(\tau|\xvec)&=\frac{\sqrt{G}}{H}=7.3655
\end{align*}
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
  \begin{block}{Marginal posterior distributions of $\mu$ and
      $\sigma$}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.48\textwidth]{gibbsNormal5}
\includegraphics[width=0.48\textwidth]{gibbsNormal6}
%\caption{Histograms of the Gibbs sampler output and the (known) marginal densities}
\label{fig:gibbsNormal5}
\end{center}
\vspace*{-3mm}
\end{figure}
\end{block}

\pause
\begin{block}{Conclusion}
\begin{itemize}
\item Posterior means match pretty closely
\item Posterior standard deviations match pretty closely
\item Marginal posterior denisties match pretty closely
 \pause
 \item They are close enough within sampling error -- but could run
  sampler for more iterations if we wanted more accurate
  estimates \pause
\item This confirms that the Gibbs sampler does indeed produce
  realisations from the correct posterior distribution
\end{itemize}
\end{block}
\end{frame}
 
\begin{frame}
\begin{block}{Accuracy of posterior estimates}
\begin{itemize}
\item Posterior output: sample mean~$\bar\mu$, standard
deviation~$s_\mu$ 
\item How accuracte are these estimates of $M=E(\mu|\xvec)$ and $\Sigma=SD(\mu|\xvec)$?

\bigskip
\item Each time we run an MCMC scheme, we obtain a different sample from the
posterior distribution
\item Suppose we have a large sample with effective sample size $N$
\item It's easiest to think of this as being a sample which
has been thinned so that it has $N$ values, say $\mu_1,\ldots,\mu_N$
(and effective sample size~$N$)
\item To quantify accuracy, we need to make an assumption about the posterior 
distribution
\item If the data sample size~$n$ is large then the posterior distribution will be
approximately normal
\item Think of the MCMC output as being a random sample from a normal distribution
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item When the random sample $\theta_1,\ldots,\theta_N$ is from the normal $N(a,1/d)$, from example 3.1, the asymptotic posterior distribution
about the mean and precision $(a,d)^T$is
\[
a|\boldsymbol{\theta}\sim N(\bar \theta,s^2/N), \qquad d|\boldsymbol{\theta}\sim
N\{1/s^2,2/(Ns^4)\},\qquad\text{independently}
\]
\item If the Gibbs outputs $\mu_1,\ldots,\mu_N$ come from a normal distribution $N(M,\Sigma^2)$, rewriting this result in terms of the MCMC sample mean~$\bar\mu$,
standard deviation~$s_\mu$ and the parameters they estimate gives
posterior distributions
\[
M\sim N(\bar\mu,s_\mu^2/N), \qquad 
\Sigma^{-2}\sim N\{1/s_\mu^2,2/(Ns_\mu^4)\},\qquad\text{independently}
\]
\item Therefore an approximate 95\% HDI for~$M$ is
\[
\bar\mu\pm z_{0.025}\frac{s_\mu}{\sqrt{N}}
\simeq\bar\mu\pm\frac{2s_\mu}{\sqrt{N}}
\]
since $z_{0.025}\simeq 2$.
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\vspace*{-3mm}
\begin{block}{}
\begin{itemize}
\item Also, from the posterior distribution for $\Sigma^{-2}$, we have
\begin{align*}
&P\left(\frac{1}{s_\mu^2}-2\sqrt{\frac{2}{Ns_\mu^4}}<\Sigma^{-2}<
\frac{1}{s_\mu^2}+2\sqrt{\frac{2}{Ns_\mu^4}}\right)\simeq 0.95\\
\Longrightarrow\quad
&P\left(\frac{1-2\sqrt{2/N}}{s_\mu^2}<\Sigma^{-2}
<\frac{1+2\sqrt{2/N}}{s_\mu^2}\right)\simeq 0.95\\
\Longrightarrow\quad
&P\left(\frac{s_\mu}{\sqrt{1+2\sqrt{2/N}}}<\Sigma<
\frac{s_\mu}{\sqrt{1-2\sqrt{2/N}}}\right)\simeq 0.95.
\end{align*}
Therefore a 95\% confidence interval for $\Sigma$ is
\[
s_\mu\left(1\pm 2\sqrt{2/N}\right)^{-1/2}
\simeq s_\mu\left(1\pm\frac{1}{2}\times 2\sqrt{2/N}\right)
=s_\mu\pm s_\mu\sqrt{\frac{2}{N}}
\]
\item It can be shown that these accuracy calculations are fairly accurate
even when the posterior distribution (from which we have the MCMC
sample) is not particularly normal
\end{itemize}
\end{block}
\end{frame}


\begin{frame}
\begin{myexample}
\label{ex:gammars}
\begin{itemize}
\item Data: 
$X_i|\alpha,\lambda\sim Ga(\alpha,\lambda)$, $i=1,2,\ldots,n$
  (independent), where the index~$\alpha$ and
  scale parameter~$\lambda$ are unknown
\item Prior: $\alpha\sim Ga(a,b)$ and $\lambda\sim Ga(c,d)$,
  independent 
\medskip
\item Determine the posterior
density for $(\alpha,\lambda)^T$ and hence the posterior conditional
densities for $\alpha|\lambda$ and $\lambda|\alpha$
\end{itemize}
\end{myexample}
\begin{soln}{}\end{soln}
\end{frame}
  
\begin{frame}
\begin{block}{Problem}
\begin{itemize}
\item Although the FCD for $\lambda$ is a standard distribution and
  easy to simulate from, the FCD for $\alpha$ is NOT!
\item Therefore we can't use a Gibbs sampler for this analysis
\item Fortunately there are other methods we can use \ldots
\end{itemize}
\end{block}
\end{frame}

%\end{document}

\slidecolour
\subsection{Metropolis--Hastings sampling}
\begin{frame}
\frametitle{{\thesubsection} Metropolis--Hastings sampling}
\begin{block}{}
\begin{itemize}
\item The Gibbs sampler is a very powerful tool
\item Only useful if the full conditional distributions (FCDs) are
  standard distributions (which are easy to simulate from)
\item Fortunately there is a class of methods which can be used when
  the FCDs are non-standard
\item These methods are known as Metropolis-Hastings schemes
\end{itemize}
\end{block}
\begin{block}{}
\begin{itemize}
\item Want to simulate realisations from the posterior density
$\pi(\thetavec|\xvec)$ 
\item All of the FCDs are non-standard 
\item Choose a \alert{proposal distribution} with density
  $q(\thetavec^*|\thetavec)$, which is easy to simulate from
\item This distribution gives us a way of proposing new
  values~$\thetavec^*$ from the current value~$\thetavec$
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Metropolis--Hastings algorithm}
\begin{enumerate}
\item Initialise the iteration counter to $j=1$, and initialise the
chain to $\thetavec^{(0)}$
\item Generate a \alert{proposed} value $\thetavec^*$ using the
  proposal distribution $q(\thetavec^*|\thetavec^{(j-1)})$
\pause
\item Evaluate the \alert{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; that is:
\begin{center}
\begin{tabular}{|c|c c|} \hline
$\thetavec^{(j)}$	&$\thetavec^*$	&$\thetavec^{(j-1)}$ \\ \hline
Probability		&$\alpha$		&$1-\alpha$		\\ \hline
\end{tabular}
\end{center}
\pause
\item Change the counter from $j$ to $j+1$ and return to step 2
\end{enumerate}
\end{block}
\end{frame}




\begin{frame}
\begin{block}{Comments}
\begin{itemize}
\item At each stage, a new value is generated from the proposal
  distribution
\item This value is either accepted, in which case the chain moves, or
  rejected, in which case the chain stays where it is
\item 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 
\item The posterior density $\pi(\thetavec|\xvec)$ only enters into the
  acceptance probability as a ratio, and so the method can be used
  when the posterior density is only 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\}
\end{align*}
\item This algorithm defines a Markov chain with
  $\pi(\thetavec|\xvec)$ as its stationary distribution
\item It holds for all possible proposal distributions (subject to
  them generating realisations from the full parameter space)
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Possible proposal distributions}
\begin{itemize}
\item Are some choices of proposal distribution better than others?
\item We now discuss some commonly used proposal distributions
\end{itemize}
\end{block}

\pause
\begin{block}{{\thesubsection}.1 Symmetric chains (Metropolis method)}
\begin{itemize}
\item Use a symmetric proposal distribution:
  $q(\thetavec^*|\thetavec)=q(\thetavec|\thetavec^*)$,
  $\forall~\thetavec,\thetavec^*$
\item The acceptance probability simplifies to
\[
\alpha(\thetavec,\thetavec^*) = 
\min\left\{ 1, \frac{\pi(\thetavec^*|\xvec)}{\pi(\thetavec|\xvec)} \right\} = \min\left\{1, \frac{\pi(\thetavec^*)\,f(\xvec|\thetavec^*)}{\pi(\thetavec)\,f(\xvec|\thetavec)}  \right\}
\]
and hence does not involve the proposal density at all 
\item Proposed moves which will take the chain to a region of higher
  posterior density are always accepted
\item Moves which take the chain to a region of lower posterior
  density are accepted with probability proportional to the ratio of
  the two densities
\item Moves which will take the chain to a region of very low density
  will be accepted with very low probability
\end{itemize}
\end{block}
\end{frame}


\begin{frame}
\begin{block}{Random walk proposals}
\begin{itemize}
\item Consider the random walk proposal in which the proposed value
$\thetavec^*$ at stage $j$ is
\[
\thetavec^*=\thetavec^{(j-1)}+\vec{w}_j
\]
where the $\vec{w}_j$ are independent and identically distributed
random $p\times 1$ vectors (completely independent of the state of the
chain)
\item Suppose that the $\vec{w}_j$ have density $f(\cdot)$, which is
  easy to simulate from, has mean~$\vec{0}$ and is symmetric about its
  mean
\item We can then simulate an \alert{innovation} $\vec{w}_j$, and set
  the \alert{proposal} value to
  $\thetavec^*=\thetavec^{(j-1)}+\vec{w}_j$
\item What distribution should we use for~$f(\cdot)$? 
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Choice of innovation distribution}
\small
\begin{itemize}
\item A distribution which is simple and easy to simulate from is
  always a good idea; for example, the uniform or normal\quad$\cdots\cdots$ \\ Normal is
  generally better, but is a bit more expensive to simulate
\item What variance should we choose for the innovation distribution?
\item This choice will affect the acceptance probability, and hence
  the overall proportion of accepted moves
\item 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 ``\alert{cold}''
\item 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
  ``\alert{hot}''
\item 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
\item Thus, the variance of the innovation distribution should be
  ``tuned'' to get an acceptance rate of around this level
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\frametitle{Example}
\begin{itemize}
\item Asked by the UK Met Office to analyse annual maximum wind speed observations (mph) at 206 locations across the UK \pause
\item Suppose $X_{i}, i=1, ..., 50$ represent these annual maximum wind observations for a particular location in the Peak District \pause
\item An appropriate model is the \textit{Generalised Extreme Value} distribution, that is, $X_{i}|\mu, \sigma, \xi \sim GEV(\mu, \sigma, \xi)$ \pause
\item After discussions with Simon Brown, a climate scientist at the Met Office, we use: 
\begin{eqnarray*}
\mu &\sim & N(55, 2) \\
\sigma &\sim& Ga(20, 10) \quad \text{and}	\\
\xi &\sim & N(0, 0.5) 
\end{eqnarray*}
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Example}
\begin{itemize}
\item No conjugate prior distributions exist for this model, and the full conditionals are not of standard form \pause
\item We need to use Metropolis-Hastings; we use Normal random walk updates:
\begin{eqnarray*}
\mu^{*}	&=& \mu^{(j-1)} + w_{\mu_{j}}, \\
\log \sigma^{*}	&=& \log \sigma^{(j-1)} + w_{\sigma_{j}} \quad \text{and}	\\
\xi^{*}	&=& \xi^{(j-1)}+w_{\xi_{j}},
\end{eqnarray*}
where 
\begin{eqnarray*}
w_{\mu_{j}} &\sim&N(0, \epsilon_{\mu}),\\
w_{\sigma_{j}} &\sim & N(0, \epsilon_{\sigma}) \quad \text{and}	\\
w_{\xi}	&\sim& N(0, \epsilon_{\xi})
\end{eqnarray*}
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Example}
\begin{center}
\includegraphics[width=8cm]{../../EXAM_PAPERS/2020/mcmcplots.pdf}
\end{center}
\end{frame}


\slidecolour
\subsection{Hybrid methods}
\begin{frame}
\frametitle{{\thesubsection} Hybrid methods}
\begin{block}{}
\begin{itemize}
\item The Metropolis-Hastings method in last section samples $\thetavec$ in a single block, i.e. all parameters are updated at the same time.
\item When the number of parameter is large, it is difficult to choose the proposal distribution. For example, when there are $5$ parameters, the proposal distribution $N(0,\Sigma)$ has $5\times5=25$ tuning parameters in $\Sigma$.
\item Recall the Gibbs sampler for how the componentwise update works.
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{{\thesubsection}.1 Componentwise transitions}
\begin{itemize}
\item Given a posterior distribution with FCDs 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
\item This is like the Gibbs sampler, but each component update is a
  Metropolis-Hastings update, rather than a direct simulation from the
  FCD
\item Each of these steps will require its own proposal distribution
\end{itemize}
\end{block}

\end{frame}

\begin{frame}
\begin{block}{The 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 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\\
The distributions $\pi(\theta_i|\thetavec_{-i}^{(j)},\xvec)$ are
just the FCDs
\end{block}
\end{frame}



\begin{frame}
\frametitle{Hybrid methods}
\begin{block}{}
\begin{itemize}
\item The Gibbs sampler can be used to sample from multivariate posterior distributions provided that we can simulate from the full conditional distributions (FCDs). The acceptance probability is $1$.
\item The Metropolis-Hastings method can be used to sample from awkward FCDs. The acceptance probability is less than $1$.
\item We can combine these, when some FCDs can be simulated from directly and some can not, to increase the acceptance probability.
\end{itemize}
\end{block}
\end{frame}


\begin{frame}
\begin{block}{{\thesubsection}.2 Metropolis within Gibbs}
\begin{itemize}
\item Given a posterior distribution with full conditional distributions,
\begin{itemize}
\item some of which may be simulated from directly, 
\item and others of which have Metropolis-Hastings updating schemes,
\end{itemize}
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
\item 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 
\item 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)$ 
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item 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,
\end{align*}
that is, we always accept the proposal from the FCD
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{myexample}
\begin{itemize}
\item Construct an MCMC scheme for the problem in
  Example~{\thesection}.\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
\item 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{itemize}
\end{myexample}
\begin{soln}{}\end{soln}
\end{frame}

\begin{frame}[fragile]
\begin{block}{R code}
\begin{itemize}
\item Use the function \R{mwgGamma} in the \R{nclbayes} library
\item Data: $n=50$, $\bar x=0.62$, $\bar x_g=0.46$ and $s=0.4$
\item Prior: $a=2$, $b=1$, $c=3$, $d=1$
\item Use a normal random walk proposal with variance $\Sigma_\alpha=0.9^2$
\item This gives a reasonable acceptance probability of 0.237
\end{itemize}

\color{red}{
\begin{verbatim}
library(nclbayes)
posterior=mwgGamma(N=20000,initial=(xbar/s)^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}
}

\begin{itemize}
\item Using $\texttt{burnin}=10$, $\texttt{thin}=20$ produces (almost)
  un-autocorrelated posterior output
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Analysis of MCMC output (all realisations)}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.85\textwidth]{metropwithingibbs}
\caption{Trace plots, autocorrelation plots and histograms of the Metropolis within Gibbs sampler output}
\end{center}
\end{figure}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{Analysis of MCMC output (with $\texttt{burnin}=10$, $\texttt{thin}=20$)}
\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.85\textwidth]{metropwithingibbs2}
\caption{Trace plots, autocorrelation plots and histograms of the Metropolis within Gibbs sampler output}
\end{center}
\end{figure}
\end{block}
\end{frame}

\slidecolour
\subsection{Summary}
\begin{frame}
\frametitle{{\thesubsection} Summary}
\begin{block}{}
\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 \R{mcmcAnalysis} and \R{mcmcProcess}
\end{enumerate}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{enumerate}[(i)]
\item[(vi)] 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[(vii)] Equi-tailed posterior confidence intervals can be
  determined from the MCMC output using \R{mcmcCi}
\end{enumerate}
\end{block}
\end{frame}

\subsection{Learning objectives}
\slidecolour
\begin{frame}
\frametitle{{\thesubsection} Learning objectives}
\vspace*{-3mm}
\begin{block}{}
By the end of this chapter, you should be able to:
\begin{itemize}
\item[1.] explain why not using a conjugate prior generally
  causes problems in determining the posterior distribution
\item[2.] 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[3.] describe the issues of processing MCMC output
  (burn-in, autocorrelation, thinning etc.) and interpret
  numerical/graphical output
\item[4.] 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
\end{itemize}
\end{block}
\end{frame}

\begin{frame}
\begin{block}{}
\begin{itemize}
\item[5.] describe a Metropolis-Hastings algorithm in general
  terms and when using symmetric random walk
  proposals or independence proposals or general proposals
\item[6.] describe the hybrid methods: \textit{componentwise
    transitions} and \textit{Metropolis within Gibbs}
\item[7.] provide a detailed description of \textbf{any} of the MCMC
  algorithms as they apply to generating realisations
  from \textbf{any} posterior distribution
\end{itemize}
\end{block}
\end{frame}


\begin{comment}
\begin{frame}
\begin{block}{}
~\vspace*{3cm}

\begin{center}
\textbf{\LARGE THE END}
\end{center}
\vspace*{3cm}~
\end{block}
\end{frame}
\end{comment}