# file "Rdisc2" -- a discrete Metropolis algorithm,
#  			to estimate E[ 2^Sqrt[ Sum[Y] ] ]
#			where Y \in {0,1}^{10 x 10}
#			and Pr[Y] \propto 2^(-Sum[neighbour pairs])

g = function(y) {
    thesum = 0
    for (i in 1:10)
      for (j in 1:9)
        thesum = thesum + y[i,j]*y[i,j+1]
    for (i in 1:9)
      for (j in 1:10)
        thesum = thesum + y[i,j]*y[i+1,j]
    return(2^(-thesum))
}

h = function(y) {
    return(2^sqrt(sum(y)))
}

M = 11000  # run length
B = 1000  # amount of burn-in
X = matrix( sample(c(0,1),100,replace=TRUE), nrow=10 )
xlist = rep(0,M)  # for keeping track of chain values at [1,1]
hlist = rep(0,M)  # for keeping track of h function values
numaccept = 0;

for (i in 1:M) {
    Y = X
    j = sample(10,1)
    k = sample(10,1)
    Y[j,k] = 1 - Y[j,k]
    alpha = g(Y) / g(X)  # for accept/reject
    if (runif(1) < alpha) {
	X = Y  # accept proposal
        numaccept = numaccept + 1;
    }
    xlist[i] = X[1,1];
    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(hlist, type='l')
# acf(xlist)
# plot(xlist, type='l')
# acf(xlist)