Static Hamiltonian Monte Carlo

ggplot2 is used for plotting, tidyr for manipulating data frames

library(ggplot2)
theme_set(theme_minimal())
library(tidyr)
library(ggforce)
library(MASS)
library(rprojroot)
library(rstan)
root<-has_file(".BDA_R_demos_root")$make_fix_file()

Parameters of a normal distribution used as a toy target distribution

y1 <- 0
y2 <- 0
r <- 0.99
S <- diag(2)
S[1, 2] <- r
S[2, 1] <- r

Draw samples from the toy distribution to visualize 90% HPD interval with ggplot's stat_ellipse()

dft <- data.frame(mvrnorm(100000, c(0, 0), S))

see BDA3 p. 85 for how to compute HPD for multivariate normal in 2d-case contour for 90% HPD is an ellipse, whose semimajor axes can be computed from the eigenvalues of the covariance matrix scaled by a value selected to get ellipse match the density at the edge of 90% HPD. Angle of the ellipse could be computed from the eigenvectors, but since the marginals are same we know that angle is pi/4 Starting value of the chain

t1 <- -2.5
t2 <- 2.5

Number of iterations.

M <- 5000

Insert your own HMC sampling here

# Allocate memory for the samples
tt <- matrix(rep(0, 2*M), ncol = 2)
tt[1,] <- c(t1, t2)    # Save starting point
# For demonstration load pre-computed values
# Replace this with your algorithm!
# tt is a M x 2 array, with M samples of both theta_1 and theta_2
load(root("demos_ch11","demo12_1b.RData"))
# Here we have intentionally used a very small step size for smooth
# simulations, but for more efficient simulations larger step size
# could be used

The rest is for illustration Creat data frame

df <- data.frame(id=rep(1,4000),
                 iter=rep(1:100, each=40), 
                 th1 = tt[1:4000, 1],
                 th2 = tt[1:4000, 2],
                 th1l = c(tt[1, 1], tt[1:(4000-1), 1]),
                 th2l = c(tt[1, 2], tt[1:(4000-1), 2]))

Take the first 1000 draws after warmup of 1

dfs <- data.frame(th1 = tt[seq(41,40001,by=40), 1], th2 = tt[seq(41,40001,by=40), 2])

Base for the plot

# Labels and frame indices for the plot
labs1 <- c('Samples', 'Steps of the sampler', '90% HPD')
p0 <- ggplot() +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '3'), level = 0.9) +
  coord_cartesian(xlim = c(-3, 3), ylim = c(-3, 3)) +
  labs(x = 'theta1', y = 'theta2') +
    scale_color_manual(values = c('red', 'forestgreen','blue'), labels = labs1) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA, NA), linetype = c(0, 1, 1)))) +
    theme(legend.position = 'bottom', legend.title = element_blank())

Plot several iterations

for (ind in seq(40,400,by=40)) {
    pp <- p0 +   geom_point(data = df[(ind-40+1):ind,], 
                      aes(th1, th2, color ='1'), alpha=0.3, size=1) +
        geom_segment(data = df[(ind-40+1):ind,], aes(x = th1, xend = th1l, color = '2',
                                                     y = th2, yend = th2l),
                     alpha=0.5) +
        geom_point(data = df[seq(1,ind+1,by=40),], 
                   aes(th1, th2, color ='1'), size=2)
    print(pp)
}