# 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=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(hlist[(B+1):M],lag.max=250)