Banana-log-target density

Metropolis-Hastings vs Hamiltonian MC

1 Banana log-target density

The target density for \((x,y)\) is given by \[ p(x,y) \propto \exp\{-0.5[x^2+(y-bx^2]^2)\}. \]

# Banana log-density
logp = function(xx) {
  b = 2
  x = xx[1]
  y = xx[2]
  return(-0.5*(x^2+(y-b*x^2)^2))
}

N    = 200
x    = seq(-3,3,length=N)
y    = seq(-3,10,length=N)
post = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N)
    post[i,j] = logp(c(x[i],y[j]))

par(mfrow=c(1,1))
contour(x,y,exp(post),drawlabels=FALSE,xlab="x",ylab="y",main="Target distribution")

2 Metropolis-Hastings sampler

metropolis = function(logp,start=c(0,0),iter=niter,warmup=burnin,proposal_sd=0.4,seed=123){
  samples     = matrix(0,iter,2)
  samples[1,] = start
  logp_curr   = logp(start)
  
  for (i in 2:iter) {
    proposal    =  rnorm(2,samples[i-1,],proposal_sd)
    logp_prop   = logp(proposal)
    accept_prob = exp(logp_prop - logp_curr)
    
    if (runif(1) < accept_prob) {
      samples[i,] = proposal
      logp_curr   = logp_prop
    } else {
      samples[i,] = samples[i-1,]
    }
  }
  samples[(warmup+1):iter,]
}

seed = 31416
set.seed(seed)
niter  = 20000
burnin = 10000
M      = niter - burnin

mh_samples = metropolis(logp)

2.1 Trace plots and ACF plots

par(mfrow=c(2,2))
ts.plot(mh_samples[,1],xlab="Iterations",ylab="x",main="Metropolis-Hastings")
acf(mh_samples[,1],main="")
ts.plot(mh_samples[,2],xlab="Iterations",ylab="y")
acf(mh_samples[,2],main="")

2.2 Target draws

par(mfrow=c(1,1))
xrange = range(mh_samples[,1],x)
yrange = range(mh_samples[,2],y)
plot(mh_samples, pch=16, cex=0.3, main="Metropolis-Hastings",xlab="x",ylab="y",
     xlim=xrange,ylim=yrange)
contour(x,y,exp(post),add=TRUE,col=2,lwd=2)

2.3 Effective sample size

ESS.mh = c(M/sum(acf(mh_samples[,1],plot=FALSE)$acf^2),M/sum(acf(mh_samples[,2],plot=FALSE)$acf^2))
c(M,round(ESS.mh))

## [1] 10000   696   375

3 Hamiltonian Monte Carlo

library(rstan)

## Loading required package: StanHeaders

## 
## rstan version 2.32.7 (Stan version 2.32.2)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)

stan_code = "
data {
  int<lower=0> N;
}
parameters {
  real x;
  real y;
}
model {
  real b = 2;
  x ~ normal(0, 1);
  y ~ normal(b * x^2, 1);
}
"

3.1 Compiling and sampling

sm = stan_model(model_code = stan_code)

fit = sampling(sm,data=list(N=1),iter=niter,warmup=burnin,chains=1,seed=seed)

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 1: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 1: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 1: Iteration:  6000 / 20000 [ 30%]  (Warmup)
## Chain 1: Iteration:  8000 / 20000 [ 40%]  (Warmup)
## Chain 1: Iteration: 10000 / 20000 [ 50%]  (Warmup)
## Chain 1: Iteration: 10001 / 20000 [ 50%]  (Sampling)
## Chain 1: Iteration: 12000 / 20000 [ 60%]  (Sampling)
## Chain 1: Iteration: 14000 / 20000 [ 70%]  (Sampling)
## Chain 1: Iteration: 16000 / 20000 [ 80%]  (Sampling)
## Chain 1: Iteration: 18000 / 20000 [ 90%]  (Sampling)
## Chain 1: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.059 seconds (Warm-up)
## Chain 1:                0.071 seconds (Sampling)
## Chain 1:                0.13 seconds (Total)
## Chain 1:

## Warning: There were 15 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

posterior = as.data.frame(fit)

3.2 Trace plots and ACF plots

par(mfrow=c(2,2))
ts.plot(posterior$x,xlab="Iterations",ylab="x",main="Hamiltonian MC")
acf(posterior$x,main="")
ts.plot(posterior$y,xlab="Iterations",ylab="y")
acf(posterior$y,main="")

3.3 Target draws

par(mfrow=c(1,1))
plot(posterior$x, posterior$y, pch=16, cex=0.3,main="HMC (Stan)",xlab="x",ylab="y")
contour(x,y,exp(post),add=TRUE,col=2,lwd=2)

3.4 Effective sample size

ESS.hmc = c(M/sum(acf(posterior$x,plot=FALSE)$acf^2),M/sum(acf(posterior$y,plot=FALSE)$acf^2))
c(M,round(ESS.hmc))

## [1] 10000  2642  2500