# file "Rmet2" -- a slowly-mixing Metropolis algorithm

g = function(x) { (1/2) * dnorm(x,0) + (1/2) * dnorm(x,10) }

h = function(x) { return( x ) }

M = 1000  # run length
B = 100  # amount of burn-in
X = c(0)  # initial Markov chain value (dim=1)
sigma = 1  # proposal scaling
xlist = NULL
hlist = NULL  # for keeping track of h function values

for (i in 1:M) {
    Y = X + sigma * rnorm(1)  # proposal value (dim=1)
    U = runif(1)  # for accept/reject
    alpha = g(Y) / g(X)  # for accept/reject
    if (U < alpha)
	X = Y  # accept proposal
    xlist = c(xlist, X)
    if (i>B)
	hlist = c(hlist, h(X))  # keep h value if burn-in over
}

cat("mean of h is about", mean(hlist), "\n")
# plot(hlist, type='l')
# plot(xlist, type='l')

se1 =  sd(hlist) / sqrt(length(hlist))
cat("iid standard error is 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")