simulation.result: Summarize one MCMC chain
In SamplerCompare: A Framework for Comparing the Performance of MCMC Samplers

simulation.result

R Documentation

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)

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.

Value

A single-row data frame of the format returned by compare.samplers.

References

Thompson, M. B. (2011), “Introduction to SamplerCompare,” Journal of Statistical Software 43(12):1-10, \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v043.i12")}.

Examples

## Not run: 
# 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)

## End(Not run)

SamplerCompare documentation built on April 24, 2023, 9:09 a.m.