\documentclass{article}

\usepackage{amstext}

% \VignetteIndexEntry{Debugging MCMC Code}

\begin{document}

<<foo,include=FALSE,echo=FALSE>>=
options(keep.source = TRUE, width = 60)
foo <- packageDescription("mcmc")
@

\title{Debugging MCMC Code}
\author{Charles J. Geyer}
\maketitle

\section{Introduction}

This document discusses debugging Markov chain Monte Carlo code
using the R contributed package \texttt{mcmc} (Version \Sexpr{foo$Version})
for examples.  It also documents the debugging output of the functions
\texttt{mcmc} and \texttt{temper}.

Debugging MCMC code if the code is taken as a black box is basically
impossible.  In interesting examples, the only thing one knows about
the equilibrium distribution of an MCMC sampler is what one learns
from the samples.  This obviously doesn't help with testing.  If the
sampler is buggy, then the only thing you know about the equilibrium
distribution is wrong, but if you don't know it is buggy, then you don't
know it is wrong.  So you don't know anything.  There is no way to tell
whether random output has the correct distribution when you don't know
anything about the distribution it is supposed to have.

The secret to debugging MCMC code lies in two principles:
\begin{itemize}
\item take the randomness out, and
\item expose the innards.
\end{itemize}
The first slogan means consider the algorithm a deterministic function
of the elementary pseudo-random numbers that are trusted (for example,
the outputs of the R random number generators, which you aren't responsible
for debugging and are also well tested).
The second slogan means output, at least for debugging purposes, enough
intermediate state so that testing is straightforward.

For a Gibbs sampler, this means outputting all of the trusted elementary
pseudo-random numbers used, the state before and after each elementary
Gibbs update, and which update is being done if a random scan is used.
Also one needs to output the initial seeds of the pseudo-random number
generator (this is true for all situations and will not be mentioned again).

For a Metropolis-Hastings sampler, this means outputting all of the trusted
elementary
pseudo-random numbers used, the state before and after each elementary
Metropolis-Hastings update, the proposal for that update, the Hastings ratio
for that update, decision (accept or reject) in that update.

For more complicated MCMC samplers, there is more ``innards'' to ``expose''
(see the discussion of the \texttt{temper} function below), but you get the
idea.  You can't output too much debugging information.

\section{The Metrop Function}

The R function \texttt{metrop} in the \texttt{mcmc} package has an argument
\verb@debug = FALSE@ that when \verb@TRUE@ causes extra debugging information
to be output.
Let \texttt{niter} be the number of iterations
\verb@nbatch * blen * nspac@, and let \texttt{d} be the dimension of the state
vector.  The result of invoking \texttt{metrop} is a list.  When
\verb@debug = TRUE@ it has the following additional components
\begin{itemize}
\item \texttt{current}, an \texttt{niter} by \texttt{d} matrix
    of mode \verb@"numeric"@, the state before iteration \texttt{i}
    is \verb@current[i, ]@
\item \texttt{proposal}, an \texttt{niter} by \texttt{d} matrix
    of mode \verb@"numeric"@, the proposal for iteration \texttt{i}
    is \verb@proposal[i, ]@
\item \texttt{z}, an \texttt{niter} by \texttt{d} matrix
    of mode \verb@"numeric"@, the vector of standard normal random variates
    used to generate the proposal for iteration \texttt{i}
    is \verb@z[i, ]@
\item \texttt{log.green}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the logarithm of the Hastings ratio
    for each iteration
\item \texttt{u}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the $\text{Uniform}(0, 1)$ random variate
    compared to the Hastings ratio for each iteration or \texttt{NA} if
    none is needed (when the log Hastings ratio is nonnegative)
\item \texttt{debug.accept}, a vector of length \texttt{niter}
    and mode \verb@"logical"@, the decision for each iteration,
    accept the proposal (\texttt{TRUE}) or reject it (\texttt{FALSE})
\end{itemize}
(The components \texttt{z} and \texttt{debug.accept} were added in
version 0.7-3 of the \texttt{mcmc} package.  Before that only the others
were output.)

