%09/03/11
\documentclass[titlepage]{article}
\usepackage{graphics}
\usepackage{natbib}
\evensidemargin 0.3cm
\oddsidemargin 0.3cm
\topmargin -1.5cm
\textwidth 15cm
\textheight 23.5cm
\columnsep 2cm

\newcommand{\E}{{\rm E}}
\newcommand{\var}{{\rm var}}
\newcommand{\covar}{{\rm covar}}
\newcommand{\eps}{\varepsilon}
\newcommand{\ux}{\underline{x}}
\newcommand{\ci}{\perp \! \perp}


\bibliographystyle{Chicago}
\begin{document}

\title{A Very Brief Introduction to MCMC}
\author{Malcolm Farrow \\
University of Newcastle upon Tyne}

\maketitle

 




\section{Monte Carlo Methods}
\subsection{Introduction}
Suppose we have a (vector) unknown (parameter) $ \theta, $ prior density $ 
p(\theta) $ and likelihood $ L(\theta; x). $ The main computational problem in 
Bayesian inference is usually the evaluation of an integral

\begin{equation}
\label{eq:int}
\int_{\theta} q(\theta) p(\theta) L(\theta;x).d \theta
\end{equation}

for some $ q(\theta). $

For example, if $ q(\theta) = 1 $ then (\ref{eq:int}) gives $ C^{-1} $ where $ 
C $ is the constant of proportionality for the posterior density. If $ 
q(\theta) $ is some function of $ \theta $ then the result of (\ref{eq:int}) 
multiplied by $ C $ gives the posterior expectation of that function.

We can evaluate the integral (\ref{eq:int}) approximately by choosing a set of 
values, $ \theta_{1}, \ldots , \theta_{n} $ of $ \theta, $ called {\em 
nodes},  and calculating

\begin{equation}
\label{eq:sum}
\sum_{i=1}^{n} w_{i} q(\theta_{i}) p(\theta_{i}) L(\theta_{i};x)
\end{equation}

where $ w_{1}, \ldots , w_{n} $ are {\em weights.} For example, suppose $ \theta 
$ is a scalar unknown such that $ 0< \theta <1. $ We might choose $ 
\theta_{i}=(i-0.5)/n $ and $ w_{i}=1/n. $ E.g., with $ n=1000, $ we would have $ 
\theta=0.0005, 0.0015, 0.0025, \ldots , 0.9995 $ and $ w_{i}=0.001. $ Here $ 
w_{i} $ is the step size. If we generalise this to a vector of parameters then $ 
w_{i} $ is the product of the step sizes. More generally, the points need not be 
equally spaced so the $ w_{i} $ need not be equal.

We might want to place more nodes in the important parts of the parameter space, 
that is where the function we are integrating is greatest. One way to do this is 
to generate nodes as random samples from a distribution with pdf $ 
g(\theta) $ and let $ w_{i} = [n g(\theta_{i})]^{-1}. $  Then

\[ \int_{\theta} q(\theta) p(\theta) L(\theta;x).d \theta = \int_{\theta} 
g(\theta) q(\theta) p(\theta) L(\theta;x)/g(\theta).d \theta \]

and we can think of this as the expected 
value of $ q(\theta) p(\theta) L(\theta;x)/g(\theta) $ over the distribution 
with pdf $ g(\theta). $ So we take random samples from the distribution with pdf 
$ g(\theta) $ and average the values of $ q(\theta) p(\theta) 
L(\theta;x)/g(\theta) $ obtained. Thus 
(\ref{eq:sum}) will tend to approximate (\ref{eq:int}) for large $ n. $ The 
variance of (\ref{eq:sum}) will be 

\[ n \, \var \left[ \frac{q(\theta) p(\theta) L(\theta;x)}{n g(\theta)} \right] 
= \frac{1}{n} \, \var \left[ \frac{q(\theta) p(\theta) 
L(\theta;x)}{g(\theta)}\right]. \]

