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

Sw ea v e Do cumen tation for “ Implemen ting Mark o v c hain Mon te Carlo: Estimating with confidence ” James M. Flegal and Galin L. Jones Ma y 18, 2022 1 In tro duction This file is the Sw eav e do cumen tation for the examples provided in Flegal and Jones (2010). 1.1 Batc h Means The follo wing function is required to calculate the v ariance estimators usin g batc h means (BM) and ov erlapping batch means (OLBM) . W e are thankful to Murali Haran whom wrote the original function to implemen t BM which we ha ve expanded here. > id <- function(x) return(x) > mcse <- function(vals, bs = "sqr oot", g = id, meth = "BM", warn = FALSE) { + N <- length(vals) + if (N < 1000) { + if (warn) + cat("WARNING: too few samples ( less than 1000)\n") + if (N < 10) + return(NA) + } + if (bs == "sqroot") { + b <- floor(sqrt(N)) + a <- floor(N/b) + } + else if (bs == "cuberoot") { + b <- floor(N^(1/3)) + a <- floor(N/b) + } + else { + stopifnot(is.numeric(bs)) 1 + b <- floor(bs) + if (b > 1) + a <- floor(N/b) + else stop("batch size invalid ( bs=", bs, ")") + } + if (meth == "BM") { + Ys <- sapply(1:a, function(k) r eturn(mean(g(vals[((k - + 1) * b + 1):(k * b)])))) + muhat <- mean(Ys) + sigmahatsq <- b * sum((Ys - muhat) ^2)/(a - 1) + bmse <- sqrt(sigmahatsq/N) + return(bmse) + } + if (meth == "OBM") { + a <- N - b + 1 + Ys <- sapply(1:a, function(k) r eturn(mean(g(vals[k:(k + + b - 1)])))) + muhat <- mean(Ys) + sigmahatsq <- N * b * sum((Ys - muhat )^2)/(a - 1)/a + bmse <- sqrt(sigmahatsq/N) + return(bmse) + } + else { + stop("method specified inval id (meth=", meth, ")") + } + } 2 Normal AR(1) Marko v Chains Consider the normal AR(1) time series defined by X n +1 = ρX n +  n (1) where the  n are i.i.d. N(0,1) and ρ < 1. This Marko v chain has in v arian t distribution N  0 , 1 / (1 − ρ 2 )  . In our example, we p erfomred calulcations for ρ ∈ { 0 . 5 , 0 . 95 } . 2.1 Mark o v Chain Sampler The follo wing ch unk of co de giv es general functions needed to sample from (1). Here w e hav e a function that provides an observ ation from the chain, an observ ation q steps ahead with a defualt of one step ahead, and p observ ations from the chain. > ar1 <- function(m, rho, tau) { + rho * m + rnorm(1, 0, tau) 2 + } > ar1.q <- function(m, rho, tau, q = 1) { + for (i in 1:q) { + m <- rho * m + rnorm(1, 0, tau) + } + m + } > ar1.gen <- function(mc, p, rho, tau, q = 1) { + loc <- length(mc) + junk <- double(p) + mc <- append(mc, junk) + for (i in 1:p) { + j <- i + loc - 1 + mc[(j + 1)] <- ar1(mc[j], rho, tau) + } + return(mc) + } 2.2 Additional F unctions The following are additional functions necessary for later calculations. The first calculates the estimated first and third quartile while the second calculates the asso ciated MCSE via subsampling. Notice this function is similar to the function necessary for OLBM. > quant <- function(input) { + quantile(input, prob = c(0.25, 0.75), type = 1) + } > subsampling <- function(vals) { + N <- length(vals) + b <- floor(sqrt(N)) + a <- N - b + 1 + Ys <- sapply(1:a, function(k) return(quant(vals[k:(k + b - + 1)]))) + muhat <- apply(Ys, 1, mean) + sigmahatsq <- N * b * apply((Ys - muhat)^2, 1, sum)/(a - + 1)/a + bmse <- sqrt(sigmahatsq/N) + return(bmse) + } 2.3 Sim ulation Settings and Calculations In this next ch unk of co de, we first give the sim ulation settings used throughout the pap er. W e then generate the t wo Mark ov chains for ρ ∈ { 0 . 5 , 0 . 95 } and calculate the corresp onding estimates and MCSEs. 3 > 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[1:k], meth = "BM"))) > obm.est1 <- sapply(1:n, function(k) return(mcse(chain1[1:k], + meth = "OBM"))) > quartile1 <- sapply(1:n, function(k) return(quant(chain1[1:k]) )) > q.mcse1 <- sapply(1:n, function(k) return(subsampling(chain1[1 :k]))) > chain2 <- ar1.gen(1, (n - 1), rho2, 1) > mean2 <- cumsum(chain2)/seq(along = chain2) > bm.est2 <- sapply(1:n, function(k) return(mcse(chain2[1:k], meth = "BM"))) > obm.est2 <- sapply(1:n, function(k) return(mcse(chain2[1:k], + meth = "OBM"))) > quartile2 <- sapply(1:n, function(k) return(quant(chain2[1:k]) )) > q.mcse2 <- sapply(1:n, function(k) return(subsampling(chain2[1 :k]))) 2.4 Initial Examination of Output The follow ing c hunk of code will create plots for the initial examination of output. The plots are also rep eated in the do cumen t. > 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 = "", + ylab = "", lwd = 2, xlim = c(0, n)) > abline(h = 0) > par(mfrow = c(1, 1), mar = c(5, 4, 4, 2)) 4 Time−Series vs. Iteration 0 500 1000 1500 2000 −4 −2 0 2 4 0 5 10 15 20 25 30 0.0 0.4 0.8 Autocorrelation vs. La g Running A verage vs. Iteration 0 500 1000 1500 2000 −0.5 0.0 0.5 1.0 > rho = rho2 > chain <- chain2 > mean <- mean2 > 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 = "", + ylab = "", lwd = 2, xlim = c(0, n)) > abline(h = 0) > par(mfrow = c(1, 1), mar = c(5, 4, 4, 2)) 5 Time−Series vs. Iteration 0 500 1000 1500 2000 −10 0 5 0 5 10 15 20 25 30 0.0 0.4 0.8 Autocorrelation vs. La g Running A verage vs. Iteration 0 500 1000 1500 2000 0 1 2 3 4 5 2.5 Running MCSEs for Exp ectations The followi ng c hunk of co de creates the plot of the running MCSEs and run- ning estimates for the exp ectations with confidence bounds. Again, the plot is con tained in the do cumen t. > rho = rho1 > chain <- chain1 > mean <- mean1 > u.obm <- mean + crit.obm * obm.est1 > l.obm <- mean - crit.obm * obm.est1 > u.bm <- mean + crit.bm * bm.est1 > l.bm <- mean - crit.bm * bm.est1 > ts.plot(mean, main = "Running Average", xlab = "Iteration", ylab = "", + ylim = c(min(l.obm[10:n]), max(u.obm[10:n])), lwd = 2, xlim = c(0, + n)) > abline(h = 0) > points(iter, u.obm, type = "l", lty = 4, lwd = 2) > points(iter, l.obm, type = "l", lty = 4, lwd = 2) 6 Running A verage Iteration 0 500 1000 1500 2000 −1.0 −0.5 0.0 0.5 1.0 > rho = rho2 > chain <- chain2 > mean <- mean2 > u.obm <- mean + crit.obm * obm.est2 > l.obm <- mean - crit.obm * obm.est2 > u.bm <- mean + crit.bm * bm.est2 > l.bm <- mean - crit.bm * bm.est2 > ts.plot(mean, main = "Running Average", xlab = "Iteration", ylab = "", + ylim = c(min(l.obm[10:n]), max(u.obm[10:n])), lwd = 2, xlim = c(0, + n)) > abline(h = 0) > points(iter, u.obm, type = "l", lty = 4, lwd = 2) > points(iter, l.obm, type = "l", lty = 4, lwd = 2) 7 Running A verage Iteration 0 500 1000 1500 2000 −1 0 1 2 3 4 5 6 2.6 Running Quartile Plots Here is the co de for the running quartile plots with and without confidence b ounds for ρ ∈ { 0 . 5 , 0 . 95 } . There are a total of four ch unks 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", + ylab = "", lwd = 2, ylim = c(max(u.sub[, 100:n]), min(l.sub[, + 100:n])), xlim = c(0, n)) > abline(h = qnorm(0.25, 0, sqrt(1/(1 - rho^2)))) > abline(h = qnorm(0.75, 0, sqrt(1/(1 - rho^2)))) 8 Running Quartiles Iteration 0 500 1000 1500 2000 1.0 0.5 0.0 −0.5 −1.0 > ts.plot(t(quartiles), main = "Running Quartiles", xlab = "Iteration", + ylab = "", lwd = 2, ylim = c(max(u.sub[, 100:n]), min(l.sub[, + 100:n])), xlim = c(0, n)) > points(iter, t(u.sub[1, ]), type = "l", lty = 4, lwd = 2) > points(iter, t(u.sub[2, ]), type = "l", lty = 4, lwd = 2) > points(iter, t(l.sub[1, ]), type = "l", lty = 4, lwd = 2) > points(iter, t(l.sub[2, ]), type = "l", lty = 4, lwd = 2) > abline(h = qnorm(0.25, 0, sqrt(1/(1 - rho^2)))) > abline(h = qnorm(0.75, 0, sqrt(1/(1 - rho^2)))) 9 Running Quartiles Iteration 0 500 1000 1500 2000 1.0 0.5 0.0 −0.5 −1.0 > rho = rho2 > chain <- chain2 > quartiles <- quartile2 > u.sub <- quartiles + rbind(crit.obm, crit.obm) * q.mcse2 > l.sub <- quartiles - rbind(crit.obm, crit.obm) * q.mcse2 > ts.plot(t(quartiles), main = "Running Quartiles", xlab = "Iteration", + ylab = "", lwd = 2, ylim = c(max(u.sub[, 100:n]), min(l.sub[, + 100:n])), xlim = c(0, n)) > abline(h = qnorm(0.25, 0, sqrt(1/(1 - rho^2)))) > abline(h = qnorm(0.75, 0, sqrt(1/(1 - rho^2)))) 10 Running Quartiles Iteration 0 500 1000 1500 2000 8 6 4 2 0 −2 −4 > ts.plot(t(quartiles), main = "Running Quartiles", xlab = "Iteration", + ylab = "", lwd = 2, ylim = c(max(u.sub[, 100:n]), min(l.sub[, + 100:n])), xlim = c(0, n)) > points(iter, t(u.sub[1, ]), type = "l", lty = 4, lwd = 2) > points(iter, t(u.sub[2, ]), type = "l", lty = 4, lwd = 2) > points(iter, t(l.sub[1, ]), type = "l", lty = 4, lwd = 2) > points(iter, t(l.sub[2, ]), type = "l", lty = 4, lwd = 2) > abline(h = qnorm(0.25, 0, sqrt(1/(1 - rho^2)))) > abline(h = qnorm(0.75, 0, sqrt(1/(1 - rho^2)))) 11 Running Quartiles Iteration 0 500 1000 1500 2000 8 6 4 2 0 −2 −4 2.7 Calculations for Exp ectations The followi ng are calculations rep orted in the pap er based on n = 2000 iterations in the chain. > signif(mean1[n], 6) [1] -0.0336436 > signif(crit.obm[n] * obm.est1[n], 6) [1] 0.0556022 > signif(mean2[n], 6) [1] -0.506755 > signif(crit.obm[n] * obm.est2[n], 6) [1] 0.450564 > signif(obm.est2[n], 6) [1] 0.351458 12 Here the implemen tation of fixed width pro cedures using ρ = 0 . 95. Here we ha ve added 198000 iterations to the chain for a total of 200000 to sp ead up the pro cessing time. If one desired a lev el of precision smaller than 0.1 it w ould probably b e necessary to add additional iterations to the chain. > chain2 <- ar1.gen(chain2, 198000, rho2, 1) > chain <- chain2 > half <- crit.obm[n] * obm.est2[n] > N <- n > while (half + 1/N > 0.1) { + N <- N + 1000 + b <- floor(sqrt(N)) + t.OBM <- qt(0.9, (N - b + 1)) + est.OBM <- mcse(chain[1:N], meth = "OBM") + half <- t.OBM * est.OBM + } > N [1] 60000 > half [1] 0.0996075 > mean(chain[1:N]) [1] -0.04415015 > signif(mean(chain[1:N]), 6) [1] -0.0441502 2.8 Calculations for Quartiles Here are the same calculations for the quartiles. First the caculations based on n = 2000 iterations for the quartiles. > signif(quartile1[, n], 6) 25% 75% -0.816573 0.777758 > signif(crit.obm[n] * q.mcse1[, n], 6) 25% 75% 0.0685046 0.0648736 > signif(quartile2[, n], 6) 13 25% 75% -2.73625 1.77970 > signif(crit.obm[n] * q.mcse2[, n], 6) 25% 75% 0.480563 0.466388 Then the fixed width calculations for quan tiles with ρ = 0 . 95. First for individual CIs, the with a Bonferonni correction. These calculations are not included in the pap er and are commented out since the computational time is appro ximately an hour. The co de is included for future reference. chain <- chain2 half <- max(crit.obm[n]*q.mc se2[,n]) N <- n while(half + 1 / N > .1){ N <- N + 2000 b <- floor(sqrt(N)) # batch size t.sub <- qt(.9, (N - b + 1)) est.sub <- subsampling(chain[1:N]) half <- max(t.sub*est.sub) } N half quant(chain[1:N]) t.sub*est.sub ##### # Then for Bonferonni Correction ##### t.sub <- t.sub <- qt(.9875, (N - b + 1)) half <- max(t.sub*est.sub) while(half + 1 / N > .1){ N <- N + 2000 b <- floor(sqrt(N)) # batch size t.sub <- qt(.9875, (N - b + 1)) est.sub <- subsampling(chain[1:N]) half <- max(t.sub*est.sub) } N half quant(chain[1:N]) t.sub*est.sub 14 3 T-Distribution Example Suppose our goal is to estimate the first t wo moments of a Students t distribution with 4 degrees of freedom and having densit y m ( x ) = 3 8  1 + x 2 4  − 5 / 2 Ob viously , there is nothing ab out this that requires MCMC since we can easily calculate that E m X = 0 and E m X 2 = 2. Nev ertheless, w e will use a data augmen tation algorithm based on the joint densit y π ( x, y ) = 4 √ 2 π y 3 2 e − y ( 2+ x 2 / 2 ) 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 lo oks lik e ( x 0 , y 0 ) → ( x, y 0 ) → ( x, y ). 3.1 Mark o v Chain Sampler The following function samples from p iterations in the Marko v c hain starting from the v alue (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:p) { + j <- i + loc - 1 + x[(j + 1)] <- rnorm(1, 0, sqrt(1/y[j])) + y[(j + 1)] <- rgamma(1, 5/2, rate = (2 + x[(j + 1)]^2/2)) + } + return(cbind(x, y)) + } 3.2 Calculaitons and Plots Using the sampler, we can run the sim ulation and create the necessary plots. > n <- 2000 > iter <- seq(1, n) > set.seed(100) > sample <- t.gen((n - 1)) > x <- sample[, 1] > y <- sample[, 2] > RB <- 1/y 15 > x.mean <- cumsum(x)/seq(along = x) > x2.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) First Moment Iteration 0 500 1000 1500 2000 0.0 0.2 0.4 0.6 0.8 1.0 Std RB > ts.plot(x2.mean, main = "Second Moment", xlab = "Iteration", + ylab = "", lwd = 2, xlim = c(0, n)) > points(iter, RB.mean, type = "l", lty = 4, lwd = 2) > abline(h = 2) > legend("bottomright", c("Std", "RB"), lty = c(1, 4), lwd = 2) 16 Second Moment Iteration 0 500 1000 1500 2000 0.5 1.0 1.5 2.0 Std RB 3.3 In terv al Estimates The follo wing c hunk of code will calculate in terv al estimates using the same data from ab o ve. > obm.x <- sapply(1:n, function(k) return(mcse(x[1:k], meth = "OBM"))) > obm.RB <- sapply(1:n, function(k) return(mcse(RB[1:k], meth = "OBM"))) > obm.x2 <- sapply(1:n, function(k) return(mcse(x[1:k]^2, meth = "OBM"))) Then we can plot the estimates for the first momen t using standard methods. Notice, there is no uncertaint y 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("topright", c("Std", "RB"), lty = c(1, 4), lwd = 2) > points(iter, upper, type = "l", lty = 1, lwd = 1) > points(iter, lower, type = "l", lty = 1, lwd = 1) 17 First Moment Iteration 0 500 1000 1500 2000 0.0 0.2 0.4 0.6 0.8 1.0 Std RB No w the plot for second moment for b oth the usual and RB estimator. > upper <- RB.mean + crit.obm * obm.RB > lower <- RB.mean - crit.obm * obm.RB > upper.x2 <- x2.mean + crit.obm * obm.x2 > lower.x2 <- x2.mean - crit.obm * obm.x2 > ts.plot(x2.mean, main = "Second Moment", xlab = "Iteration", + ylab = "", lwd = 2, xlim = c(0, n), ylim = c(min(lower.x2[10:n]), + max(upper.x2[10:n]))) > points(iter, RB.mean, type = "l", lty = 4, lwd = 2) > abline(h = 2) > legend("bottomright", c("Std", "RB"), lty = c(1, 4), lwd = 2) > points(iter, upper, type = "l", lty = 4, lwd = 1) > points(iter, lower, type = "l", lty = 4, lwd = 1) > points(iter, upper.x2, type = "l", lty = 1, lwd = 1) > points(iter, lower.x2, type = "l", lty = 1, lwd = 1) 18 Second Moment Iteration 0 500 1000 1500 2000 0.5 1.0 1.5 2.0 2.5 Std RB 4 Estimating Marginals Example Suppose Y i | µ, θ ∼ N( µ, θ ) independen tly for i = 1 , . . . , m where m ≥ 3 and prior ν ( µ, θ ) ∝ θ − 1 / 2 . The target is the p osterior density π ( µ, θ | y ) ∝ θ − ( m +1) / 2 e − m 2 θ ( s 2 +( ¯ y − µ ) 2 ) where s 2 is the usual biased sample v ariance. It is easy to see that µ | θ , y ∼ N( ¯ y , θ/m ) and that θ | µ, y ∼ IG(( m − 1) / 2 , m [ s 2 + ( ¯ y − µ ) 2 ] / 2). W e will conisder the Gibbs sampler that up dates µ then θ so that a one-step transition is given b y ( µ 0 , θ 0 ) → ( µ, θ 0 ) → ( µ, θ ). W e will use this Gibbs sampler to estimate the marginal densities of µ and θ . 4.1 Mark o v Chain Sampler The first function provides an observ ation one step ahead using the Gibbs sam- pler. The second function results in p observ ations from the Marko v chain. > ex2.gibbs <- function(m, t, y.bar = 1, s2 = 4, n = 11) { + alpha <- (n - 1)/2 + beta <- n * (s2 + (y.bar - m)^2)/2 + t <- rgamma(1, alpha, beta) 19 + t <- 1/t + m <- rnorm(1, y.bar, sqrt(t/n)) + cbind(m, t) + } > ex2.gen <- function(mc, p = 100, q = 1, y.bar = 1, s2 = 4, n = 11) { + if (is.matrix(mc) == TRUE) { + loc <- dim(mc)[1] + } + else { + loc <- 1 + mc <- t(as.matrix(mc)) + } + junk <- matrix(rep(NA, p * 2), ncol = 2) + mc <- rbind(mc, junk) + for (i in 1:p) { + j <- i + loc - 1 + mc[(j + 1), ] <- ex2.gibbs(mc[j, 1], mc[j, 2]) + } + return(mc) + } 4.2 Calculations and Settings No w supp ose m = 11, ¯ y = 1 and s 2 = 4. W e simulated 2000 realizations of the Gibbs sampler starting from ( µ 1 , λ 1 ) = (1 , 1). The marginal density plots were created using the default settings for the density function in R. The biv ariate densit y plot w as created using R functions kde2d and persp . The resulting p osterior is simple, so it is no surprise that the Gibbs sampler has b een sho wn to conv erge in just a few iterations (Jones and Hob ert, 2001). > n <- 2000 > set.seed(100) > start <- c(1, 1) > sample <- ex2.gen(start, (n - 1)) > mu <- sample[, 1] > theta <- sample[, 2] No w we can plot the marginal densities and corresp onding histograms for µ and θ . > hist(mu, main = "", xlab = "", freq = F) > junk <- density(mu) > points(junk$x, junk$y, type = "l") 20 Density −3 −2 −1 0 1 2 3 4 0.0 0.1 0.2 0.3 0.4 0.5 > hist(theta, main = "", xlab = "", freq = F, ylim = c(0, 0.18)) > junk <- density(theta) > points(junk$x, junk$y, type = "l") 21 Density 0 5 10 15 20 25 30 35 0.00 0.05 0.10 0.15 No w we can plot the biv ariate densit y plot. > library(MASS) > junk <- kde2d(mu, theta, n = 50, lims = c(-1.5, 3.5, 1, 15)) > persp(junk, theta = 120, phi = 30, xlab = "mu", ylab = "theta", + zlab = "", expand = 0.6) 22 mu theta 4.3 Alternativ e Curv e Estimation A clever technique for estimating a marginal is based on the same idea as RB estimators (W ei and T anner, 1990). The followi ng ch unk of co de implemen ts this for the current example for the marginal of µ . > hist(mu, main = "", xlab = "", freq = F) > junk <- density(mu) > points(junk$x, junk$y, type = "l") > evals <- seq(-3, 4, 0.01) > norm.evals <- dnorm(evals, 1, sqrt(mean(theta)/11)) > points(evals, norm.evals, type = "l", lty = 4, lwd = 2) 23 Density −3 −2 −1 0 1 2 3 4 0.0 0.1 0.2 0.3 0.4 0.5 References Flegal, J. M. and Jones, G. L. (2010). Implementing Marko v c hain Mon te Carlo: Estimating with confidence. In Bro oks, S., Gelman, A., Jones, G., and Meng, X., editors, Handb o ok of Markov Chain Monte Carlo . Chapman & Hall/CRC Press. Jones, G. L. and Hob ert, J. P . (2001). Honest exploration of int ractable proba- bilit y distributions via Mark ov chain Monte Carlo. Statistic al Scienc e , 16:312– 334. W ei, G. C. G. and T anner, M. A. (1990). A Monte Carlo implementation of the EM algorithm and the p o or man’s data augmentation algorithms. Journal of the Americ an Statistic al Asso ciation , 85:699–704. 24