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

\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("%.6f", x)
    )
    paste(res, collapse = ", ")
  }
}
knit_hooks$set(inline = inline_hook)
@


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

\vspace{0.1 in} \textsc{Metropolis-Hastings for Continuous State Spaces}
\end{center}

We wish to sample from the $Gamma(\gamma,1)$ distribution
\[
\pi(x) = \frac{1}{\Gamma(\gamma)} x^{\gamma-1} e^{-x } \qquad x >0
\]
using the Metropolis-Hastings algorithm.  To implement the algorithm we need to
\begin{enumerate}[1.]
  \item pick a starting value $x_0$
  \item for steps $t=0,1,2,\ldots,N,\ldots$
  \begin{enumerate}[(i)]
    \item propose a new value $z$ from some transition density $q(x_t,z)$, some conditional density given $x_t$,
    \item set $x_{t+1}=z$ with probability
    \[
    \alpha(x_t,z) = \min \left\{1,\frac{\pi(z) q(z,x_t)}{\pi(x_t) q(x_t,z)} \right\}
    \]
    otherwise set $x_{t+1}=x_t$.
  \end{enumerate}
\end{enumerate}
The selected transition density should leave the resulting chain irreducible, aperiodic and positive recurrent, so that $\pi(x)$ is the stationary distribution of the chain.  It is sufficient to achieve a $q(x,z)$ that proposes values on the support of $\pi(x)$.  In the algorithms below, the proposal that sets
\[
z = |x + \delta| \qquad \delta \sim Normal(0,\sigma_\gamma^2)
\]
for some $\sigma_\gamma >0$ is used.  With this choice it follows that $q(x,z) = q(z,x)$, and
    \[
    \alpha(x,z) = \min \left\{1,\frac{\pi(z)}{\pi(x)} \right\}.
    \]
Different results are obtained for different settings.  Performance can be assessed using
\begin{itemize}
\item trace plots of sampled $\{x_t\}$;

\item the autocorrelation function (acf) of the sampled values;

\item the computed effective sample size (ESS).
\end{itemize}
In these examples, $\gamma=2.5$.  First, a run starting from $x_0 = 0$ with $\sigma_\gamma=1$.
<<MH00,comment='+', fig.width=7, fig.height=4.0, fig.align='center'>>=
set.seed(4532)
N<-10000
gam<-2.5
sigg<-1
X<-rep(0,N)
X[1]<-0
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigg))
	al<-min(1,dgamma(Z,gam)/dgamma(x,gam))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}

par(mar=c(4,4,2,0))
plot(X[1:200],type="s",xlab="t",ylab="X",main='Trace plot')
par(mar=c(4,4,2,0))
hist(X,br=seq(0,12,by=0.25),xlab=expression(gamma),main="Histogram of sampled values",
    col="black",border="white",ylim=range(0,800));box()
xv<-seq(0,12,by=0.01);yv<-dgamma(xv,gam);lines(xv,yv*N*0.25,col="red",lwd=2)
par(mar=c(4,4,2,0))
acf(X,lag.max=100,main=' ');title("Autocorrelation for MH sampler",line=1)
library(coda)
effectiveSize(X)
@
Here the effective sample size is \Sexpr{coda::effectiveSize(X)}, which is quite low given the run length.
We now run the algorithm again with $\sigma_\gamma=0.1$

<<MH01,comment='+', fig.width=7, fig.height=4.25, fig.align='center'>>=
set.seed(4532)
N<-10000
gam<-2.5
sigg<-0.1
X<-rep(0,N)
X[1]<-0
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigg))
	al<-min(1,dgamma(Z,gam)/dgamma(x,gam))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}
par(mar=c(4,4,2,0))
plot(X[1:200],type="s",xlab="t",ylab="X",main='Trace plot')
par(mar=c(4,4,2,0))
hist(X,br=seq(0,12,by=0.25),xlab=expression(gamma),main="Histogram of sampled values",
    col="black",border="white",ylim=range(0,800));box()
xv<-seq(0,12,by=0.01);yv<-dgamma(xv,gam)
lines(xv,yv*N*0.25,col="red",lwd=2)
par(mar=c(4,4,2,0))
acf(X,lag.max=100,main=' ');title("Autocorrelation for MH sampler",line=1)
library(coda)
effectiveSize(X)
@
Here the effective sample size is \Sexpr{coda::effectiveSize(X)}, which is extremely poor. Now with $\sigma_\gamma=3$.
<<MH02,comment='+', fig.width=7, fig.height=4.25, fig.align='center'>>=
set.seed(48532)
N<-10000
gam<-2.5; sigg<-3
X<-rep(0,N);X[1]<-0
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigg))
	al<-min(1,dgamma(Z,gam)/dgamma(x,gam))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}

par(mar=c(4,4,2,0))
plot(X[1:200],type="s",xlab="t",ylab="X",main='Trace plot')
par(mar=c(4,4,2,0))
hist(X,br=seq(0,12,by=0.25),xlab=expression(gamma),main="Histogram of sampled values",
    col="black",border="white",ylim=range(0,800));box()
xv<-seq(0,12,by=0.01);yv<-dgamma(xv,gam)
lines(xv,yv*N*0.25,col="red",lwd=2)
par(mar=c(4,4,2,0))
acf(X,lag.max=100,main=' ');title("Autocorrelation for MH sampler",line=1)
effectiveSize(X)
@
Here the effective sample size is \Sexpr{coda::effectiveSize(X)}, which is much better.

\medskip

Finally we start the chain with $x_0 = 20$, and use $\sigma_\gamma=0.1$.  In the plots, the central 0.95 high probability region of the target distribution is marked by the horizontal red lines; the chain should stabilize to spend 95\% of the time in that region.

<<MH03,comment='+', fig.width=7, fig.height=4.25, fig.align='center'>>=
set.seed(48532)
N<-20000
gam<-2.5; sigg<-0.1
X<-rep(0,N);X[1]<-20
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigg))
	al<-min(1,dgamma(Z,gam)/dgamma(x,gam))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}
par(mar=c(4,4,3,0))
plot(X[1:1000],type="s",xlab="t",ylab="X",ylim=range(0,25),main='First 1000 iterations')
abline(h=qgamma(c(0.025,0.975),gam),lty=2,lwd=2,col="red")
par(mar=c(4,4,3,0))
plot(X[1:10000],type="s",xlab="t",ylab="X",ylim=range(0,25),main='First 10000 iterations')
abline(h=qgamma(c(0.025,0.975),gam),lty=2,lwd=2,col="red")
par(mar=c(4,4,3,0))
plot(X[1:N],type="s",xlab="t",ylab="X",ylim=range(0,25),main='First 20000 iterations')
abline(h=qgamma(c(0.025,0.975),gam),lty=2,lwd=2,col="red")
par(mar=c(4,4,3,0))
plot(15001:N,X[15001:N],type="s",xlab="t",ylim=range(0,10),main='Last 5000 iterations')
abline(h=qgamma(c(0.025,0.975),gam),lty=2,lwd=2,col="red")
@
Eventually, the chain recovers from the extreme starting value.  To regard the collected sample as a correlated sample from the target, we should discard the initial iterations.  However, the effective sample sizes are very low due to the high sample autocorrelation in the sampled values.
<<MH05,comment='+', fig.width=7, fig.height=4.25, fig.align='center'>>=
effectiveSize(X[1:1000])  #out of 1000
effectiveSize(X[1:10000]) #out of 10000
effectiveSize(X[1:N])     #out of 20000
effectiveSize(X[15001:N]) #out of 5000
@

\end{document}