9.4 Discussion

library(tidyverse)
library(scales)

9.4 Implement a random walk Metropolis sampler for generating the standard Laplace distribution (see Exercise 3.2). For the increment, simulate from a normal distribution. Compare the chains generated when different variances are used for the proposal distribution. Also, compute the acceptance rates of each chain.

The standard Laplace distribution which has density \(f(x) = \frac{1}{2}e^{-|x|},x \in I\!R\)

laplace_dist <- function(x){return((1/2)*exp(-abs(x)))}

RW.Metro.Laplace <- function(sigma, x.0, N) { 
                      x <- numeric(N) 
                      x[1] <- x.0 
                      u <- runif(N) 
                      k <- 0 
                      for (i in 2:N) { 
                        y <- rnorm(1, x[i-1], sigma) #increment
                        if (u[i] <= (laplace_dist(y) / laplace_dist(x[i-1]))) #proposal 
                            x[i] <- y 
                        else { 
                          x[i] <- x[i-1] 
                             k <- k + 1 } 
                      } 
                      return(list(x = x, k = k)) 
                    }

N <- 200
sigma <- c(1/3, 2/3, 1, 3)
x.0 <- 100
results <- NULL 

for (i in 1:length(sigma)){
  a <- c(.05, seq(.1, .9, .1), .95)
  results[[i]] <- c(sigma_used = sigma[i], N = N, result = 
                    list(RW.Metro.Laplace(sigma = sigma[i], x.0 = x.0, N = N))) %>% 
                  c(., accepted.perc = scales::percent(1-unlist(.[[3]]['k'])/N),
                       quantile = list(quantile(unlist(.[[3]]['x']), a)))
}

results %>% map(., function(x) c(x[1], x[4], x[5]))

## [[1]]
## [[1]]$sigma_used
## [1] 0.3333333
## 
## [[1]]$accepted.perc
## [1] "92%"
## 
## [[1]]$quantile
##        5%       10%       20%       30%       40%       50%       60% 
##  91.72473  91.91172  92.46968  92.80247  93.19715  93.70292  94.68265 
##       70%       80%       90%       95% 
##  96.83940  98.19781  99.64213 100.17239 
## 
## 
## [[2]]
## [[2]]$sigma_used
## [1] 0.6666667
## 
## [[2]]$accepted.perc
## [1] "78%"
## 
## [[2]]$quantile
##       5%      10%      20%      30%      40%      50%      60%      70% 
## 74.03986 74.24794 74.93441 76.01350 78.08940 82.41742 85.24023 90.09921 
##      80%      90%      95% 
## 91.63237 94.70841 96.88536 
## 
## 
## [[3]]
## [[3]]$sigma_used
## [1] 1
## 
## [[3]]$accepted.perc
## [1] "79.5%"
## 
## [[3]]$quantile
##       5%      10%      20%      30%      40%      50%      60%      70% 
## 54.26417 56.22853 63.09038 70.25908 76.53846 82.16091 87.13873 90.29959 
##      80%      90%      95% 
## 94.41561 98.18941 99.04662 
## 
## 
## [[4]]
## [[4]]$sigma_used
## [1] 3
## 
## [[4]]$accepted.perc
## [1] "54%"
## 
## [[4]]$quantile
##         5%        10%        20%        30%        40%        50% 
## -1.6309785 -0.9127513 -0.6537852  0.1531675  0.7101430  2.0380109 
##        60%        70%        80%        90%        95% 
##  3.5040530 31.3634156 44.1652177 82.5001811 84.8880256

Rizzo; Maria L.. Statistical Computing with R (Chapman & Hall/CRC The R Series) (Page 254). Taylor and Francis CRC ebook account. Kindle Edition.