# file "Rmwg" -- Metropolis-within-Gibbs algorithm, systematic scan g = function(x) { if ( (x[1]<0) || (x[1]>5) || (x[2]<0) || (x[2]>4) ) return(0) else return( abs( cos(sqrt(x[1]*x[2])) ) ) } h = function(x) { return( exp(x[1]) + x[2]^2 ) } M = 2100 # run length (number of scans) B = 100 # amount of burn-in X = c(runif(1,0,5),runif(1,0,4)) # overdispersed starting distribution sigma = 1 # proposal scaling x1list = x2list = hlist = rep(0,2*M) for (i in 1:M) { for (coord in 1:2) { Y = X Y[coord] = X[coord] + sigma * rnorm(1) # propose in direction "coord" U = runif(1) # for accept/reject alpha = g(Y) / g(X) # for accept/reject if (U < alpha) X = Y # accept proposal x1list[2*i-2+coord] = X[1] x2list[2*i-2+coord] = X[2] hlist[2*i-2+coord] = h(X) } } cat("mean of h is about", mean(hlist[(2*(B+1)):(2*M)]), "\n") plot(x1list, x2list, type='l') # plot(hlist, type='l') # plot(x1list, type='l') # plot(x2list, type='l') se1 = sd(hlist[(2*(B+1)):(2*M)]) / sqrt(2*(M-B)) cat("iid standard error would be about", se1, "\n") varfact <- function(xxx) { 2 * sum(acf(xxx, plot=FALSE)$acf) - 1 } cat("true standard error is about", se1 * sqrt( varfact(hlist[(B+1):M]) ), "\n")