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)