Sweave Documentation for 'Implementing Markov chain Monte Carlo: Estimating with confidence'

The following function is required to calculate the variance estimators using batch means (BM) and overlapping batch means (OLBM). We are thankful to Murali Haran whom wrote the original function to implement BM which we have expanded here. Consider the normal AR(1) time series defined by

where the n are i.i.d. N(0,1) and ρ < 1. This Markov chain has invariant distribution N 0, 1/(1 -ρ 2 ) .

In our example, we perfomred calulcations for ρ ∈ {0.5, 0.95}.

The following chunk of code gives general functions needed to sample from (1). Here we have a function that provides an observation from the chain, an observation q steps ahead with a defualt of one step ahead, and p observations from the chain.

The following are additional functions necessary for later calculations. The first calculates the estimated first and third quartile while the second calculates the associated MCSE via subsampling. Notice this function is similar to the function necessary for OLBM.

In this next chunk of code, we first give the simulation settings used throughout the paper. We then generate the two Markov chains for ρ ∈ {0.5, 0.95} and calculate the corresponding estimates and MCSEs.

n <-2000 > iter <-seq(1, n) > crit.bm <-qt(0.9, (sqrt(iter) -1)) > crit.obm <-qt(0.9, (iter -sqrt(iter))) > set.seed(1976) > rho1 <-0.5 > rho2 <-0.95 > chain1 <-ar1.gen(1, (n -1), rho1, 1) > mean1 <-cumsum(chain1)/seq(along = chain1) > bm.est1 <-sapply(1:n, function(k) return(mcse(chain1

The following chunk of code will create plots for the initial examination of output. The plots are also repeated in the document.

rho = rho1 > chain <-chain1 > mean <-mean1 > par(mfrow = c(3, 1), mar = c(3, 4, 4, 2)) > ts.plot(chain, main = “Time-Series vs. Iteration”, xlab = “”, + ylab = “”, xlim = c(0, n)) > abline(h = 2 * sqrt(1/(1 -rho^2))) > abline(h = -2 * sqrt(1/(1 -rho^2))) > acf(chain, main = “Autocorrelation vs. Lag”, ylab = “”, xlab = “”) > ts.plot(mean, main = “Running Average vs. Iteration”, xlab = “”,

acf(chain, main = “Autocorrelation vs. Lag”, ylab = “”, xlab = “”) > ts.plot(mean, main = “Running Average vs. Iteration”, xlab = “”,

The following chunk of code creates the plot of the running MCSEs and running estimates for the expectations with confidence bounds. Again, the plot is contained in the document.

Here is the code for the running quartile plots with and without confidence bounds for ρ ∈ {0.5, 0.95}. There are a total of four chunks of code here, one for each plot.

rho = rho1 > chain <-chain1 > quartiles <-quartile1 > u.sub <-quartiles + rbind(crit.obm, crit.obm) * q.mcse1 > l.sub <-quartiles -rbind(crit.obm, crit.obm) * q.mcse1 > ts.plot(t(quartiles), main = “Running Quartiles”, xlab = “Iteration”,

The following are calculations reported in the paper based on n = 2000 iterations in the chain.

[1] -0.0336436

Here the implementation of fixed width procedures using ρ = 0.95. Here we have added 198000 iterations to the chain for a total of 200000 to spead up the processing time. If one desired a level of precision smaller than 0.1 it would probably be necessary to add additional iterations to the chain.

[1] -0.04415015

[1] -0.0441502

Here are the same calculations for the quartiles. First the caculations based on n = 2000 iterations for the quartiles.

Suppose our goal is to estimate the first two moments of a Students t distribution with 4 degrees of freedom and having density

Obviously, there is nothing about this that requires MCMC since we can easily calculate that E m X = 0 and E m X 2 = 2. Nevertheless, we will use a data augmentation algorithm based on the joint density

so that the full conditionals are X|Y ∼ N(0, y -1 ) and Y |X ∼ Gamma(5/2, 2 + x 2 /2). Consider the Gibbs sampler that updates X then Y so that a one-step transition looks like (x , y ) → (x, y ) → (x, y).

The following function samples from p iterations in the Markov chain starting from the value (1, 1).

t.gen <-function(p, x = 1, y = 1) { + loc <-length(x) + junk <-double(p) +

x <-append(x, junk) + y <-append(y, junk) + for (i in 1

Using the sampler, we can run the simulation and create the necessary plots.

.mean <-cumsum(x^2)/seq(along = x) > RB.mean <-cumsum(1/y)/seq(along = y) > ts.plot(x.mean, main = “First Moment”, xlab = “Iteration”, ylab = “”, + lwd = 2, xlim = c(0, n)) > abline(h = 0) > abline(h = 0, lwd = 2, lty = 4) > legend(“topright”, c(“Std”, “RB”), lty = c(1, 4), lwd = 2)

The following chunk of code will calculate interval estimates using the same data from above. Then we can plot the estimates for the first moment using standard methods. Notice, there is no uncertainty with the RB estimator in this setting.

upper <-x.mean + crit.obm * obm.x > lower <-x.mean -crit.obm * obm.x > ts.plot(x.mean, main = “First Moment”, xlab = “Iteration”, ylab = “”, + lwd = 2, xlim = c(0, n)) > abline(h = 0) > abline(h = 0, lwd = 2, lty = 4) > legend(“t

…(Full text truncated)…