HMC.one.step: One step of Hamiltonian Monte Carlo

In dereksonderegger/STA578: Functions for a statistical computing course

Description

One step of Hamiltonian Monte Carlo

Usage

HMC.one.step(U, current_q, Eps, L, m = 1)

Arguments

U	A function taking a single argument, which is the position.
Eps	The step size, epsilon.
L	The number of leapfrog steps.
m	The mass of the particle.
q	The current position.

Examples

dtarget <- function(x){
  dmvnorm(x, mean=c(3,10), sigma=matrix(c(1,0,0,1), nrow=2))
}
x1 <- seq(0,6,length=101)
x2 <- seq(7,13, length=101)
contour.data <- expand.grid(x1=x1, x2=x2)
contour.data$Z <- apply(contour.data, MARGIN=1, dtarget)
target.map <- ggplot(contour.data, aes(x=x1, y=x2)) +
  stat_contour(aes(z=Z))
target.map

U <- function(x){ return( -log(dtarget(x))) }
steps <- HMC.one.step(U, c(runif(1,1,5),runif(1,8,12)), Eps=.4, L=20, m=1)
plot_HMC_proposal_path(target.map, steps)


# A little harder...
dtarget <- function(x){
  dmvnorm(x, mean=c(3,10), sigma=matrix(c(3,3,3,7), nrow=2))
}
x1 <- seq(-6,12,length=101)
x2 <- seq(-11,31, length=101)
contour.data <- expand.grid(x1=x1, x2=x2)
contour.data$Z <- apply(contour.data, MARGIN=1, dtarget)
target.map <- ggplot(contour.data, aes(x=x1, y=x2)) +
  stat_contour(aes(z=Z))
target.map

U <- function(x){ return( -log(dtarget(x))) }
steps <- HMC.one.step(U, c(runif(1,-2,8),runif(1,4,16)), Eps=.4, L=40, m=1)
plot_HMC_proposal_path(target.map, steps)


# Now a hard example!
dtarget <- function(x){
  B <- .05
  exp( -x[1]^2 / 200 - (1/2)*(x[2]+B*x[1]^2 -100*B)^2 )
}
x1 <- seq(-20,20, length=201)
x2 <- seq(-15,10, length=201)
contour.data <- expand.grid(x1=x1, x2=x2)
contour.data$Z <- apply(contour.data, MARGIN=1, dtarget)
target.map <- ggplot(contour.data, aes(x=x1, y=x2)) +
  stat_contour(aes(z=Z))
target.map

U <- function(x){ return( -log(dtarget(x))) }
steps <- HMC.one.step(U, c(runif(1,-10,10),runif(1,0,10)), Eps=.15, L=1000, m=1)
plot_HMC_proposal_path(target.map, steps)

dereksonderegger/STA578 documentation built on May 15, 2019, 3:58 a.m.