# file "Rnorm" -- a simple geometric Metropolis algorithm in one dimension g = function(y) { dnorm(y, 5, 4) } h = function(y) { return(y^2) } # mean should be 5^2 + 4^2 = 41 M = 21000 # run length B = 1000 # amount of burn-in X = runif(1,0,20) # overdispersed starting distribution X = 5 sigma = 5 # 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[(B+1):M]) ), "\n") plot(xlist, type='l', ylim=c(-5,15)) # acf(xlist) # acf(hlist)