# file "Rheavy" -- a non-geometric Metropolis algorithm in one dimension g = function(y) { 1 / (1+y^2) / pi } # cauchy distribution h = function(y) { if ( abs(y) < 10 ) return(1) else return(0) } 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"); cat("mean of h is about", mean(hlist[(B+1):M]), "\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 } cat("true standard error is about", se1 * sqrt( varfact(hlist) ), "\n") plot(xlist, type='l')