Two components of the list returned by the \texttt{metrop} function always
(whether \verb@debug = TRUE@ or \verb@debug = FALSE@) are also necessary
for debugging.  They are
\begin{itemize}
\item \texttt{initial.seed} the value of the variable \texttt{.Random.seed}
    that contains the seeds of the R random number generator system before
    invocation of the \texttt{metrop} function
\item \texttt{final}, a vector of length \texttt{d} and mode \verb@"numeric"@,
    the state after the last iteration
\end{itemize}

All of the files in the \texttt{tests} directory of the source for the
package (not installed but found in the source tarball on CRAN) test
the \texttt{metrop} function except those beginning \texttt{temp},
which test the \texttt{temper} function.  Since these tests were written
many years ago, are spread out over many files, and are not commented,
we will not describe them in detail.  Suffice it to say that they check
every aspect of the functioning of the \texttt{metrop} function.

\section{The Temper Function}

The R function \texttt{temper} in the \texttt{mcmc} package has an argument
\verb@debug = FALSE@ that when \verb@TRUE@ causes extra debugging information
to be output.
Let \texttt{niter} be the number of iterations
\verb@nbatch * blen * nspac@, let \texttt{d} be the dimension of the state
vector, and let \texttt{ncomp} be the number of components of the tempering
mixture.
The result of invoking \texttt{temper} is a list.  When
\verb@debug = TRUE@ and \verb@parallel = TRUE@ it has the following additional
components
% which
% unif.which
% state
% log.hastings
% unif.hastings
% proposal
% acceptd
% norm
% unif.choose
% coproposal
\begin{itemize}
\item \texttt{which}, a vector of length \texttt{niter}
    and mode \verb@"logical"@ the type of update for each iteration,
    within component (\texttt{TRUE}) or swap components (\texttt{FALSE}).
\item \texttt{unif.which}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the $\text{Uniform}(0, 1)$ random variate
    used to decide which type of update is done.
\item \texttt{state}, an \texttt{niter} by \texttt{ncomp} by \texttt{d}
    array of mode \verb@"numeric"@, the state before iteration \texttt{i}
    is \verb@state[i, , ]@
\item \texttt{proposal}, an \texttt{niter} by \verb@d + 1@
    matrix of mode \verb@"numeric"@, the proposal for iteration \texttt{i}
    is \verb@proposal[i, ]@ (explanation below)
\item \texttt{coproposal}, an \texttt{niter} by \verb@d + 1@
    matrix of mode \verb@"numeric"@, the proposal for iteration \texttt{i}
    is \verb@coproposal[i, ]@ (explanation below)
\item \texttt{log.hastings}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the logarithm of the Hastings ratio for
    each iteration
\item \texttt{unif.hastings}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the $\text{Uniform}(0, 1)$ random variate
    compared to the Hastings ratio for each iteration or \texttt{NA} if
    none is needed (when the log Hastings ratio is nonnegative)
\item \texttt{acceptd}, a vector of length \texttt{niter}
    and mode \verb@"logical"@, the decision for each iteration,
    accept the proposal (\texttt{TRUE}) or reject it (\texttt{FALSE})
\item \texttt{norm}, an \texttt{niter} by \texttt{d} matrix
    of mode \verb@"numeric"@, the vector of standard normal random variates
    used to generate the proposal for iteration \texttt{i} is \verb@z[i, ]@
    unless none are needed (for swap updates) when it is \texttt{NA} 
\item \texttt{unif.choose}, an \texttt{niter} by 2 matrix
    of mode \verb@"numeric"@, the vector of $\text{Uniform}(0, 1)$
    random variates used to choose the components to update in iteration
    \texttt{i} is \verb@unif.choose[i, ]@; in a swap update two are used;
    in a within-component update only one is used and the second is \texttt{NA}
\end{itemize}

In a within-component update, one component say \texttt{j} is chosen for
update.  The \emph{coproposal} is the current value of the state for this
component, which is a vector of length \verb@d + 1@, the first
component of which is \texttt{j} and the rest of which is \verb@state[i, j, ]@
if we are in iteration \texttt{i}.
The \emph{proposal} is a similar vector, the first
component of which is again \texttt{j} and the rest of which is a multivariate
normal random vector centered at \verb@state[i, j, ]@.
The coproposal is the current state; the proposal is the possible value
(if accepted) of the state at the next time.