This will be small if $ g(\theta) $ has a similar shape 
to $ |q(\theta) p(\theta) L(\theta;x)| $ in which case sampling in this way is 
called {\em importance 
sampling.}

\subsection{Example}
Suppose we take $ n+1 $ observations $ y_{0}, \ldots ,y_{n} $ on a stationary 
first order autoregressive process

\[ Y_{i} - \mu = \phi(Y_{i-1} - \mu) + \eps_{i} \]

where $ \eps \sim N(0, \sigma^{2}). $

The likelihood is

\begin{eqnarray*}
L(\theta;y) & = & (2 \pi)^{-(n+1)/2} \left( \frac{\sigma^{2}}{1- \phi^{2}} 
\right)^{-1/2} \exp\left[ - \frac{1}{2 \sigma^{2}/(1-\phi^{2})} ( y_{0} - 
\mu)^{2} \right] \\
            &   & \times ( \sigma^{2})^{-n/2} \exp\left[ - \frac{1}{2 
\sigma^{2}} \sum_{i=1}^{n} \{ y_{i} - \mu  - \phi(y_{i-1} - \mu) \}^{2} \right]
\end{eqnarray*}

For large $ n $ and reasonably vague prior the posterior density is dominated 
by 

\[ ( \sigma^{2})^{-n/2} \exp\left[ - \frac{1}{2 \sigma^{2}} \sum_{i=1}^{n} \{ 
y_{i} - \mu  - \phi(y_{i-1} - \mu) \}^{2} \right] . \]

We can approximate this shape by giving $ p = (\sigma^{2})^{-1} $ a gamma 
distribution, gamma$(1 + n/2, (1/2) \sum[y_{i} - \mu^{\star} - \phi^{\star} 
y_{i-1}]^{2}) $ and, independently, $ \phi $ and $ \tilde{\mu} = \mu(1-\phi) $ 
a normal distribution with mean $ \phi^{\star}, \mu^{\star}, $ the least 
squares estimates of $ \phi $ and $ \tilde{\mu} $ based on regressing $ y_{i} $ 
on $ y_{i-1}, $ and variance matrix given by that of the least squares 
estimates.

Now suppose, for example, we want to find the marginal posterior distribution 
of $ \phi. $ Then

