SamplerCompare (version 1.3.4)

simulation.result: Summarize one MCMC chain

Description

Summarize one MCMC chain in the format used by compare.samplers

Usage

simulation.result(target.dist, sampler.name, X,
                  evals=NULL, grads=NULL, tuning=NULL, cpu=NULL,
                  burn.in=0.2, y=NULL,
                  sampler.expr=sprintf("plain('%s')", sampler.name),
                  aborted=NA)

Value

A single-row data frame of the format returned by

compare.samplers.

Arguments

target.dist

A distribution object of the sort generated by make.dist representing the distribution sampled from. This is used to obtain the dimension and name of the distribution; the log density function does not need to be specified.

sampler.name

The name of the sampler that generated this simulation. If generated by SamplerCompare, this would usually be the name attribute of the sampler function.

X

A matrix (or object that can be coerced to a matrix) containing the simulation results. It should have one row per iteration and one column for each component of the state space. Corresponds to the X element of the list returned by a sampler.

evals

The total number of log density evaluations used in the simulation; corresponds to the evals element of the list returned by a sampler.

grads

The total number of log density gradient evaluations used in the simulation; corresponds to the grads element of the list returned by a sampler.

tuning

The scalar tuning parameter passed to the sampler.

cpu

The processor time used to generate the simulation in seconds.

burn.in

Initial fraction of X to discard before computing autocorrelation times.

y

A vector with the same number of elements as X has rows containing the log densities at the states represented by those rows.

sampler.expr

The name of the sampler that generated this simulation in plotmath format. If generated by SamplerCompare, this would usually be the name.expression attribute of the sampler function.

aborted

A logical scalar indicating whether the simulation was prematurely aborted.

Details

This function summarizes a simulation into a single-row data frame by computing the autocorrelation time of its slowest-mixing component and, if possible, the autocorrelation time of the log density and the error in the sample mean. The autocorrelation time of the slowest-mixing component can always be estimated, but is more accurate if the true mean is specified in target.dist. The autocorrelation time of the log density can be estimated if either the log density function is specified in target.dist or an explicit vector of log densities is passed as y. The error in the sample mean can be computed if the mean is specified in target.dist.

This function is intended to be called once per simulation for a variety of simulations. The results are to be combined with rbind and can be visualized with comparison.plot. While the evals and tuning arguments are optional, the result cannot be used with comparison.plot if it is not set. simulation.result is normally called internally by compare.samplers but is exported so that simulations run in external systems such as JAGS can be analyzed with SamplerCompare. See the “Examples” section for an example of this usage.

References

Thompson, M. B. (2011), “Introduction to SamplerCompare,” Journal of Statistical Software 43(12):1-10, tools:::Rd_expr_doi("10.18637/jss.v043.i12").

See Also

compare.samplers, comparison.plot, “Introduction to SamplerCompare” (vignette)

Examples

if (FALSE) {
# An example generated with the following JAGS model:
#
# model {
#   mu[1] <- 0
#   mu[2] <- 0
#   Sigma[1,1] <- 1
#   Sigma[2,2] <- 1
#   Sigma[1,2] <- 0.7
#   Sigma[2,1] <- 0.7
#   x ~ dmnorm(mu, inverse(Sigma))
# }
#
# and the following JAGS script:
#
# model in "mv.7.model"
# compile, nchains(1)
# initialize
# update 1000
# monitor x
# update 10000
# coda *

# Load data written by JAGS

library(coda)
X <- read.coda('CODAchain1.txt', 'CODAindex.txt')

# Dummy distribution object.

N2.dist <- make.dist(2, '2D Normal, cor=0.7', mean=c(0,0))

# Compute simulation result.  evals and tuning are hacks; they
# are undefined with Gibbs sampling.  JAGS can do its own burn-in,
# so set burn.in to zero.

sim.result <- simulation.result(N2.dist, 'JAGS', X,
                                evals=nrow(X)*ncol(X), tuning=1,
                                burn.in=0)
}