library(tidyverse)
library(VGAM)

Attempt to implement a sampler for 9.3 or 9.4 in Rizzo’s text. Post your code / issues here, and discuss with your peers.

The CDF of the Laplace Distribution is:

\(F(x)=\begin{cases} \frac{1}{2}e^{\frac{x-\mu}{b}}, & x < \mu \\ 1-\frac{1}{2}e^{-\frac{x-\mu}{b}}, & x \geq \mu \end{cases}\)

Setting this equal to a random number generated via \(U[0,1]=R\), I get:

\[X_R = \mu-b \frac{R-\frac{1}{2}}{|R-\frac{1}{2}|} \log(1-2|R-\frac{1}{2}|)\]

my_laplace <- function(x, b = 0, mu = 1){
  #R <- runif(1)
  #val <- mu - b * (R-.5)/abs(R-.5) * log(1 - 2 * abs(R-.5)) 
  cntrl <- dlaplace(x, location = b, scale = mu)
  #return(list( val = val, c = cntrl))
  return(cntrl)
}

ttt <- runif(5)
test <- my_laplace(ttt)

test
## [1] 0.3285851 0.2910063 0.2849646 0.3828346 0.2083211

Now to run the walk:

my_metro <- function(sig, x_0, N, name){
  x <- numeric(N)
  x[1] <- x_0
  reject <- 0                         # index for rejections
  
  for(i in 2:N){
   Y <- rnorm(1, x[i-1], sig)         # incrament via Normal
   U <- runif(1)                      # generate U
   
   if(U <= my_laplace(Y)/my_laplace(x[i-1])){
     x[i] <- Y
   }else{
     reject <- reject + 1 
     x[i] <- x[i-1]
   }
  }
  
  rej_perc <- round(reject/N * 100,   # calc % of rejects
                    digits = 2)
  rslt <- list(walk=data.frame(walk = x, 
                               nm = name,
                               index = c(1:N),
                               stringsAsFactors = F), 
               rejects = rej_perc)
  
  return(rslt)

}

Now I’ll put them together and compare plots.

n = 500
x_0 <- 10
a <- my_metro(0.1,x_0,n, "sig=0.1")
b <- my_metro(0.5,x_0,n, "sig=0.5")
c <- my_metro(1.0,x_0,n, "sig=1")
d <- my_metro(1.5,x_0,n, "sig=1.5")
#rr$rejects
combo <- rbind(a$walk,b$walk,c$walk,d$walk)
ggplot(combo, aes(x = index, y = walk, fill=nm))+
  geom_line(alpha =.5) + 
  facet_wrap(~nm, ncol = 2)

Notes:

The proposal distribution must be chosen so that the generated chain will converge to a stationary distribution - that is: the target distrbution \(f\). Required conditions for the generated chain are irreducibility, postive recurrence, and aperiodicity. A proposal distributuion with the same support distribution will usually satisfy these regularity conditions [pg248].

The Metropols-Hastings Sampler generates the Markov chain \(\{X_0, X_1,...\}\).

  1. Choose proposal distribution \(g(\cdot|X_t)\)
  2. Generated \(X_0\) via \(g\).
  3. Repeat until for \(n\) times.