In a swap update, two components say \texttt{j1} and \texttt{j2} are chosen for
update.  Strictly, speaking the coproposal is the pair of vectors
\verb@c(j1, state[i, j1, ])@ and \verb@c(j2, state[i, j2, ])@
and the proposal is these swapped, that is, the pair of vectors
\verb@c(j2, state[i, j1, ])@ and \verb@c(j1, state[i, j2, ])@
if we are in iteration \texttt{i}.
Since, however, there is a lot of redundant information here,
the vector \verb@c(j1, state[i, j1, ])@ is output as \verb@coproposal[i, ]@
and the vector \verb@c(j2, state[i, j2, ])@ is output as \verb@proposal[i, ]@.

When \verb@debug = TRUE@ and \verb@parallel = FALSE@
the result of invoking \texttt{temper} is a list having
the following additional components
% which
% unif.which
% state
% log.hastings
% unif.hastings
% proposal
% acceptd
% norm
% unif.choose
\begin{itemize}
\item \texttt{which}, a vector of length \texttt{niter}
    and mode \verb@"logical"@ the type of update for each iteration,
    within component (\texttt{TRUE}) or jump from one component to
    another (\texttt{FALSE}).
\item \texttt{unif.which}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the $\text{Uniform}(0, 1)$ random variate
    used to decide which type of update is done.
\item \texttt{state}, an \texttt{niter} by \verb@d + 1@
    matrix of mode \verb@"numeric"@, the state before iteration \texttt{i}
    is \verb@state[i, ]@
\item \texttt{proposal}, an \texttt{niter} by \verb@d + 1@
    matrix of mode \verb@"numeric"@, the proposal for iteration \texttt{i}
    is \verb@proposal[i, ]@
\item \texttt{log.hastings}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the logarithm of the Hastings ratio for
    each iteration
\item \texttt{unif.hastings}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the $\text{Uniform}(0, 1)$ random variate
    compared to the Hastings ratio for each iteration or \texttt{NA} if
    none is needed (when the log Hastings ratio is nonnegative)
\item \texttt{acceptd}, a vector of length \texttt{niter}
    and mode \verb@"logical"@, the decision for each iteration,
    accept the proposal (\texttt{TRUE}) or reject it (\texttt{FALSE})
\item \texttt{norm}, an \texttt{niter} by \texttt{d} matrix
    of mode \verb@"numeric"@, the vector of standard normal random variates
    used to generate the proposal for iteration \texttt{i} is \verb@z[i, ]@
    unless none are needed (for jump updates) when it is \texttt{NA} 
\item \texttt{unif.choose}, a vector of length \texttt{niter}
    and mode \verb@"numeric"@, the $\text{Uniform}(0, 1)$
    random variates used to choose the component to update in iteration
    \texttt{i} is \verb@unif.choose[i, ]@; in a jump update one is used;
    in a within-component update none is used and \texttt{NA} is output
\end{itemize}

All of the files in the \texttt{tests} directory of the source for the
package (not installed but found in the source tarball on CRAN)
beginning \texttt{temp} test the \texttt{temper} function.
They check every aspect of the functioning of the \texttt{temper} function.

In the file \texttt{temp-par.R} in the \texttt{tests} directory, the following
checks are made according to the comments in that file
\begin{enumerate}
\item check decision about within-component or jump/swap
\item check proposal and coproposal are actually current state or part thereof
\item check hastings ratio calculated correctly
\item check hastings rejection decided correctly
\item check acceptance carried out or not (according to decision) correctly
\item check within-component proposal
\item check swap proposal
\item check standard normal and uniform random numbers are as purported
\item check batch means
\item check acceptance rates
\item check scale vector
\item check scale matrix
\item check scale list
\item check outfun
\end{enumerate}
In the file \texttt{temp-ser.R} in the \texttt{tests} directory, the all of
the same checks are made according to the comments in that file except for
check number 2 above, which would make no sense because there is no
\texttt{coproposal} component in the serial (\verb@parallel = FALSE@) case.

\end{document}