# file "Ropt" -- optimising a high-dimensional RWM algorithm, in dim=20

# specify the target mean (targMean) and covariance (targSigma)
dim = 20
load("targ20")
identity = diag(dim)

# specify the proposal increment covariance (choose one)
# propSigma = 0.025 * identity
propSigma = (2.38)^2/dim * targSigma

library(mvtnorm)

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

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

M = 5100  # run length
B = 100  # amount of burn-in
X = runif(dim,5,15)  # 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 = numzeroes = 0;

for (i in 1:M) {
    Y = X + rmvnorm(1, sigma = propSigma)  # proposal value
    U = runif(1)  # for accept/reject
    # Use log scale 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);
    # Output progress report.
    if ((i %% 100) == 0)
      cat (" ...",i);
}

cat("\n\nRan Metropolis algorithm for", M, "iterations, with burn-in", B, "\n");
# cat("zeroes rate =", numzeroes/M, "\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=500)$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(hlist[(B+1):M],lag.max=500)