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

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

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

mydnorm = function( yourx, yourmean, yourvar, log=FALSE ) {
  thedim = length(yourx)
  thelog = -thedim/2 * log(2*3.1415926) -
        1/2*log(det(yourvar)) -
        1/2 * ( (yourx-yourmean) %*% solve(yourvar) %*% (yourx-yourmean) )
  if (log)
    return( thelog )
  else
    return( exp(thelog) )
}

library(MASS)

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

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

M = 11000  # run length
B = 1000  # 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 + mvrnorm(1, mu=rep(0,dim), 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)