THE R/C PACKAGE 'AMCMC', BY JEFFREY S. ROSENTHAL, 2007 (updated March 2009)
Description:
Estimate the expected value of a functional with respect to a
multi-dimensional density function, by performing the "adaptive
Metropolis-within-Gibbs" MCMC algorithm of Section 3 of Roberts
and Rosenthal (2006).
Usage:
amcmc(densfn=dnorm, functional=firstcoordsq,
numbatches=1000, batchlength=10, logflag=FALSE, adapt=TRUE,
burninfrac = 0.1, initval = 1.0, verboselevel=0, cdensity=cfns,
cfunctional=cfns, cfns=FALSE, vectorlength=0, write=FALSE)
Arguments:
densfn: the target density function
functional: the real-valued function whose mean is sought
numbatches: the number of batches the algorithm will perform
batchlength: the number of full scans performed in each batch
logflag: if TRUE, densfn is treated as the logarithm of the density
function, rather than as the density function itself
adapt: if TRUE, the algorithm adapts the proposal variances as it runs
burninfrac: the fraction of initial output values to be discared
initval: the initial value of each variable (at time zero)
verboselevel: the amount of output (integer from 0 to 5)
cdensity: if TRUE, densfn replaced by pre-compiled C version
cfunctional: if TRUE, functional replaced by pre-compiled C version
cfns: used to jointly set both cdensity and cfunctional flags
vectorlength: if a positive integer, specifies that arguments to
the R functions densfn and functional should be given as
a single vector of length vectorlength, rather than as a
sequence of separate real numbers (when vectorlength=0)
(vectorlength is ignored if cdensity=TRUE)
write: if TRUE, create separate files "amcmcvals" and "amcmcsigmas"
(which can then be sourced into R), saving the chain values
and log proposal variances after each batch
Details:
Must first type "source('amcmc')" to load the package. (This will
automatically compile the C code, unless a file "amcmc.so" already
exists; compiling can always be forced by typing "cfns()".)
densfn and functional must each take the same number of arguments
(excluding arguments with pre-defined default values); that common
number of arguments, called K, is the algorithm's dimension.
Certain specific densities and functionals are pre-defined, including
baseballlogdens and baseballfirstfunct for the baseball data of Efron
and Morris (1975) and Morris (1983), as analysed in Rosenthal (1996).
The package works by calling the C program from within R; the C
program in turn normally uses R to evaluate the function calls.
If cdensity and/or cfunctional (or cfns) are TRUE, the C program
instead evaluates the corresponding function internally (much
faster). This requires first creating an auxiliary C file NAME.c
which defines mydim, mydensity, and myfunctional as C objects,
and then typing "cfns('NAME')" from within R. (For examples,
see the included files "simple.c" and "baseball.c".)
For handling such auxiliary C files, the packages provides two
additional R functions in addition to cfns(). The function
editcfns('NAME') invokes the editor getOption('editor') on the
file NAME.c (beginning with a dummy version if the file does
not already exist), and then automatically runs cfns('NAME').
[To change editors, type e.g. options(editor="vi").] And, showcfns()
shows the base name of the auxiliary C file last compiled.
Value:
Returns the estimate of the expected value of functional, with
respect to the density function densfn.
References:
B. Efron and C. Morris (1975), Data analysis using Stein's estimator
and its generalizations. J. Amer. Stat. Assoc. 70, 311-319.
C. Morris (1983), Parametric empirical Bayes confidence intervals.
Scientific Inference, Data Analysis, and Robustness, 25-50.
G.O. Roberts and J.S. Rosenthal (2006), Examples of Adaptive MCMC.
J. Comp. Graph. Stat., to appear.
J.S. Rosenthal (1996), Analysis of the Gibbs Sampler for a model
related to James-Stein Estimators. Stat. and Comput. 6, 269-275.
Examples:
# load the amcmc package
source('amcmc')
# estimate E[X^2] where X ~ N(0,1)
amcmc()
# do the same, more accurately and more verbosely
amcmc(numbatches=5000, batchlength=100, verboselevel=2)
# do the same, but much faster (after downloading "simple.c")
cfns("simple")
amcmc(cfns=TRUE)
# do the same, more accurately
cfns("simple")
amcmc(numbatches=5000, batchlength=100, initval=0, cfns=TRUE)
# estimate E[X^4] where X ~ Uniform[0,1]
myfn <- function(x) return(x^4)
amcmc(dunif, myfn)
# Note: with old versions of R, may need quotes: amcmc("dunif", "myfn")
# estimate E[2(X-Y+Z)^2] where
# X ~ N(0,1), Y ~ Uniform[0,1], Z ~ Exp(1) (all indepedent)
myfn <- function(x,y,z) return(2*(x-y+z)^2)
mydens <- function(x,y,z) return(dnorm(x)*dunif(y)*dexp(z))
amcmc(mydens, myfn)
# do the same, but using R functions in vector form
myfnvect <- function(theargs)
return(2*(theargs[1]-theargs[2]+theargs[3])^2)
mydensvect <- function(theargs)
return(dnorm(theargs[1])*dunif(theargs[2])*dexp(theargs[3]))
amcmc(mydensvect, myfnvect, vectorlength=3)
# estimate the mean of the first theta value in the baseball model
amcmc(baseballlogdens, baseballfirstfunct, logflag=TRUE, verboselevel=3)
# do the same, but much faster (after downloading "baseball.c")
cfns("baseball")
amcmc(logflag=TRUE, cfns=TRUE, verboselevel=3)
# do the same, but using R functions in vector form
amcmc(baseballlogdensvect, baseballfirstfunctvect, logflag=TRUE,
verboselevel=3, vectorlength=20)
# do the same, writing output to files
amcmc(baseballlogdens, baseballfirstfunct, logflag=TRUE,
verboselevel=3, write=TRUE)
source('amcmcvals');
plot(amcmcvals[1,], type='l'); # plot the theta1 values
source('amcmcsigmas');
plot(amcmcsigmas[3,], type='l'); # plot the log(sigma3) values
Entire package and instructions available for free download at:
http://probability.ca/amcmc/
Author:
Jeffrey S. Rosenthal, University of Toronto (probability.ca/jeff)