Computation for the Introduction to MCMC Chapter of Handbook of Markov Chain Monte Carlo

This technical report does the computation for the“Introduction to MCMC”chapter of Brooks, Gelman, Jones and Meng (forthcoming). All analyses are done in R (R Development Core Team, 2008) using the Sweave function so this entire technical report and all of the analyses reported in it are exactly reproducible by anyone who has R with the mcmc package (Geyer, 2005) installed and the R noweb file specifying the document. 1 R Package MCMC We use R statistical computing environment (R Development Core Team, 2008) in our analysis. It is free software and can be obtained from http://cran.r-project.org. Precompiled binaries are available for Windows, Macintosh, and popular Linux distributions. We use the contributed package mcmc. If R has been installed, but this package has not yet been installed, do install.packages("mcmc") from the R command line (or do the equivalent using the GUI menus if on Apple Macintosh or Microsoft Windows). This may require root or administrator privileges. Assuming the mcmc package has been installed, we load it > library(mcmc) The version of R used to make this document is 2.11.1. The version of the mcmc package used to make this document is 0.8. This entire document and all of the calculations shown were made using the R command Sweave and hence are exactly reproducible by anyone who has R and the R noweb (RNW) file from which it was created. Both the RNW file and and the PDF document produced from it are available at http://www.stat.umn.edu/geyer/mcmc. For further details on the use of Sweave and R see the web site http://www/stat.umn.edu/~charlie/Sweave/. Not only can one exactly reproduce the results in the printable document, one can also modify the parameters of the simulation and get different results. Finally, we set the seeds of the random number generator so that we obtain the same results every time. To get different results, obtain the RNW file, change this statement, and reprocess using Sweave and LATEX. > set.seed(42) 2 AR(1) Simulation Simulate an AR(1) time series. > rho <0.99 > nsim <10000 > x <0 > for (i in 2:nsim) x out sqrt((1 + rho)/(1 rho)/nsim) [1] 0.1410674 In order to get a 10% relative error, we would need sample size > ceiling((1 + rho)/(1 rho)/0.1^2) [1] 19900 In order to get a 10% relative error, we would need sample size > ceiling((1 + rho)/(1 rho)/0.01^2) [1] 1990000 Variance inflation factor (VIF) > (1 + rho)/(1 rho) [1] 199 Effective sample size > nsim/((1 + rho)/(1 rho)) [1] 50.25126 Figure 1 (page 3) is made by the following R statements > par(mar = c(5, 4, 1, 1) + 0.1) > plot(out, ylab = "x") Figure 2 (page 4) is made by the following R statements > foo par(mar = c(5, 4, 1, 1) + 0.1) > plot(as.ts(foo), ylab = expression(hat(mu)[n]), xlab = expression(n)) Figure 3 (page 5) is made by the following R statements > lag.max <500 > par(mar = c(5, 4, 1, 1) + 0.1) > plot(acf(out, lag.max = lag.max), ci.col = "black", + main = "") > curve(rho^x, from = 0, to = lag.max, add = TRUE, + lty = 3, lwd = 2)