HMCsample.R

D <- 2
tmax <- 100
dt <- 0.1
v <- 0.4
m <- 1

a <- c(1,1)
r <- 0.3

potential <- function(x){
  sum(x**2)+r*sum((x-a)**2)
}

dpotential <- function(x){
  2*x+r*(2*x-2*a)
}

hamiltonian <- function(x,p){
  potencial(x)+p**2/(2*m)
}

prob <- function(x){
  exp(-potential(x))
}

hprob <- function(x,p){
  exp(-hamiltonian(x,p))
}

Gibbsstep <- function(x){
  xn <- x+rnorm(D,0,v)
  x <- ifelse(runif(1,0,1)<prob(xn)/prob(x),xn,x)
}

HMCstep <- function(xp){
  pn <- xp$p +dt*dpotential(xp$x)+rnorm(D,0,v)
  xn <- xp$x -dt*pn+rnorm(D,0,v)
  if(runif(1,0,1)<hprob(xn,pn)/hprob(xp$x,xp$p)){
    xp <-list(x=xn,p=pn)
  }
  xp
}

run <- function(){
  x <- c(4,4)
  p <- c(0.1,0.1)
  xp <- list(x=x0,p=p0)
#Gibbs
  for(t in (1:tmax)){x <- rbind(x,Gibbsstep(x[t,]))}
#HMC
  for(t in (1:tmax)){xp <- rbind(xp,HMCstep(xp[t,]))}

  image(grid,grid,outer(grid,grid,function(a,b){potential(c(a,b))}))
  par(new=T)
  plot(x[,1],x[,2],col="red")
  xh<-numeric()
  for(t in (1:tmax))xh<-rbind(xh,xp$x)
  par(new=T)
  plot(xh[,1],xh[,2],col="blue")

}