By implementing the Metropolis-Hastings algorithm, one will be able to generate a sequence of random numbers to understand the use of each parameter.

#setting parameters

N <- 10000
s <- 1
x0 <- 0

#define distribution
f <- function(x) {
  1 / 2 * exp(-abs(x))
}

#Initialize chain
x <- numeric(N)
x[1] <- x0

#start metropolis-hasting coding
for (i in 2:N) {

  x_prop <- rnorm(1, mean = x[i - 1], sd = s)

  r <- f(x_prop) / f(x[i - 1])
#accepting or rejecting value
  u <- runif(1)
  if (u < r) {
    x[i] <- x_prop
  } else {
    x[i] <- x[i - 1]
  }
}

#plotting histogram and kernel density 
hist(x, freq = FALSE, main = "Histogram and Kernel Density Estimate", xlab = "x")
lines(density(x), col = "red")

#shape of graph
curve(f, add = TRUE, col = "blue")

#sample mean and s.d.
mean(x)
## [1] 0.0004132084
sd(x)
## [1] 1.530075
#display graph
plot(x, type = "l", main = "Generated Samples", xlab = "Iteration", ylab = "x")

By running the algorithm multiple times, one will be able to see if the results are consistent with different starting values.

#setting parameter
N <- 2000
J <- 4
s_values <- seq(0.001, 1, length.out = 100)

#define distribution
f <- function(x) {
  1 / 2 * exp(-abs(x))
}

#run metropolis-hasting and calculate rhat
calculate_Rhat <- function(s) {
  #start chain
chains <- matrix(nrow = N, ncol = J)

#run metropolis-hastings for each chain 
for (j in 1:J) {
    x <- numeric(N)
    x[1] <- runif(1)
    for (i in 2:N) {
      x_prop <- rnorm(1, mean = x[i - 1], sd = s)
      r <- f(x_prop) / f(x[i - 1])
      #accept or reject value
      u <- runif(1)
      if (u < r) {
         x[i] <- x_prop
      } else {
        x[i] <- x[i - 1]
      }
    }
    
    chains[, j] <- x
}
#count rhat
 M_j <- colMeans(chains)
  V_j <- apply(chains, 2, var)
  W <- mean(V_j)
  M <- mean(M_j)
  B <- var(M_j)
  Rhat <- sqrt((B + W) / W)

  return(Rhat)
}
Rhat_values <- sapply(s_values, calculate_Rhat)
plot(s_values, Rhat_values, type = "l", main = "Rhat vs. Standard Deviation (s)", 
     xlab = "Standard Deviation (s)", ylab = "Rhat")
abline(h = 1.05, col = "red", lty = 2) 

print(paste("Rhat for s = 0.001:", calculate_Rhat(0.001)))
## [1] "Rhat for s = 0.001: 26.5778676921239"