\documentstyle[12pt,fullpage]{article}

\frenchspacing

\long\def\comment#1{}

\def\vfill{\vspace*{\fill}}
\def\gap{\hspace*{0.5em}}
\def\vgap{\vspace*{0.5em}}

\let\vfill\relax

\def\bul{$\bullet$ }
\let\ul\underline
\def\beq{$$}
\def\eeq{$$}
\def\bea{\begin{eqnarray*}}
\def\eea{\end{eqnarray*}}
\def\implies{\Rightarrow}
\def\tendsto{\rightarrow}

\let\alias\newcommand

\def\Cov{\mbox{Cov}\,}
\def\Var{\mbox{Var}\,}
\alias{\x}[1]{x_{{#1}1}, \ldots, x_{{#1}J}}
\alias{\y}[1]{Y_{{#1}1}, \ldots, Y_{{#1}J_1}}
\alias{\z}[1]{Z_{{#1}1}, \ldots, Z_{{#1}J_2}}
\alias{\X}[1]{X_{{#1}1}, \ldots, X_{{#1}J}}
\alias{\th}{\theta}
\alias{\Th}{\Theta}
\alias{\w}{\omega}
\alias{\ua}{{\ul\alpha}}
\alias{\combos}[2]{\renewcommand{\arraystretch}{.75}%
        \left(\begin{array}{@{}c@{}}#1\\#2\end{array}\right)%
                \renewcommand{\arraystretch}{1}}
\def\sumsum{\mathop{\sum\sum}}

\def\bs{\beta\mbox{'s}}
\def\ts{\theta\mbox{'s}}

\def\th{\theta}
\def\bt{\beta}
\def\tone{\theta^{1}}
\def\bone{\beta^{1}}
\def\lone{\lambda^1}
\def\tzero{\theta^{0}}
\def\bzero{\beta^{0}}
\def\lzero{\lambda^0}
\def\iprime{{i^{'}}}

\def\nu{k}

\def\M{{\cal M}}

\def\rest{\mbox{\em rest}}

\parskip=0.6em
\parindent=0pt



\begin{document}

\begin{center}

{\large\bf Some notes on Markov Chain Monte Carlo: Metropolis-Hastings
  and Gibbs Steps for Exploring Posterior Distributions
}

\end{center}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Practical Bayesian Inference}

In the Bayesian paradigm we build a {\em joint model} for 
{\em data} $\ul y = (y_1, \ldots, y_N)$, and 
{\em parameters} $\tau$:

\[
\begin{array}[t]{rcl}
        p(\ul y,\tau) & = & p(\ul y|\tau)\cdot p(\tau)
        \end{array}
\]
and we always want to know something about the {\em posterior density}
\[
\begin{array}[t]{rcl}
        p(\tau|\ul y) & = & \displaystyle\frac{p(\ul y|\tau)\cdot p(\tau)}
        {\displaystyle\int p(\ul y|t)\cdot p(t)\,\, dt} \\
        & \propto & p(\ul y|\tau)\cdot p(\tau) ~~ = ~~ p(\ul y,\tau)
        \end{array}
\]
For example, 
\begin{itemize}
\item {\em Posterior mode}: $argmax_{\tau} \, p(\tau|\ul y)$
\item {\em Posterior mean}: $E[\tau|\ul y] = \int \tau \, p(\tau|\ul y) \,
  d\tau$ (or some other risk-minimizer)
\item {\em Credible interval}: Find a set $A$ such that $p(\tau\in A |
  \ul y) = 0.95$ 
\item {\em Graph or ``shape''} of $p(\tau|\ul y)$
\end{itemize}

{\ul{All} of these operations
% \begin{itemize}
% \item Point Maximization
% \item Integration
% \item Location of modal regions
% \item Graphing
% \end{itemize}
  require that you know a lot about the shape of $p(\tau|\ul y)$.  {\em
    How do you learn enough about the shape of $p(\tau|\ul y)$?}

\begin{itemize}

\item If $p(\tau)$ is {\em conjugate} to $p(\ul y|\tau)$, then $p(\tau|\ul y)$
will belong to the same family as $p(\tau)$.

\item If not, then you have a numerical integration problem.  These
problems are difficult because they depend on the ``interesting''
features of $p(\tau|\ul y)$.  
\begin{itemize}
\item Difficult but do-able if $p(\tau|\ul y)$ is {\em smooth} and {\em
low-dimensional}.
\item Even if $p(\tau|\ul y)$ is smooth,  as
  $d$ $=$ dim($\tau$) grows, space becomes too large for deterministic
  methods. {\em Curse of dimensionality:} For $n$ ``evaluations'' of
  $p(\tau|\ul y)$, learn at rate $O(1/n^{k/d})$.
\end{itemize}

\item Monte Carlo (MC) methods have long been used to help with the
curse.  With a random sample $\tau^{(1)}, \tau^{(2)}, \ldots,
\tau^{(n)}$, can plot histograms (shape), compute averages
(integration), etc.
\begin{itemize}
\item Improves convergence to $O(1/\sqrt{n})$
\item Less smoothing of interesting features as $n$ grows
\item Curse returns: how do you sample from
$p(\tau|\ul y)$?
\end{itemize}

\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{MCMC: The Cliff's Notes}


{\em The problem:}
\begin{itemize}
\item Learn about $\pi(\tau)=\pi(\theta,\beta,\lambda)$, where
$\tau=(\theta,\beta,\lambda)$ is some high-dimensional set of
variables (parameters).
\end{itemize}

{\em The essential idea:}
\begin{itemize}

\item Define a (stationary) Markov chain $\M_0, \M_1, \M_2, \ldots$ 
with states $\M_\nu=(\theta^\nu,\beta^\nu,\lambda^\nu)$.

\item Simulate the $\M_{\nu}$'s; under regularity conditions
(e.g., Tierney, 1994), $\M_{\nu}$ will converge in distribution to a
{\em stationary distribution} $\pi(\theta,\beta,\lambda)$:
% , as $\nu$ grows.
\begin{eqnarray*}
        \pi(\tone,\bone,\lone)  & = &  
%       \lefteqn{\pi(\tone,\bone,\lone)}  \\ & = & 
        \int P[\M_{\nu+1} = (\tone, \bone, \lone)~|~\M_{\nu}
        = (\tzero, \bzero, \lzero)] ~\times~       % \\ & & ~~~~~~\times~
        \,\pi(\tzero,\bzero,\lzero)\, d(\tzero,\bzero,\lzero)  
\end{eqnarray*}

\item For Bayes, design chain so that $\pi(\theta, \beta, \lambda)$
turns out to be the posterior, $p(\theta, \beta,
\lambda|\ul y)$.

% \item The behavior of the Markov chain is determined by its transition
% kernel 
% \[
% %     t[(\tzero, \bzero, \lzero), (\tone, \bone, \lone)] =
%       P[\M_{\nu+1} = (\tone, \bone, \lone)~|~\M_{\nu} = (\tzero,
%       \bzero, \lzero)]
% \]

\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\em Typical MCMC paradigm:}

\begin{itemize}

\item Write $p(\theta, \beta,\lambda| \ul y) = p(\tau_1,\tau_2, \ldots,
\tau_M|\ul y)$,
where $(\tau_1, \ldots, \tau_M)$ is a disjoint partition of the
original model parameters $(\theta,\beta,\lambda)$.

\item General theory of MCMC (e.g., Tierney, 1994; Chib and Greenberg,
1995): construct state $\M_\nu$ $=$ $(\tau_1^{(\nu)}, \ldots,
\tau_M^{(\nu)})$, by sampling each $\tau_m$ from its ``complete
conditional'' distribution:

\hfil\fbox{\parbox{6in}{\small\em
To step from $\M_{\nu-1}=(\tau_1^{(\nu-1)}, \ldots, \tau_M^{(\nu-1)})$ to
$\M_{\nu}=(\tau_1^{(\nu)}, \ldots, \tau_M^{(\nu)})$:
\begin{enumerate}
\item $\tau_1^{(\nu)}$ $\sim$ $p(\tau_1 | \tau_2^{(\nu-1)}, \ldots,
  \tau_M^{(\nu-1)}, \ul y)$;
\item $\tau_2^{(\nu)}$ $\sim$  $p(\tau_2 | \tau_1^{(\nu)},
  \tau_3^{(\nu-1)}, \ldots, \tau_M^{(\nu-1)},\ul y)$;
\item $\tau_3^{(\nu)}$ $\sim$ $p(\tau_3 | \tau_1^{(\nu)},
  \tau_2^{(\nu)}, \tau_4^{(\nu-1)}, \ldots, \tau_M^{(\nu-1)}, \ul y)$;

\hfil\vdots\hfil\vdots\hfil\vdots\hfil\vdots

\item[M.] $\tau_M^{(\nu)}$ $\sim$ $p(\tau_M | \tau_1^{(\nu)},
  \tau_2^{(\nu)}, \tau_3^{(\nu)}, \ldots, \tau_{M-1}^{(\nu)}, \ul y,)$;
\end{enumerate}
}}

This is the basic idea behind all so-called ``Gibbs samplers.''  

\end{itemize}

% \begin{itemize}
% 
% \item 
Markov Chain Monte Carlo (MCMC) vs. MC:
\begin{itemize}
\item Replace draws from $p(\tau|\ul y)$ with a sequence of
low-dimensional draws from $p_1(\tau_1)$, $p_2(\tau_2)$, etc.
\item For Bayes, often can exploit ``partial''
conjugacy to simplify problem of drawing from  $p(\tau_{(k)})$.
\item {\em Costs:} dependent rather than iid samples, approximate rather
  than exact sampling of $p(\tau|\ul y)$.
\end{itemize}


% \item Costs: 
% \begin{itemize}
% \item $\tau$'s from MC are iid, often {\em directly} from $p(\tau|\ul
%   y)$.  Easy statistics, ``fast'' convergence.
% \item $\tau$'s from MCMC are dependent, often from {\em approximation} to
%   $p(\tau|\ul y)$ .  Harder statistics, ``slow''
% convergence.
% \end{itemize} 
% \end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Implementation Details}

% \begin{itemize}
% \item The $\tau_m$'s may be univariate or multivariate
% \item The order in which the $\tau_m$'s are sampled is usually
% irrelevant
% \end{itemize}

\subsection*{Sampling schemes}

Once a set of complete conditionals is decided upon:
\begin{itemize}
\item {\em If you can} sample $\tau_m$ directly from its complete
  conditional distribution, do so---this is a \ul{\em Gibbs step};
\item {\em Otherwise} perform a \ul{\em Metropolis step}:
\begin{itemize}
\item Sample a ``proposal value'' $\tau_m^*$ from any convenient
  proposal distribution $q_\nu(\tau_m|\tau_m^{(\nu-1)})$
\item Compute the ``acceptance probability'' $\alpha_m^*$ from
  $q_\nu(\tau_m|\tau_m^{(\nu-1)})$ and the complete conditional
  $p(\tau_m | \rest)$ $\equiv$ $p(\tau_m$ $|$ $\tau_1^{(\nu)},$
  $\ldots,$ $\tau_{m-1}^{(\nu)},$ $\tau_{m+1}^{(\nu-1)},$ $\ldots,$
  $\tau_M^{(\nu-1)},$ $\ul y)$:
\[
        \alpha_n^* = \min\left\{
            \frac{p(\tau_n^*|\rest)q_\nu(\tau_n^{(\nu-1)}|\tau_n^*)}
            {p(\tau_n^{(\nu-1)}|\rest)q_\nu(\tau_n^*|\tau_n^{(\nu-1)})}
            , 1 \right\} ;
\]

\item Accept $\tau_n^{(\nu)} = \tau_n^*$ with probability
  $\alpha_n^*$; otherwise set $\tau_n^{(\nu)} = \tau_n^{(\nu-1)}$.
\end{itemize}
\end{itemize}
This is often called ``Metropolis-Hastings within Gibbs''.


\subsection*{Identifying the complete conditionals}


Suppose the $\ts$, $\bs$ and $\lambda$'s have prior distributions
$p(\theta|\lambda)$, $p(\beta)$ and $p(\lambda)$.  Then the ``complete
conditional'' for $\theta$, for example, is

\vfill

\begin{eqnarray*}
        p(\theta|\rest) ~=~
        p(\theta|\ul y,\beta,\lambda) 
        & = & \frac{p(\ul y,\theta,\beta,\lambda)}
        {\int p(\ul y,t,\beta,\lambda)\,dt}\\
        & = & \frac{p(\ul y|\theta,\beta)\,p(\theta|\lambda)
        \,\,p(\beta)\,\,p(\lambda)}
        {\int p(\ul y|t,\beta)\,p(t|\lambda)\,p(\beta)\,p(\lambda)\,dt}\\
        & \propto & p(\ul y|\theta,\beta)\,p(\theta|\lambda)
\end{eqnarray*}

\hfil\fbox{\parbox{6in}{\bf Key observation: \em The shape of
$p(\theta|\ul y,\beta,\lambda)$ is determined by just the parts of the
likelihood that depend explicitly on $\theta$.}}

\clearpage

\subsection*{Convergence}

\[
      \underbrace{\M_0, \M_1, \M_2, \ldots, \M_k,}_{\mbox{``Burn-in'' segment}}
      ~~~~
      \underbrace{\M_{k+1}, \ldots, \M_{k+n}}_{\mbox{Usable MCMC sample}}
\]

How large should $k$ be?
\begin{itemize}
\item Use {\tt acf()} in Splus/R to check when
  $Corr(\M_{\ell},\M_{\ell+k})$ small;
\item Use time series plots to see when the chain has ``stabilized''.
\item Run several chains from different starting points; let $k$ be so
  large that {\em (within-chain variation) $<$ (between chain variation)}.
\item {\tt CODA} subroutine package for Splus/R offers a menu of such
  checks \\ (www.mrc-bsu.cam.ac.uk/bugs/).
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Example: The hierarchical beta-binomial model}

\subsection*{Non-hierarchical version}

Basic beta-binomial model, for examinees $i=1,\ldots, N$, each getting $y_i$
of $n_i$ questions right.  [$\ul y=(y_1,\ldots,y_N)$]
\begin{enumerate}
\item $y_i | \theta_i \sim Bin(\theta_i,n_i)$;
\item $\theta_i | \alpha, \beta\sim Beta(\alpha,\beta)$
\end{enumerate}
Full model $p(\ul y, \ul \theta|\alpha,\beta)$ is
\[
        \prod_{i=1}^N \combos{n_i}{y_i} \theta_i^{y_i}(1-\theta_i)^{n_i-y_i}
        \prod_{i=1}^N \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} 
        \theta_i^{\alpha-1} (1-\theta_i)^{\beta-1}
\]
Use the ``key observation'' to calculate
\[
        p(\theta_i|y_i,\alpha,\beta) \propto \theta_i^{\alpha+y_i-1} 
        (1-\theta_i)^{\beta+n_i-y_i-1}
\]
i.e. $\theta_i|y_i \sim Beta(\alpha+y_i,\beta+n_i-y_i)$.

Choose fixed values for $\alpha$, $\beta$, then use
$Beta(\alpha+y_i,\beta+n_i-y_i)$ for posterior mode, posterior mean,
credible interval for $\theta$, etc. 

No need for anything fancy yet.


\subsection*{Hierarchical version}


What if we want to estimate $\alpha$, $\beta$ (estimate the shape of
the population / latent distribution)?  Add a third modeling assumption:
\begin{enumerate}
\item[3.] $\alpha,\beta \sim p(\alpha,\beta)$
\end{enumerate}
Now the full model $p(\ul y, \ul \theta, \alpha, \beta)$ is
{\small
\[
        \prod_{i=1}^N \combos{n_i}{y_i} \theta_i^{y_i}(1-\theta_i)^{n_i-y_i}
        \prod_{i=1}^N \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} 
        \theta_i^{\alpha-1} (1-\theta_i)^{\beta-1}
        p(\alpha,\beta)
\]
}%
Can again use the ``key observation'' to calculate
{\small
\[
\begin{array}{r@{~}c@{~}l}
        p(\theta_i|y_i,\alpha,\beta) & = & Beta(\theta_i|\alpha+y_i,\beta+n_i-y_i)\\
        p(\alpha,\beta|\ul y) & \propto  & p(\alpha,\beta)
        \displaystyle\prod_{i=1}^N 
        \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} 
        \frac{\Gamma(\alpha+y_i)\Gamma(\beta+n_i-y_i)}{\Gamma(\alpha+\beta+n_i)}
\end{array}
\]
}% 
% Gelman et al. (1995, pp. 131--133) suggest computing
% $p(\alpha,\beta|\ul y)$---or actually
% $p(\log(\alpha/\beta),\log(\alpha+\beta)|\ul y)$---on a grid, using
% trial and error to place the grid over the ``interesting'' part of the
% density.

\subsection*{An MCMC solution}

From the full model $p(\ul y, \ul \theta, \alpha, \beta)$ is
{\small
\[
        \prod_{i=1}^N \combos{n_i}{y_i} \theta_i^{y_i}(1-\theta_i)^{n_i-y_i}
        \prod_{i=1}^N \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} 
        \theta_i^{\alpha-1} (1-\theta_i)^{\beta-1}
        p(\alpha,\beta)
\]
}%
the ``key observation'' gives the complete conditionals
{\small
\begin{eqnarray*}
        p(\theta_i|\rest) & = & Beta(\theta_i|\alpha + y_i,\beta + n_i-y_i) \\
        p(\alpha|\rest) & \propto & 
        \left[\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)}
        \right]^N \prod_{i=1}^N \theta_i^\alpha p(\alpha,\beta) \\
        p(\beta|\rest) & \propto & 
        \left[\frac{\Gamma(\alpha+\beta)}{\Gamma(\beta)}
        \right]^N \prod_{i=1}^N (1-\theta_i)^\beta p(\alpha,\beta) 
\end{eqnarray*}
}%
\begin{itemize}
\item Usually take $p(\alpha,\beta) = p(\alpha) p(\beta)$; in this
  case we could take independent Gamma or Inv-Gamma densities, e.g.
\item Canned routines sample $\theta_i\sim Beta(\cdots)$
\item Can set up Metropolis steps for $\alpha$ and $\beta$ using
Normal proposal draws (``random walk M-H'').  Normal variances are
``tuning parameters''. 
\end{itemize} 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% \clearpage
% 
% {\em A brief comparison:}
% 
% \vfill
% 
% \begin{itemize}
% 
% \item Direct, grid-based contour plots and sampling (Gelman et
% al. 1995, pp. 128ff.)
% 
% \vfill
% 
% \item MCMC
% 
% \end{itemize}
% 
% \vfill
% 
% 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
