% From SamplerCompare, (c) 2010-2011 Madeleine Thompson

%\VignetteIndexEntry{Introduction to SamplerCompare}

\documentclass[nojss,article]{jss}
%\usepackage{Sweave}
\usepackage{array}
\usepackage{verbatim}

\def\svndate$#1: #2-#3-#4 #5 #6 (#7) ${\def\fmtdate{#2--#3--#4}}
\svndate$Date: 2011-08-23 21:20:33 -0400 (Tue, 23 Aug 2011) $

\title{Introduction to \pkg{SamplerCompare}}
\author{Madeleine B. Thompson\\University of Toronto}
\Plainauthor{Madeleine B. Thompson}
\date{\fmtdate}

\Abstract{
  \pkg{SamplerCompare} is an \proglang{R} package for comparing the
  performance of Markov chain Monte Carlo samplers.  It samples
  from a collection of distributions with a collection of MCMC
  methods over a range of tuning parameters.  Then, using log density
  evaluations per uncorrelated observation as a figure of merit,
  it generates a grid of plots showing the results of the simulation.
  It comes with a collection of predefined distributions and samplers
  and provides \proglang{R} and \proglang{C} interfaces for defining
  additional ones.  It also provides the means to import simulation
  data generated by external systems.  This document provides
  background on the package and demonstrates the basics of running
  simulations, visualizing results, and defining distributions and
  samplers in \proglang{R}.
}
\Keywords{MCMC, visualization}

%% \Volume{...}
%% \Issue{...}
%% \Month{...}
%% \Year{...}
%% \Submitdate{...}
%% \Acceptdate{...}

