Introduction to SamplerCompare

SamplerCompare is an R package for comparing the performance of Markov chain Monte Carlo (MCMC) 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 ﬁgure of merit, it generates a grid of plots showing the results of the simulation. It comes with a collection of predeﬁned distributions and samplers and provides R and C interfaces for deﬁning 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 deﬁning distributions and samplers in R .


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. SamplerCompare (Thompson 2011a) is an R (R Development Core Team 2011) 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 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 effectiveSize in the coda (Plummer et al. 2006) 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.

Installation and documentation
Binary packages for MacOS and Windows and a platform-independent source package can be obtained from the Comprehensive R Archive Network at http://CRAN.R-project.org/ package=SamplerCompare. To install it, one must first install the mvtnorm package (Genz et al. 2011). To use SamplerCompare's graphics, one must install the ggplot2 package (Wickham 2009). To use multithreading, one must install the multicore (Urbanek 2011) and synchronicity (Kane 2011) packages. These two are not available for Windows, so Windows users are limited to single-threaded simulations.
Alternatively, one can use the install.packages R command: R> install.packages("SamplerCompare", + dependencies = c("Depends", "Suggests")) More information on SamplerCompare is available in the R online help for the package and Thompson (2011b). After the package is installed, a list of online help topics and vignettes can be found by typing library(help = "SamplerCompare"). Vignettes can be read with the vignette command, e.g., vignette("glue"). PDF copies can be found in the doc directory of the installed package. Further information on the mathematical background of the comparisons and analysis of the plots is available in Thompson (2010).

Simulations with included samplers and distributions
The three central types of objects in SamplerCompare are distributions (which have the class dist), sampler functions, and simulation results. The function 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. SamplerCompare comes with a collection of predefined samplers (listed in Table 1) and distributions (listed in Table 2).
Suppose we would like to compare adaptive Metropolis (adaptive.metropolis.sample) and adaptive rejection Metropolis (arms.sample) with the tuning parameters 0.1, 1, 10, and 100 on two-dimensional Gaussian (make.gaussian) and Gamma (make.mv.gamma.dist) distributions. We can do this with compare.samplers using the R following code and generating the trace below, with one line for each simulation: R> library("SamplerCompare") R> gauss.cor7 <-make.gaussian(mean = c(1, 2), rho = 0.7) R> gamma.shape23 <-make.mv.gamma.dist(shape = c(2, 3))  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 compare.samplers (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 act column by the evals and cpu columns: Notice that the first product is the same as the one reported in the corresponding line from the compare.samplers trace output (the seventh, N2,rho=0.7 Adaptive Metropolis . . . tuning=1). See the R help page for sampler.comparison for a full list of the columns of its return value.
The returned data is also incrementally saved to the file named by completed.file. It can be loaded in a second R process with load to see how the simulation is progressing. If unset, a temporary file is created for this purpose and deleted on successful completion of the simulation.

Visualizing results
To visually compare the efficiency of a collection of simulations, one can use the comparison. plot function. It has a single required argument, a data frame containing results from compare.samplers or simulation.result. The previous section has an example of the use of compare.samplers. The online help for simulation.result contains an example of its use to load a simulation generated by JAGS (Plummer 2003).
comparison.plot returns a ggplot2 plot object. One can call print on this object to view the plot; it can also be edited with the grid package (Murrell 2005, Chapter 5-6). To plot the results from the previous example (see Figure 1), one would type: 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 ar.act for more information on how it is computed.
In Figure 1, 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 Thompson (2010, Section 5).

Defining a sampler
MCMC samplers are specified by functions that have the signature: sampler(target.dist, x0, sample.size, tuning) They must also have a name attribute, a length-one character vector. The target.dist parameter specifies the target distribution; see the R help for make.dist for details on its structure. x0 specifies the start state for the simulation, sample.size specifies the sample size, and tuning specifies a scalar tuning parameter.
A sampler function should return a list with two elements: X, a matrix with one row per observation, and evals, a count of the number of times it evaluated the log density (with target.dist$log.density). If the sampler evaluates the gradient of the log density (with target.dist$grad.log.density), the list should contain a grads element, indicating the number of times it did this.
The following code specifies a Metropolis sampler with multivariate proposals: return(list(X = X, evals = evals)) + } R> attr(metropolis.sample, "name") <-"Metropolis" See the R help for compare.samplers for more information on writing samplers in R. See the R help for wrap.c.sampler and Thompson (2011b) for more information on writing samplers in C.

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

Limitations
SamplerCompare was created to support my own research; I am releasing it with the hope that others find it useful. Some current limitations include: Distributions are assumed to be continuous and to be of a constant dimension.
Samplers are assumed to have exactly one scalar tuning parameter.
All samplers in a given invocation of compare.samplers are run with the same simulation length and set of tuning parameters.
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.
Multithreading is not supported on Windows.