# TEST2.R - Test of short cut Metropolis on a multivariate Gaussian. source("../short-met.r") # Description of the multivariate Gaussian distribution. ldim <- 2 # Large dimensions sdim <- 5 # Small dimensions sd <- c (rep(1,ldim), rep(0.1,sdim)) # Numbers of iterations & seed. seed <- 1 cnt <- 1500 cntsc.r <- 2.72*cnt cntsc.ra <- 2.00*cnt cntsc.ra2 <- 2.48*cnt blk <- 600 # Log probability density function for the multivariate Gaussian. gauss.lpr <- function (x) { sum (dnorm(x,0,sd,log=T)) } # Proposal generation. gauss.pf <- function () { rnorm(ldim+sdim,0,1) } # Initial state. initial.x <- sd # Standard Metropolis, w=0.02. set.seed(seed) lpr.evals <- 0 system.time (res.stand.0.02 <- standard.metropolis (initial.x, gauss.lpr, gauss.pf, 0.02, cnt*blk) ) R.stand.0.02 <- res.stand.0.02$states rej.stand.0.02 <- res.stand.0.02$rejected evals.stand.0.02 <- lpr.evals # Standard Metropolis, w=0.1. set.seed(seed) lpr.evals <- 0 system.time (res.stand.0.1 <- standard.metropolis (initial.x, gauss.lpr, gauss.pf, 0.1, cnt*blk) ) R.stand.0.1 <- res.stand.0.1$states rej.stand.0.1 <- res.stand.0.1$rejected evals.stand.0.1 <- lpr.evals # Standard Metropolis, w=0.5. set.seed(seed) lpr.evals <- 0 system.time (res.stand.0.5 <- standard.metropolis (initial.x, gauss.lpr, gauss.pf, 0.5, cnt*blk) ) R.stand.0.5 <- res.stand.0.5$states rej.stand.0.5 <- res.stand.0.5$rejected evals.stand.0.5 <- lpr.evals # Standard Metropolis, with w=0.02 for 200 updates, w=0.1 for 200 updates, # and w=0.5 for 200 updates. set.seed(seed) lpr.evals <- 0 system.time (res.stand.all <- standard.sequences (initial.x, gauss.lpr, gauss.pf, rep(c(0.02,0.1,0.5),cnt), rep(c(200,200,200),cnt)) ) R.stand.all <- res.stand.all$states rej.stand.all <- res.stand.all$rejected step.stand.all <- res.stand.all$stepsize evals.stand.all <- lpr.evals # Short cut Metropolis, with w=0.02 for 60 updates, w=0.1 for 150 updates, # and w=0.5 for 390 updates, reversing on all rejections (except for w=0.02). set.seed(seed) lpr.evals <- 0 system.time (res.short.r <- short.cut.sequences (initial.x, gauss.lpr, gauss.pf, rep(c(0.02,0.1,0.5),cntsc.r), rep(c(10,10,10),cntsc.r), rep(c(6,15,39),cntsc.r), c(0,0,0), c(10,9,9)) ) R.short.r <- res.short.r$states1 rej.short.r <- res.short.r$rejected copied.short.r <- res.short.r$copied index.short.r <- res.short.r$index step.short.r <- res.short.r$stepsize evals.short.r <- lpr.evals # Short cut Metropolis, with w=0.02 for 200 updates, w=0.1 for 200 updates, # and w=0.5 for 200 updates, reversing on all rejections (except for w=0.02) # or on all acceptances (except for w=0.1). set.seed(seed) lpr.evals <- 0 system.time (res.short.ra <- short.cut.sequences (initial.x, gauss.lpr, gauss.pf, rep(c(0.02,0.1,0.5),cntsc.ra), rep(c(10,10,10),cntsc.ra), rep(c(20,20,20),cntsc.ra), c(1,1,0), c(10,9,9)) ) R.short.ra <- res.short.ra$states1 rej.short.ra <- res.short.ra$rejected copied.short.ra <- res.short.ra$copied index.short.ra <- res.short.ra$index step.short.ra <- res.short.ra$stepsize evals.short.ra <- lpr.evals # Short cut Metropolis, with w=0.02 for 200 updates, w=0.1 for 200 updates, # and w=0.5 for 200 updates, reversing on all rejections (except for w=0.02), # or on all or all but one acceptances (except for w=0.1). set.seed(seed) lpr.evals <- 0 system.time (res.short.ra2 <- short.cut.sequences (initial.x, gauss.lpr, gauss.pf, rep(c(0.02,0.1,0.5),cntsc.ra), rep(c(10,10,10),cntsc.ra2), rep(c(20,20,20),cntsc.ra2), c(2,2,0), c(10,9,9)) ) R.short.ra2 <- res.short.ra2$states1 rej.short.ra2 <- res.short.ra2$rejected copied.short.ra2 <- res.short.ra2$copied index.short.ra2 <- res.short.ra2$index step.short.ra2 <- res.short.ra2$stepsize evals.short.ra2 <- lpr.evals # Compute autocorrelation times. ac.time <- function (x, lag.max=8000) { 2 * sum(acf(as.numeric(x),plot=F,lag.max=lag.max)$acf) - 1 } ac.times.m <- c(ac.time(R.stand.0.02[,1],lag.max=12000), ac.time(R.stand.0.1[,1]), ac.time(R.stand.0.5[,1],lag.max=12000), ac.time(R.stand.all[,1]), ac.time(R.short.r[,1]), ac.time(R.short.ra[,1]), ac.time(R.short.ra2[,1])) # Compute and print table of results. results <- data.frame (list ( method = c("stand.0.02", "stand.0.1", "stand.0.5", "stand.all", "short.r", "short.ra", "short.ra2"), evals = c(evals.stand.0.02, evals.stand.0.1, evals.stand.0.5, evals.stand.all, evals.short.r, evals.short.ra, evals.short.ra2), rej = round(c(mean(rej.stand.0.02), mean(rej.stand.0.1), mean(rej.stand.0.5), mean(rej.stand.all), mean(rej.short.r), mean(rej.short.ra), mean(rej.short.ra2)),3), mean = round(c(mean(R.stand.0.02[,1]), mean(R.stand.0.1[,1]), mean(R.stand.0.5[,1]), mean(R.stand.all[,1]), mean(R.short.r[,1]), mean(R.short.ra[,1]), mean(R.short.ra2[,1])),4), ac.m = round(ac.times.m,2), s.e.m = round( sqrt(1/(blk*c(cnt,cnt,cnt,cnt,cntsc.r,cntsc.ra,cntsc.ra2)/ac.times.m)), 4) )) print(results) # Function to compute fractions of updates that are copied, separately # for each sort of short-cut sequence. copy.frac <- function (copied, cnt) { ix <- c() for (i in 1:length(cnt)) { ix <- c(ix,rep(i,cnt[i])) } ix <- rep(ix,length=length(copied)) r <- rep(NA,length(cnt)) for (i in 1:length(cnt)) { r[i] <- mean(copied[ix==i]) } r }