\Address{
  Madeleine B. Thompson\\
  E-mail: \email{madeleineth@gmail.com}\\
  URL: \url{https://mdln.net}
}

\SweaveOpts{eval=TRUE}

\begin{document}

<<echo=FALSE, results=hide>>=
dir.create('Figures', showWarnings=FALSE)
capture.output(source('Code/ex-compare-samplers.Rfrag'),
  file='Figures/ex-compare-samplers.txt')
capture.output(source('Code/ex-compare-results.Rfrag'),
  file='Figures/ex-compare-results.txt')

pdf('Figures/ex-comparison-plot.pdf', height=4, width=4)
print(comparison.plot(sampler.comparison, base_size=8))
dev.off()

source('Code/ex-metropolis.Rfrag')
source('Code/ex-beta.Rfrag')
source('Code/ex-final.Rfrag')
capture.output(source('Code/ex-final-compare.Rfrag'),
  file='Figures/ex-final-compare.txt')
@

\noindent A version of this document was published as  % remove-for-jss
\citet{sc-intro-jss}.                                  % remove-for-jss

\section{Purpose of package}

When a researcher develops a new Markov chain Monte Carlo (MCMC)
method, they will wish to determine how it compares to existing methods on
a representative set of distributions.  Similarly, when a statistician
specifies a new distribution, they may want to know which common
MCMC methods are most efficient at sampling from it.  \pkg{SamplerCompare}
\citep{samplercompare} is an \proglang{R} \citep{r} package that
automates these tasks.  It draws samples from a collection of
probability distributions with a collection of MCMC samplers, with a
range of tuning parameters, and presents the results of such simulations
graphically.  These comparisons allow researchers to better understand
which MCMC methods perform best in which circumstances.

The main goal of \pkg{SamplerCompare} is to generate a grid of plots
in which an individual grid cell corresponds to a single MCMC sampler
and a single distribution.  In each grid cell, the efficiency of a
simulation is summarized by a plot of the product of the number of
log density evaluations per iteration multiplied by the autocorrelation
time of the slowest-mixing component of the state space for a range
of scale tuning parameter values.  Autocorrelation time, like
effective sample size, which can be computed by \code{effectiveSize} in
the \pkg{coda} \citep{plummer10} package, accounts for the
often-substantial correlation between observations in MCMC-generated
samples.  Log density evaluations accounts for elapsed processor
time in a machine-independent way.  By viewing the cost measure in
a grid of plots, a user can see patterns in the performance more
easily than they could if they viewed the same results as numbers
in a table.

\section{Installation and documentation}

Binary packages for MacOS and Windows and a platform-independent
source package can be obtained from CRAN at:
\begin{center}
\url{http://cran.r-project.org/web/packages/SamplerCompare/index.html}
\end{center}
To install it, one must first install the \pkg{mvtnorm} package
\citep{mvtnorm}.  To use \pkg{SamplerCompare}'s graphics, one must
install the \pkg{ggplot2} package \citep{wickham09}.  To use
multithreading, one must install the
\pkg{synchronicity} \citep{synchronicity} package.  It
is not available for Windows, so Windows users are limited to
single-threaded simulations.

Alternatively, one can use the \code{install.packages} \proglang{R} command:
\begin{verbatim}
R> install.packages("SamplerCompare", dependencies=c("Depends","Suggests"))
\end{verbatim}

More information on \pkg{SamplerCompare} is available in the
\proglang{R} online help for the package and \citet{rcglue}.  After
the package is installed, a list of online help topics and vignettes
can be found by typing:
\begin{verbatim}
R> library(help='SamplerCompare')
\end{verbatim}
Vignettes can be read with the \code{vignette} command.  For example:
\begin{verbatim}
R> vignette('glue')
\end{verbatim}
PDF copies can be found in the \texttt{doc} directory of the installed
package.  Further information on the mathematical background of the
comparisons and analysis of the plots is available in
\citet{corlen}.

\section{Simulations with included samplers and distributions}

\begin{table}[p]
\begin{tabular}{>{\tt}lp{3.40in}}
\rm \proglang{R} function & Sampler \\ \hline
multivariate.metropolis.sample & Metropolis--Hastings with spherically
  symmetric Gaussian proposals \\
univar.metropolis.sample & Metropolis--Hastings with single-coordinate
  updates \\
adaptive.metropolis.sample & Adaptive Metropolis--Hastings \citep{roberts09} \\
arms.sample & Adaptive Rejection Metropolis \citep{gilks95} \\
stepout.slice.sample & slice sampler with stepping out \citep[\S4]{neal03} \\
interval.slice.sample & slice sampler without stepping out
  \citep[\S4]{neal03} \\
univar.eigen.sample & adaptive slice sampler with univariate steps
  along eigenvectors of covariance matrix \citep[ch.~3]{thompson-thesis} \\
hyperrectangle.sample & slice sampler with hypercube for initial
  slice approximation, shrinkage using gradient \citep[\S5.1]{neal03} \\
nograd.hyperrectangle.sample & slice sampler with hypercube for
  initial slice approximation, shrinkage in all dimensions
  \citep[\S5.1]{neal03} \\
oblique.hyperrect.sample & adaptive slice sampler with hyperrectangle
  for initial slice approximation \citep[ch.~3]{thompson-thesis} \\
nonadaptive.crumb.sample & slice sampler with Gaussian crumbs
  \citep[\S5.2]{neal03} \\
cov.match.sample & covariance-matching slice sampler \citep[\S4]{thompson10} \\
shrinking.rank.sample & shrinking rank slice sampler \citep[\S5]{thompson10} \\
\end{tabular}
\caption{Predefined samplers; see the \proglang{R} help for the
sampler's \proglang{R} function for more information on an individual
method.}
\label{samplers}
\end{table}

\begin{table}[p]
\begin{tabular}{>{\tt}lp{4in}}
\rm \proglang{R} symbol & Distribution \\ \hline
N2weakcor.dist & weakly correlated two-dimensional Gaussian \\
N4poscor.dist & strongly positively correlated four-dimensional Gaussian \\
N4negcor.dist & strongly negatively correlated four-dimensional Gaussian \\
schools.dist & ten-dimensional multilevel model
  \citep[pp.~138--145]{gelman04} \\
funnel.dist & ten-dimensional distribution with funnel-shaped
  marginals \citep[p.~732]{neal03} \bigskip\\
\rm \proglang{R} function & Distributions generated\\ \hline
make.gaussian & multivariate Gaussians \\
make.cone.dist & distributions with cone-shaped log density
  \citep{roberts02} \\
make.multimodal.dist & mixtures of standard Gaussians \\
make.mv.gamma.dist & distributions with uncorrelated gamma marginals \\
\end{tabular}
\caption{Predefined distributions and functions that generate
distributions; see the \proglang{R} help for a symbol for more
information on an individual distribution or generator.}
\label{dists}
\end{table}

The three central types of objects in \pkg{SamplerCompare} are
distributions (which have the class \code{scdist}), sampler functions, and
simulation results.  The function \code{compare.samplers} runs a list
of samplers on a list of distributions with a set of tuning parameters
and returns a data frame containing simulation results.  Sampler
functions are assumed to have a single scalar tuning parameter.  If
they have more, wrapper functions can represent a single sampler
with a varying tuning parameter as multiple samplers.  \pkg{SamplerCompare}
comes with a collection of predefined samplers (listed in
Table~\ref{samplers}) and distributions (listed in Table~\ref{dists}).

Suppose we would like to compare Adaptive Metropolis
(\code{adaptive.metropolis.sample}) and Adaptive Rejection Metropolis
(\code{arms.sample}) with the tuning parameters 0.1, 1, 10, and 100 on
two-dimensional Gaussian (\code{make.gaussian}) and Gamma (\code{make.mv.gamma.dist})
distributions.  We can do this with \code{compare.samplers} using the
\proglang{R} code:
\verbatiminput{Code/ex-compare-samplers.Rfrag}

The call to \code{compare.samplers} generates the following trace, with
one line for each simulation:
\verbatiminput{Figures/ex-compare-samplers.txt}
Each line in the trace has the distribution name, the sampler name,
the number of evaluations per uncorrelated observation with 95\%
confidence interval in parentheses, the tuning parameter, and the
autocorrelation time of the log density.

The return value of \code{compare.samplers} (\code{sampler.comparison} in
this example) is a data frame with one row per simulation.  To see,
for example, how many evaluations and how many processor seconds
Adaptive Metropolis needed to generate an uncorrelated observation
on the two-dimensional Gaussian, one would multiply the \code{act} column
by the \code{evals} and \code{cpu} columns:
\verbatiminput{Code/ex-compare-results.Rfrag}
This generates the output:
\verbatiminput{Figures/ex-compare-results.txt}
Notice that the first product is the same as the one reported in
the corresponding line from the \code{compare.samplers} trace output
(the seventh, \texttt{N2,rho=0.7 Adaptive Metropolis}~.~.~.
\texttt{tuning=1}).  See the \proglang{R} help page for
\code{sampler.comparison} for a full list of the columns of its return
value.

\section{Visualizing results}
\label{visualizing}

To visually compare the efficiency of a collection of simulations,
one can use the \texttt{comparison.\discretionary{}{}{}plot} function.
It has a single required argument, a data frame containing results
from \code{compare.samplers} or \code{simulation.result}.  The previous section
has an example of the use of \code{compare.samplers}.  The online help
for \code{simulation.result} contains an example of its use to load a
simulation generated by \pkg{JAGS} \citep{jags220}.

\code{comparison.plot} returns a \pkg{ggplot2} plot object.  One can
call \code{print} on this object to view the plot; it can also be edited
with the \pkg{grid} package \citep[ch.~5--6]{murrell05}.  To plot
the results from the previous example, one would type:

\begin{verbatim}
R> print(comparison.plot(sampler.comparison))
\end{verbatim}
The results are shown in Figure~\ref{ex-comparison-plot}.

\begin{figure}
\begin{center}
\includegraphics{Figures/ex-comparison-plot.pdf}
\end{center}
\caption{A comparison between Adaptive Metropolis and ARMS on
two-dimensional Gaussians and gammas.  See section~\ref{visualizing}
for discussion.}
\label{ex-comparison-plot}
\end{figure}

\noindent In this figure, the columns of plots represent the
samplers, and the rows of plots represent the distributions.  The
scale tuning parameter is plotted on the horizontal
axis; the number of log density evaluations per uncorrelated
observation is plotted as a dot on the vertical axis with a bar for the 95\%
confidence interval.  Log density evaluations per uncorrelated
observation is computed by multiplying the average number of log
density evaluations per simulation iteration by the autocorrelation
time of the slowest-mixing component of the simulation.  The
autocorrelation time is the ratio of the sample size to the effective
sample size. It accounts for linear dependence between successive
states.  See the help for \code{ar.act} for more information on how it
is computed.

In Figure~\ref{ex-comparison-plot}, one can see that on the two
distributions compared, when sampling with ARMS, the cost measure
does not vary much with the tuning parameter.  Adaptive Metropolis seems
more sensitive to the tuning parameter. However, when the components
of the target distribution are correlated, it performs better than
ARMS when the tuning parameter is well chosen (in this case, equal
to one).  For more discussion of the interpretation of these plots,
see \citet[\S5]{corlen}.

\section{Defining a sampler}
\label{sec-sampler}

MCMC samplers are specified by functions that have the signature:
\begin{verbatim}
sampler(target.dist, x0, sample.size, tuning)
\end{verbatim}
They must also have a \code{name} attribute, a length-one character
vector.  The \code{target.dist} parameter specifies the target distribution;
see the \proglang{R} help for \code{make.dist} for details on its
structure.  \code{x0} specifies the start state for the simulation,
\code{sample.size} specifies the sample size, and \code{tuning} specifies a
scalar tuning parameter.

A sampler function should return a list with two elements: \code{X}, a
matrix with one row per observation, and \code{evals}, a count of the number
of times it evaluated the log density (with \code{target.dist$log.density}).
If the sampler evaluates the gradient of the log density (with
\code{target.dist$grad.log.density}), the list should contain a \code{grads}
element, indicating the number of times it did this.

The following code specifies a Metropolis sampler with multivariate
proposals:
\verbatiminput{Code/ex-metropolis.Rfrag}

See the \proglang{R} help for \code{compare.samplers} for more information
on writing samplers in \proglang{R}.  See the \proglang{R} help for
\code{wrap.c.sampler} and \citet{rcglue}
for more information on writing samplers in \proglang{C}.

\section{Defining a distribution}
\label{sec-dist}
\nobreak
\code{make.dist} can be used to specify a distribution whose log density
is expressed in \proglang{R}.  (See the \proglang{R} help for
\code{make.c.dist} and \citet{rcglue} for more information on specifying
distributions in \proglang{C}.)  Its most important arguments are
\code{ndim}, \code{name}, and \code{log.density}.  \code{ndim} specifies the dimension
of the distribution and \code{name} names the distribution.  \code{log.density}
is a function of one vector argument of length \code{ndim} that returns
the log density at that point; it should return \code{-Inf} if the point
is outside the support of the distribution.  The log density does
not need to be normalized.

The following \proglang{R} code defines a Beta(2,3) distribution:
\verbatiminput{Code/ex-beta.Rfrag}
The optional \code{mean} argument to \code{make.dist} makes the autocorrelation
time computation in \code{compare.samplers} more accurate, so it is advisable
to specify it when the mean is known.

\section{Comparing user-defined distributions and samplers}

User-defined samplers and distributions can be used just like the
included samplers and distributions.  To use the Metropolis sampler
defined in section~\ref{sec-sampler} to sample from the beta
distribution defined in section~\ref{sec-dist}, one can run the
code from those two sections and then run:
\verbatiminput{Code/ex-final.Rfrag}
User-defined samplers can also be used with \code{compare.samplers}:
\verbatiminput{Code/ex-final-compare.Rfrag}
The call to \code{print(subset(...))} shows some of the columns of the
result object:
\verbatiminput{Figures/ex-final-compare.txt}
One can see that since the evaluations per iteration (\code{evals}) and
processor-seconds per iteration (\code{cpu}) are similar for each
simulation, and the autocorrelation time (\code{act}) is lowest for a
tuning parameter of $1.0$, that choice would seem to be better than
the other two.  However, the plots produced by \code{comparison.plot}
are easier to interpret when more than a few chains are run.

\section{Limitations}

\pkg{SamplerCompare} was created to support my own research; I am
releasing it with the hope that others find it useful.  Some current
limitations include:
\begin{itemize}
\item Distributions are assumed to be continuous and to be of a
  constant dimension.
\item Samplers are assumed to have exactly one scalar tuning parameter.
\item All samplers in a given invocation of \code{compare.samplers} are
  run with the same simulation length and set of tuning parameters.
\item Distributions are defined entirely in terms of their log
  density; there is no way to specify that a distribution is
  unimodal or that a particular parameter is always positive.
\item Multithreading is not supported on Windows.
\end{itemize}

\bibliography{jss668}

\end{document}