# 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) )
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')