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