# file "RMH" -- a Metropolis-Hastings algorithm

# Proposal distribution: MVN(X_{n-1}, sigma^2 (1+|X_{n-1}|^2) I)

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 ) }

qq = function(x,y) {
    return( 1/(1+sum(x^2))^2 *
		exp( - sum((y-x)^2) / ( 2 * sigma^2 * (1+sum(x^2))^2 ) ) )
}

M = 5000  # run length
B = 100  # amount of burn-in
X = c(2.5,2)  # initial Markov chain value (dim=2 ... central values ...)
sigma = 0.1  # proposal scaling
x1list = x2list = NULL
hlist = NULL  # for keeping track of h function values

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

cat("mean of h is about", mean(hlist), "\n")

# OTHER POSSIBLE COMMANDS:

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

se1 =  sd(hlist) / sqrt(length(hlist))
cat("iid standard error is about", se1, "\n")

# acf(x1list)

varfact <- function(xxx) { 2 * sum(acf(xxx, plot=FALSE)$acf) - 1 }
cat("true standard error is about", se1 * sqrt( varfact(hlist) ), "\n")