\[ q(\theta) = I_{u} = \left\{ \begin{array}{cc} 1 & \phi < u \\ 0 & \phi \geq 
u \end{array} \right. \]

may be evaluated (for each of a range of values of $ u) $ along with $ w_{i}, 
p(\theta_{i}) $ and $ L(\theta_{i};y) $ and the sum (\ref{eq:sum}) formed by 
sampling $ p, \phi, \tilde{\mu}. $

\section{The Gibbs Sampler}
\subsection{Introduction}
Fortunately it is no longer necessary to consider apparently complicated 
approximation methods as in the previous section every time we want to make an 
inference. Use of the Gibbs sampler (and other {\em Markov chain Monte Carlo} 
methods) has been developed so that the computation for models of many types 
can now be done relatively straightforwardly. In particular software is 
available to do this and, in this course, we will be using the software package 
``BUGS'' (Bayes Using Gibbs Sampling).

Suppose we have two unknowns, $ X, Y. $ It may be difficult to write down their 
joint density function or integrate over it but suppose that for any value $ y 
$ we can sample from the conditional distribution of $ X|Y=y $ and vice versa. 
Suppose we want to find the expected value of $ q(X,Y). $ (This includes most, 
if not all, of the things we might want to find). Then the Gibbs sampling 
algorithm works as follows.

\begin{enumerate}
\item Choose starting values $ (x_{0}, y_{0}). $ (Let $ i=1). $

\item \label{step2} Take a sample $ Y_{i} $ from the conditional distribution 
of $ Y|X=x_{i-1}. $

\item Take a sample $ X_{i} $ from the conditional distribution of $ X|Y=y_{i}. 
$

\item Go to \ref{step2}. (Let $ i $ become $ i+1). $
\end{enumerate}

Actually the way a program like BUGS works involves a lot of details to do with 
how the model is specified and how samples are taken from the conditional 
distributions. For these you are referred to the BUGS manual (\cite{Spieg95}, available from the BUGS Web site), \cite{Gilks96} and the BUGS web site, at the following address.
\begin{verbatim}
  http://www.mrc-bsu.cam.ac.uk/bugs/ 
\end{verbatim}

``Eventually'' the sample $ (X_{i}, Y_{i}) $ are from the joint distribution of 
$ X,Y. $ Thus, having allowed for a suitable run-in period, we may accumulate 
values of $ q(X,Y) $ and find an average (or even plot a histogram).

In fact $ X,Y $ may each be vectors (i.e.\ sets of unknowns) and we may have 
more than two such sets and cycle around them.

For further comments see \cite{Shaw93}, \cite{Gilks96} etc. In 
particular there are important issues to do with 
convergence.

\subsection{Examples}
\subsubsection{Example 1 (very simple)}
This example is so simple that we would not really use Gibbs sampling in this 
case. However it illustrates how it works.

We have observations $ y_{1}, \ldots ,y_{n} $ from a normal distribution with 
mean $ \mu $ and variance $ \sigma^{2}=p^{-1}. $ The observations are 
independent given $ \mu $ and $ p. $ We have independent priors for $ \mu $ and 
$ p $ such that $ \mu \sim N(\mu_{0},\pi_{0}^{-1}) $ and $ p \sim {\rm 
gamma}(\alpha, \beta) $ (an ``inverse gamma'' prior for $ \sigma^{2}). $ Then 
the joint prior density is proportional to 

\[ p^{\alpha-1} \exp(-\beta p) \exp\{-\pi_{0}(\mu-\mu_{0})^{2}/2\}. \]

the likelihood is proportional to

\[ p^{n/2} \exp\left\{-p\sum(y_{i}-\mu)^{2}/2\right\}. \]

The posterior density is proportional to

\[ p^{\alpha-1+n/2} \exp\left\{-\beta p -\left[\pi_{0}(\mu-\mu_{0})^{2} + p 
\sum(y_{i}-\mu)^{2}\right]/2\right\} . \]

Suppose we fix a value for $ \mu, $ say $ \mu_{1}. $ Then the conditional 
posterior density for $ p $ is proportional to

\[ p^{\alpha-1+n/2} \exp\left\{-\left[\beta  + 
\sum(y_{i}-\mu)^{2}/2\right]p\right\} . \]

This is a gamma$(\alpha+n/2, \beta+ \sum[y_{i}-\mu]^{2}/2) $ distribution and we 
sample a value, say $ p_{1}, $ from this.

Next we fix $ p=p_{1}. $ Then the conditional posterior density for $ \mu $ is 
proportional to

\[ \exp\left\{ -\left[\pi_{0}(\mu-\mu_{0})^{2} + p_{1} 
\sum(y_{i}-\mu)^{2}\right]/2 \right\} . \]

That is

\[ \exp\left\{-\left[\mu^{2}(\pi_{0}+np_{1}) - 2 \mu\left( \pi_{0} \mu_{0} + 
p_{1} \sum y_{i}\right) + \pi_{0} \mu_{0}^{2} + p_{1} \sum 
y_{i}^{2}\right]/2\right\} . \]

This is proportional to

\[ \exp\left\{-(\pi_{0} + np_{1})\left[\mu - \left(\pi_{0} \mu_{0} + p_{1} \sum 
y_{i}\right)(\pi_{0} + np_{1})^{-1}\right]^{2}/2 \right\}. \]

This is a $ N([\pi_{0} \mu_{0} + p_{1} \sum y_{i}][\pi_{0} + n p_{1}]^{-1}, 
[\pi_{0} + n p_{1}]^{-1}) $ distribution and we sample a value, say $ \mu_{2}, $ 
from this.

We can then sample a new value for $ p, $ given $ \mu=\mu_{2} $ and so on.

\subsubsection{Example 2 (Taken from \cite{Gelf90})}
\label{subsubs:gelfand}

Suppose we have $ k $ samples of observations and observation $ j $ in sample $ 
i $ is

\[ Y_{ij} = \theta_{i} + \eps_{ij} \]

for $ j = 1, \ldots, J, $ where $ \eps_{ij} \sim N(0,\sigma^{2}_{\eps}) $ and $ 
\theta_{i} \sim N(\mu, \sigma^{2}_{\theta}) $ (all independent).

Suppose we have independent priors for the three parameters:

\begin{eqnarray*}
\mu & \sim & N(\mu_{0}, \sigma^{2}_{0}) \\
\sigma^{2}_{\theta} & \sim & {\rm IG}(a_{1}, b_{1}) \\
\sigma^{2}_{\eps} & \sim & {\rm IG}(a_{2}, b_{2})
\end{eqnarray*}

where IG stands for ``inverse gamma'' (i.e.\ $ (\sigma^{2})^{-1} $ has a gamma 
distribution).

Then

\begin{eqnarray*}
\sigma^{2}_{\theta} | Y, \mu, \theta, \sigma^{2}_{\eps} & \sim & {\rm IG}( 
a_{1} + k/2, b_{1} + (1/2) \sum[ \theta_{i} - \mu]^{2}) \\
\sigma^{2}_{\eps} | Y, \mu, \theta, \sigma^{2}_{\theta} & \sim & {\rm IG}( 
a_{2} + kJ/2, b_{2} + (1/2) \sum \sum [Y_{ij} - \theta_{i}]^{2}) \\
\mu | Y, \theta, \sigma^{2}_{\theta}, \sigma^{2}_{\eps} & \sim & N \left( 
\frac{ \sigma^{2}_{\theta} \mu_{0} + \sigma^{2}_{0} \sum \theta_{i}}{ 
\sigma^{2}_{\theta} + k \sigma^{2}_{0}} , \frac{\sigma^{2}_{\theta} 
\sigma^{2}_{0}}{\sigma^{2}_{\theta} + k \sigma^{2}_{0}} \right) \\
\theta | Y, \mu, \sigma^{2}_{\theta}, \sigma^{2}_{\eps} & \sim & N \left( 
\frac{J \sigma^{2}_{\theta}}{J \sigma^{2}_{\theta} + \sigma^{2}_{\eps}} \bar{Y} 
+ \frac{ \sigma^{2}_{\eps}}{J \sigma^{2}_{\theta} + \sigma^{2}_{\eps}} \mu 
\underline{1} , \frac{\sigma^{2}_{\theta} \sigma^{2}_{\eps}}{J 
\sigma^{2}_{\theta} + \sigma^{2}_{\eps}} I \right)
\end{eqnarray*}

where $ \underline{1} $ is a $ k \times 1 $ vector of 1's and $ I $ is a $ k 
\times k $ identity matrix.

\subsection{Convergence}
One of the main difficulties with the Gibbs sampler and other Markov chain Monte 
Carlo methods is convergence. We wish to obtain a sample of values from the 
posterior distribution. We know that the equilibrium distribution of the chain 
is the required posterior distribution but how do we know that the chain has 
reached equilibrium?

Of course we use a ``burn-in'' period when we take samples but discard them 
before we start collecting samples. In some straightforward models a burn-in of 
a few hundred samples may be reasonably assumed to be sufficient. In more 
complicated cases, involving posterior distributions with ``diagonal'' ridges or 
more than one local maximum, the chain can take a long time to move around the 
whole parameter space and a much longer burn-in may be required. We may also 
require a much larger sample to be sure that the chain has not just been 
temporarily ``stuck'' in one part of the parameter space.

Because of these difficulties much work has been done on {\em convergence 
diagnostics} and other output analysis. For further information, see, e.g., the 
CODA manual \citep{Best95}, Section 4.

\bibliography{examplerefs}




\end{document}

