\documentclass[notitlepage]{article}
\usepackage{../../Math598}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumerate}
\usepackage{multirow}
\usepackage{bbm}
\usepackage{bbold}
\usepackage{bm}



\usepackage{amsfonts,amsmath,amssymb}

\setlength{\parindent}{0pt}


\lstloadlanguages{R}

\definecolor{keywordcolor}{rgb}{0,0.6,0.6}
\definecolor{delimcolor}{rgb}{0.461,0.039,0.102}
\definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432}

\lstdefinestyle{Rsettings}{
  basicstyle=\ttfamily,
  breaklines=true,
  showstringspaces=false,
  keywords={if, else, function, theFunction, tmp}, % Write as many keywords
  otherkeywords={},
  commentstyle=\itshape\color{Rcommentcolor},
  keywordstyle=\color{keywordcolor},
  moredelim=[s][\color{delimcolor}]{"}{"},
}

\lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}}

\begin{document}

<<setup, cache=FALSE, echo=FALSE>>=
library(knitr)
# global chunk options
opts_chunk$set(cache=TRUE, autodep=TRUE)
options(scipen=999)
options(repos=c(CRAN="https://cloud.r-project.org/"))
inline_hook <- function (x) {
  if (is.numeric(x)) {
    # ifelse does a vectorized comparison
    # If integer, print without decimal; otherwise print two places
    res <- ifelse(x == round(x),
      sprintf("%6i", x),
      sprintf("%.3f", x)
    )
    paste(res, collapse = ", ")
  }
}
knit_hooks$set(inline = inline_hook)
@


\begin{center}
\textsclarge{MATH 598: TOPICS IN STATISTICS}

\vspace{0.1 in} \textsc{Hamiltonian MCMC: 1d example}
\end{center}

\textbf{Example:}  Suppose $Y \equiv Gamma(\alpha,1)$, and consider the pdf of $X = \log Y$.  Then for $x \in \R$
\[
f_X(x) = e^x f_Y(e^x) = \frac{1}{\Gamma(\alpha)} \exp\{ \alpha x - e^x\}.
\]
We set $\pi_X(x) \equiv f_X(x)$ so that
\begin{align*}
\log \pi_X(x) & = \alpha x - e^x - \log \Gamma(\alpha)
\end{align*}
and hence
\[
U(x) = - \alpha x + e^x + \text{constant}
\]
with the constant chosen to that $U(x) \geq 0$ for all $x$.  As 
\[
\dot{U}(x) = -\alpha + e^x
\]
the function $U(x)$ is minimized at $x=\log \alpha$, so we define
\[
U(x) = - \alpha x + e^x + \alpha \log \alpha - \alpha.
\]
Define also
\[
K(v) = \frac{1}{2} v^2.
\]
We consider the leapfrog algorithm for proposing moves on the joint pdf
\[
\pi_{X,V}(x,v) \varpropto \exp\left\{-U(x) - K(v)  \right\} \varpropto \exp\left\{\alpha x - e^x - \frac{1}{2} v^2  \right\}
\]
We take $\alpha=10$, and implement the Leapfrog algorithm from the same fixed starting point $x=1$, $v=0$ for different settings of $\delta$ and $L=500$.  The path traced out by the leapfrog algorithm is plotted in white.

<<H1,comment='+', fig.width=7, fig.height=5.0, fig.align='center', echo=TRUE>>=
set.seed(23984)
al<-10

library(mvnfast,fields,quietly=TRUE)
xvec<-seq(-1,4,by=0.01)
yvec<-seq(-4,4,by=0.01)
dfunc<-function(xv,vv,av){
    dval<-av*xv -exp(xv) - lgamma(av) - 0.5*vv^2
}
f <- Vectorize(dfunc,vectorize.args=c("xv","vv"))
zmat<-outer(xvec,yvec,f,av=al)
par(pty='s',mar=c(4,3,2,2))
colfunc <- colorRampPalette(c("blue","lightblue","white","yellow","orange","red"))
image.plot(xvec,yvec,zmat,col=colfunc(200),
           xlab=expression(x),ylab=expression(v),cex.axis=0.8)
contour(xvec,yvec,zmat,add=T,levels=seq(-50,0,by=4))
title(expression(paste('Log ',pi[paste(X,',',V,sep='')](x,v))))
#Leapfrog run from fixed starting position
L<-200
del<-0.01
x<-1
v<-0
xpath<-matrix(0,nrow=L,ncol=2)
Ham<-U<-K<-rep(0,L)
xpath[1,]<-c(x,v)
U[1]<--al*x+exp(x)+al*log(al)-al
K[1]<-0.5*v^2
Ham[1]<-U[1]+K[1]
for(l in 1:(L-1)){
    vstar<-v-0.5*del*(-1)*(al-exp(x))
    xnew<-x+del*vstar
    vnew<-vstar-0.5*del*(-1)*(al-exp(x))
    x<-xnew
    v<-vnew
    xpath[l+1,]<-c(x,v)
    U[l+1]<--(al*x-exp(x)-lgamma(al))+al*log(al)-al
    K[l+1]<-0.5*v^2
    Ham[l+1]<-U[l+1]+K[l+1]
}
for(i in 2:L){lines(xpath[c(i-1,i),1],xpath[c(i-1,i),2],lty=2,col='white')}
points(xpath[1,1],xpath[1,2],col='green',pch=20)
points(xpath[L,1],xpath[L,2],col='white',pch=20);

@

<<H2,comment='+', fig.width=7, fig.height=4.0, fig.align='center', echo=FALSE>>=
par(mar=c(4,3,1,0))
plot(Ham,type='l',ylim=range(0,30),xlab='Leapfrog step')
lines(U,lty=2,col='red')
lines(K,lty=2,col='blue')
legend(15,30,c(expression(H(x,v)),expression(U(x)),expression(K(v))),
        lty=c(1,2,2),col=c('black','red','blue'))
