# SCRIPT FOR PROGRAMMING ASSIGNMENT 2. source("hmc.r") cat("\nDATA1:\n\n") data1 <- as.matrix(read.table("data1")) s1 <- c(0.4,0.8,0.1,0.1,0.1) set.seed(1) iters <- 1100 accept_count <- 0 prev <- rep(0,5) res1 <- matrix(0,iters,5) for (i in 1:iters) { e <- runif(1,0.9,1.1)*0.3 prev <- res1[i,] <- HMC (log_pi, grad_log_pi, 1/s1^2, e, 20, prev, data1) } cat("acceptance rate:",round(accept_count/iters,2),"\n") pdf("plots1.pdf") discard <- 1:100 plot(as.data.frame(res1[-discard,,drop=FALSE]),pch=".") dev.off() pdf("plots1b.pdf") par(mfrow=c(2,1)) plot(res1[-discard,5],pch=".",ylim=c(-0.8,-1.3)) plot(res1[-discard,1],pch=".") dev.off() cat("\nDATA2:\n\n") data2 <- as.matrix(read.table("data2")) s2 <- c(0.5,0.6,0.2,0.1,0.1) set.seed(1) iters <- 1100 accept_count <- 0 prev <- rep(0,5) res2 <- matrix(0,iters,5) for (i in 1:iters) { e <- runif(1,0.9,1.1)*0.04 prev <- res2[i,] <- HMC (log_pi, grad_log_pi, 1/s2^2, e, 500, prev, data2) } cat("acceptance rate:",round(accept_count/iters,2),"\n") pdf("plots2.pdf") discard <- 1:100 plot(as.data.frame(res2[-discard,,drop=FALSE]),pch=".") dev.off() pdf("plots2b.pdf") par(mfrow=c(2,1)) plot(res2[-discard,5],pch=".",ylim=c(1,11)) plot(res2[-discard,1],pch=".") dev.off() cat("\nDATA3:\n\n") data3 <- as.matrix(read.table("data3")) s3 <- c(0.5,0.5,0.4,0.3,0.3) set.seed(1) iters <- 1100 accept_count <- 0 prev <- rep(0,5) res3 <- matrix(0,iters,5) for (i in 1:iters) { e <- runif(1,0.9,1.1)*0.15 prev <- res3[i,] <- HMC (log_pi, grad_log_pi, 1/s3^2, e, 20, prev, data3) } cat("acceptance rate:",round(accept_count/iters,2),"\n") pdf("plots3.pdf") discard <- 1:100 plot(as.data.frame(res3[-discard,,drop=FALSE]),pch=".") dev.off() pdf("plots3b.pdf") par(mfrow=c(2,1)) plot(res3[-discard,5],pch=".",ylim=c(-0.5,2.0)) plot(res3[-discard,1],pch=".") dev.off()