# file "Rheavy" -- a non-geometric Metropolis algorithm in one dimension g = function(y) { 1 / (1+y^4) } h = function(y) { y^2 } M = 21000 # run length B = 1000 # amount of burn-in X = runif(1,-10,10) # overdispersed starting distribution sigma = 1 # proposal scaling xlist = rep(0,M) # for keeping track of chain values hlist = rep(0,M) # for keeping track of h function values numaccept = 0; for (i in 1:M) { Y = X + sigma * rnorm(1) # proposal value U = runif(1) # for accept/reject alpha = g(Y) / g(X) # for accept/reject if (U < alpha) { X = Y # accept proposal numaccept = numaccept + 1; } xlist[i] = X; hlist[i] = h(X); } cat("ran Metropolis algorithm for", M, "iterations, with burn-in", B, "\n"); cat("acceptance rate =", numaccept/M, "\n"); estmean = mean(hlist[(B+1):M]) cat("mean of h is about", estmean, "\n") se1 = sd(hlist[(B+1):M]) / sqrt(M-B) cat("iid standard error would be about", se1, "\n") varfact <- function(xxx) { 2 * sum(acf(xxx, plot=FALSE)$acf) - 1 } se2 = se1 * sqrt( varfact(hlist[(B+1):M]) ) cat("true standard error is about", se2, "\n") cat("95% confidence interval is about (", estmean - 1.96*se2, ",", estmean + 1.96*se2, ")\n") plot(xlist, type='l') # acf(xlist) # acf(hlist)