# file "Rmet2" -- a simple Metropolis algorithm in two dimensions

g = function(x) {
    if ( (x[1]<0) || (x[1]>5) || (x[2]<0) || (x[2]>4) )
	return(0)
    else
	return( abs( cos(sqrt(x[1]*x[2])) ) )
}

h = function(x) { return( exp(x[1]) + x[2]^2 ) }

M = 11000  # run length
B = 1000  # amount of burn-in
X = c(5*runif(1),4*runif(1))  # overdispersed starting distribution (dim=2)
sigma = 1  # proposal scaling
x1list = x2list = hlist = rep(0,M)  # for keeping track of values
numaccept = 0;

for (i in 1:M) {
    Y = X + sigma * rnorm(2)  # proposal value (dim=2)
    U = runif(1)  # for accept/reject
    alpha = g(Y) / g(X)  # for accept/reject
    if (U < alpha) {
	X = Y  # accept proposal
        numaccept = numaccept + 1;
    }
    x1list[i] = X[1];
    x2list[i] = X[2];
    hlist[i] = h(X);
}

cat("ran Metropolis algorithm for", M, "iterations, with burn-in", B, "\n");
cat("acceptance rate =", numaccept/M, "\n");
u = mean(hlist[(B+1):M])
cat("mean of h is about", u, "\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 }
thevarfact = varfact(hlist[(B+1):M]) 
se = se1 * sqrt( thevarfact )
cat("varfact = ", thevarfact, "\n")
cat("true standard error is about", se, "\n")
cat("approximate 95% confidence interval is (", u - 1.96 * se, ",",
                                                u + 1.96 * se, ")\n\n")

# plot(x1list, type='l')
# plot(x2list, type='l')
# plot(x1list, x2list, type='l')
# plot(x1list, x2list, type='p')
# acf(hlist)