# file "Rnorm2" -- a simple geometric Metropolis algorithm in one dimension, # with confidence intervals g = function(y) { dnorm(y, 5, 4) } h = function(y) { return(y^2) } # mean should be 5^2 + 4^2 = 41 M = 11000 # run length B = 1000 # amount of burn-in X = runif(1,0,20) # 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')