# file "Ropt" -- optimising a high-dimensional RWM algorithm

# specify the target covariance and mean
dim = 5
tmpmat = matrix(
c( 0.458731768, -1.138985820,  0.784868032, -1.017293716, -0.723827196,
 0.564371814, -0.155166457,  0.902601504,  0.186805644,  0.952120668,
 1.486382516,  1.458366992,  0.709626067, -1.138458393,  0.897889709,
 0.431802958, -0.721054075, -0.656762999,  1.940817630, -1.015356510,
 0.014969681,  0.422744911,  0.073882617, -0.007217515,  0.664031927),
								nrow=dim )
targSigma = t(tmpmat) %*% tmpmat
targMean = c( 0.2133594,  0.5641370, -0.8031936, -1.1818973, -0.7307054 )
identity = diag(dim)

# specify the proposal increment covariance (choose one)
propSigma = identity
# propSigma = 0.1 * identity
# propSigma = targSigma
# propSigma = 1.4 * targSigma

library(mvtnorm)

logpi = function(X) { dmvnorm(X, targMean, targSigma, log=TRUE) }

h = function(X) { return(X[1]^2) }

M = 4000  # run length
B = 1000  # amount of burn-in
X = runif(dim,-5,5)  # overdispersed starting distribution
x1list = rep(0,M)  # for keeping track of chain values
x2list = rep(0,M)  # for keeping track of chain values
hlist = rep(0,M)  # for keeping track of functional values
numaccept = 0;

for (i in 1:M) {
    Y = X + rmvnorm(1, sigma = propSigma)  # proposal value
    U = runif(1)  # for accept/reject
    # Use log scale for accept/reject to avoid zeroes.
    if (log(U) < logpi(Y) - logpi(X)) {
	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");
estmean = mean(hlist[(B+1):M]);
cat("estimated mean of h is ", estmean, "\n")
cat("true mean of h is ", targMean[1]^2 + targSigma[1,1], "\n");

varfact <- function(xxx) { 2 * sum(acf(xxx, plot=FALSE, lag.max=250)$acf) - 1 }
thevarfact = varfact(hlist[(B+1):M])
cat("varfact =", thevarfact, "\n")

se1 =  sd(hlist[(B+1):M]) / sqrt(M-B)
cat("iid standard error would be about", se1, "\n")

se2 = se1 * sqrt( thevarfact );
cat("true standard error is about", se2, "\n")

cat("95% confidence interval is about (", 
      estmean - 1.96*se2, ",", estmean + 1.96*se2, ")\n")

plot(x1list, type='l')
# plot(x1list[(B+1):M], type='l')
# acf(x1list[(B+1):M])
# acf(x1list[(B+1):M], lag.max=400)
# acf(hlist[(B+1):M], lag.max=250)