HMC2: Functions for simple HMC simulations

In rmcelreath/rethinking: Statistical Rethinking book package

Functions for simple HMC simulations

Description

Conduct simple Hamiltonian Monte Carlo simulations.

Usage

HMC2(U, grad_U, epsilon, L, current_q , ... )
HMC_2D_sample( n=100 , U , U_gradient , step , L , start=c(0,0) , 
  xlim=c(-5,5) , ylim=c(-4,4) , xlab="x" , ylab="y" , 
  draw=TRUE , draw_contour=TRUE , nlvls=15 , adj_lvls=1 , ... )

Arguments

U	Function to return log density
grad_U	Function to return gradient of U
epsilon	Size of leapfrog step
L	Number of leapfrog steps
current_q	Initial position
n	Number of transitions to sample
step	Size of leapfrog step
start	Initial position
xlim	Horizontal boundaries of plot, for when draw is TRUE
ylim	Vertical boundaries of plot, for when draw is TRUE
draw	Whether to draw the samples and trajectories
draw_contour	Whether to draw contour of density
nlvls	Number of contour levels
adj_lvls	Factor to multiply levels by
...	Optional arguments to pass to density and gradient functions, typically optional parameters

Details

These functions provide simple Hamiltonian Monte Carlo simulations.

HMC2 is based on Neals's 2010 code (see citation below), but returns the trajectories.

HMC_2D_sample simulates multiple sequential trajectories and can also plot the simulation. See examples below.

To use either function, the user must supply custom density and gradient functions.

Author(s)

Richard McElreath

See Also

Radford M. Neal, 2010. MCMC using Hamiltonian dynamics. The Handbook of Markov Chain Monte Carlo.

Examples

# Devil's Funnel from Chapter 13
U_funnel <- function( q , s=3 ) {
    v <- q[2]
    x <- q[1]
    U <- sum( dnorm(x,0,exp(v),log=TRUE) ) + dnorm(v,0,s,log=TRUE)
    return( -U )
}
U_funnel_gradient <- function( q , s=3 ) {
    v <- q[2]
    x <- q[1]
    Gv <- (-v)/s^2 - length(x) + exp(-2*v)*sum( x^2 ) #dU/dv
    Gx <- -exp(-2*v)*x #dU/dx
    return( c( -Gx , -Gv ) ) # negative bc energy is neg-log-prob
}
HMC_2D_sample( n=3 , U=U_funnel , U_gradient=U_funnel_gradient , 
  step=0.2 , L=10 , ylab="v"  , adj_lvls=1/12 )

# Same, but with non-centered parameterization

U_funnel_NC <- function( q , s=3 ) {
    v <- q[2]
    z <- q[1]
    U <- sum( dnorm(z,0,1,log=TRUE) ) + dnorm(v,0,s,log=TRUE)
    return( -U )
}
U_funnel_NC_gradient <- function( q , s=3 ) {
    v <- q[2]
    z <- q[1]
    Gv <- (-v)/s^2 #dU/dv
    Gz <- (-z) #dU/dz
    return( c( -Gz , -Gv ) ) # negative bc energy is neg-log-prob
}
HMC_2D_sample( n=4 , U=U_funnel_NC , U_gradient=U_funnel_NC_gradient , 
  step=0.2 , L=15 , ylab="v" , xlab="z" , xlim=c(-2,2) , nlvls=12 , adj_lvls=1 )

rmcelreath/rethinking documentation built on Sept. 18, 2023, 2:01 p.m.