@

Now with $\delta=0.02$.

<<H3,comment='+', fig.width=7, fig.height=5.0, fig.align='center', echo=FALSE>>=
par(pty='s',mar=c(4,3,2,2))
image.plot(xvec,yvec,zmat,col=colfunc(200),
           xlab=expression(x),ylab=expression(v),cex.axis=0.8)
contour(xvec,yvec,zmat,add=T,levels=seq(-50,0,by=4))
title(expression(paste('Log ',pi[paste(X,',',V,sep='')](x,v))))
#Leapfrog run from fixed starting position
L<-200
del<-0.02
x<-1
v<-0
xpath<-matrix(0,nrow=L,ncol=2)
Ham<-U<-K<-rep(0,L)
xpath[1,]<-c(x,v)
U[1]<--al*x+exp(x)+al*log(al)-al
K[1]<-0.5*v^2
Ham[1]<-U[1]+K[1]
for(l in 1:(L-1)){
    vstar<-v-0.5*del*(-1)*(al-exp(x))
    xnew<-x+del*vstar
    vnew<-vstar-0.5*del*(-1)*(al-exp(x))
    x<-xnew
    v<-vnew
    xpath[l+1,]<-c(x,v)
    U[l+1]<--(al*x-exp(x)-lgamma(al))+al*log(al)-al
    K[l+1]<-0.5*v^2
    Ham[l+1]<-U[l+1]+K[l+1]
}
for(i in 2:L){lines(xpath[c(i-1,i),1],xpath[c(i-1,i),2],lty=2,col='white')}
points(xpath[1,1],xpath[1,2],col='green',pch=20)
points(xpath[L,1],xpath[L,2],col='white',pch=20);
@
These plots demonstrate that if moving quickly around the parameter space is a priority, then the choice of $L$ makes a big difference; in the next four plots we demonstrate where the proposal lands with $L=50,100,150$ and 200.

<<H4,comment='+', fig.width=7, fig.height=6.0, fig.align='center', echo=FALSE>>=
par(pty='s',mar=c(2,0,2,0),mfrow=c(2,2))
for(j in 1:4){
    sL<-j*(L/4)
    image(xvec,yvec,zmat,col=colfunc(200),
           xlab=expression(x),ylab=expression(v),cex.axis=0.8)
    contour(xvec,yvec,zmat,add=T,levels=seq(-50,0,by=4))
    for(i in 2:sL){lines(xpath[c(i-1,i),1],xpath[c(i-1,i),2],lty=2,col='white')}
    points(xpath[1,1],xpath[1,2],col='green',pch=20)
    points(xpath[sL,1],xpath[sL,2],col='white',pch=20)
    title(substitute(paste('L = ',lval),list(lval=sL)))
}
@

\medskip

The leapfrog proposal traverses the contours of $\pi_{X,V}(x,v)$: with $L=50$ or $150$, the landing point is distant from the starting point, whereas for $L=100$ or $200$ the landing point is quite close to the starting point.

\pagebreak

Of course, this is starting point dependent: if we start from $(x,v)=(1.5,1)$ then the leapfrog trajectory is different.
<<H5,comment='+', fig.width=7, fig.height=6.0, fig.align='center', echo=FALSE>>=
par(pty='m',mar=c(4,3,2,2))
#image.plot(xvec,yvec,zmat,col=colfunc(200),
#           xlab=expression(x),ylab=expression(v),cex.axis=0.8)
#contour(xvec,yvec,zmat,add=T,levels=seq(-50,0,by=4))
#title(expression(paste('Log ',pi[paste(X,',',V,sep='')](x,v))))
#Leapfrog run from fixed starting position
L<-200
del<-0.02
x<-1.5
v<-1
xpath<-matrix(0,nrow=L,ncol=2)
Ham<-U<-K<-rep(0,L)
xpath[1,]<-c(x,v)
U[1]<--al*x+exp(x)+al*log(al)-al
K[1]<-0.5*v^2
Ham[1]<-U[1]+K[1]
for(l in 1:(L-1)){
    vstar<-v-0.5*del*(-1)*(al-exp(x))
    xnew<-x+del*vstar
    vnew<-vstar-0.5*del*(-1)*(al-exp(x))
    x<-xnew
    v<-vnew
    xpath[l+1,]<-c(x,v)
    U[l+1]<--(al*x-exp(x)-lgamma(al))+al*log(al)-al
    K[l+1]<-0.5*v^2
    Ham[l+1]<-U[l+1]+K[l+1]
}
#for(i in 2:L){lines(xpath[c(i-1,i),1],xpath[c(i-1,i),2],lty=2,col='white')}
#points(xpath[1,1],xpath[1,2],col='green',pch=20)
#points(xpath[L,1],xpath[L,2],col='white',pch=20);
par(pty='s',mar=c(2,0,2,0),mfrow=c(2,2))
for(j in 1:4){
    sL<-j*(L/4)
    image(xvec,yvec,zmat,col=colfunc(200),
           xlab=expression(x),ylab=expression(v),cex.axis=0.8)
    contour(xvec,yvec,zmat,add=T,levels=seq(-50,0,by=4))
    for(i in 2:sL){lines(xpath[c(i-1,i),1],xpath[c(i-1,i),2],lty=2,col='white')}
    points(xpath[1,1],xpath[1,2],col='green',pch=20)
    points(xpath[sL,1],xpath[sL,2],col='white',pch=20)
    title(substitute(paste('L = ',lval),list(lval=sL)))
}
@

\end{document}