# file "Rcompwise" -- Componentwise Metropolis algorithm, systematic scan

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 = 2100  # run length (number of scans)
B = 100  # amount of burn-in
X = c(runif(1,0,5),runif(1,0,4))  # overdispersed starting distribution
sigma = 1  # proposal scaling
x1list = x2list = hlist = rep(0,2*M)

for (i in 1:M) {
  for (coord in 1:2) {
    Y = X
    Y[coord] = X[coord] + sigma * rnorm(1)  # propose in direction "coord"
    U = runif(1)  # for accept/reject
    alpha = g(Y) / g(X)  # for accept/reject
    if (U < alpha)
	X = Y  # accept proposal
    x1list[2*i-2+coord] = X[1]
    x2list[2*i-2+coord] = X[2]
    hlist[2*i-2+coord] = h(X)
  }
}

cat("mean of h is about", mean(hlist[(2*B+1):(2*M)]), "\n")

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

se1 =  sd(hlist[(2*B+1):(2*M)]) / sqrt(2*(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[(2*B+1):(2*M)]) ), "\n")