# TEST1.R - Test of short cut Metropolis on a 1D mixture of normals. source("../short-met.r") # Numbers of iterations & seed. seed <- 3 cnt <- 10000 cntsc.r <- 1.65*cnt cntsc.ra <- 1.80*cnt blk <- 120 # Log probability density function for the mixture distribution. mix.lpr <- function (x) { log (0.5 * dnorm(x,0,10) + 0.5 * dnorm(x,10,1)) } # Mean, variance, and probability of value <10 for the mixture distribution. mix.mean <- 5 mix.var <- 5^2 + 0.5*(100+1) mix.pr <- 0.5*pnorm(10,0,10) + 0.5*pnorm(10,10,1) # Proposal generation. mix.pf <- function () { rnorm(1,0,1) } # Standard Metropolis, w=2. set.seed(seed) lpr.evals <- 0 res.stand.2 <- standard.metropolis (0, mix.lpr, mix.pf, 2, cnt*blk) R.stand.2 <- res.stand.2$states rej.stand.2 <- res.stand.2$rejected evals.stand.2 <- lpr.evals # Standard Metropolis, w=20. set.seed(seed) lpr.evals <- 0 res.stand.20 <- standard.metropolis (0, mix.lpr, mix.pf, 20, cnt*blk) R.stand.20 <- res.stand.20$states rej.stand.20 <- res.stand.20$rejected evals.stand.20 <- lpr.evals # Naive adaptive Metropolis, w0=2 and w1=20. set.seed(seed) lpr.evals <- 0 res.naive <- naive.adaptive.metropolis (0, mix.lpr, mix.pf, 2, 20, 10, 0.5, cnt*blk) R.naive <- res.naive$states rej.naive <- res.naive$rejected step.naive <- res.naive$stepsize evals.naive <- lpr.evals # Short cut Metropolis, with w=2 for 30 updates and w=20 for 90 updates, # reversing on all rejections. set.seed(seed) lpr.evals <- 0 res.short.r <- short.cut.sequences (0, mix.lpr, mix.pf, rep(c(2,20),cntsc.r), rep(c(5,5),cntsc.r), rep(c(6,18),cntsc.r)) 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=2 for 60 updates and w=20 for 60 updates, # reversing on all rejections or all acceptances. set.seed(seed) lpr.evals <- 0 res.short.ra <- short.cut.sequences (0, mix.lpr, mix.pf, rep(c(2,20),cntsc.ra), rep(c(5,5),cntsc.ra), rep(c(12,12),cntsc.ra), 1) 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 # Compute autocorrelation times. ac.time <- function (x, lag.max=500) { 2 * sum(acf(as.numeric(x),plot=F,lag.max=lag.max)$acf) - 1 } ac.times.m <- c(ac.time(R.stand.2), ac.time(R.stand.20), ac.time(R.naive), ac.time(R.short.r), ac.time(R.short.ra)) ac.times.pr <- c(ac.time(R.stand.2<10), ac.time(R.stand.20<10), ac.time(R.naive<10), ac.time(R.short.r<10), ac.time(R.short.ra<10)) # Compute and print table of results. results <- data.frame (list ( method = c("stand.2", "stand.20", "naive", "short.r", "short.ra"), evals = c(evals.stand.2, evals.stand.20, evals.naive, evals.short.r, evals.short.ra), rej = round(c(mean(rej.stand.2), mean(rej.stand.20), mean(rej.naive), mean(rej.short.r), mean(rej.short.ra)),3), mean = round(c(mean(R.stand.2), mean(R.stand.20), mean(R.naive), mean(R.short.r), mean(R.short.ra)),4), ac.m = round(ac.times.m,2), s.e.m = round( sqrt(mix.var/(blk*c(cnt,cnt,cnt,cntsc.r,cntsc.ra)/ac.times.m)),4), pr.lt.10 = round(c(mean(R.stand.2<10), mean(R.stand.20<10), mean(R.naive<10), mean(R.short.r<10), mean(R.short.ra<10)),4), ac.pr = round(ac.times.pr,2), s.e.pr = round( sqrt(mix.pr*(1-mix.pr)/(blk*c(cnt,cnt,cnt,cntsc.r,cntsc.ra)/ac.times.pr)),4) )) print(results) # Produce some postscript plots. plt.portion <- function (res, r2, L, blks) { r1 <- 1+L*(range(r2)-1) r1 <- r1[1]:r1[2] s1 <- as.vector(res$states1[r1,]) s2 <- as.vector(res$states2[r2,]) plot (range(r1)-1, range(c(s1,s2)), type="n", ylim=c(-20,25), xlab="", ylab="", xaxs="i", xpd=NA) abline (v=seq(r1[1]-1,r1[length(r1)]-1,by=5), col="gray") abline (v=blks, col="black") points (r1[2:length(r1)]-1, s1[2:length(r1)], pch=20, xpd=NA, col=c("gray45","gray80")[1+res$copied[r1[2:length(r1)]-1]]) points (L*(r2-1), s2, pch=1, xpd=NA, col="black") abline(h=c(8,12)) } postscript("test1r.ps",pointsize=9,width=6.9,height=3.5,horiz=F) plt.portion(res.short.r,1:49,5,c(0,30,120,150,240)) dev.off() postscript("test1ra.ps",pointsize=9,width=6.9,height=3.5,horiz=F) plt.portion(res.short.ra,145:193,5,720+c(0,60,120,180,240)) dev.off()