Simple example of sampling using Metropolis algorithm

metropolis-algorithm-example.R

set.seed(123)

x <- rnorm(7)

dx <- density(x)
plot(dx)

f <- approxfun(dx, yleft = 0, yright = 0)

integrate(f, lower = -6, upper = 6)

xx <- seq(-6, 6, by = 0.001)
plot(xx, f(xx), type = "l")
y <- sample(xx, 1e5, prob = f(xx), replace = TRUE)

hist(y, freq = FALSE, 100)
lines(xx, f(xx), col = "red")


R <- 1e5
y <- numeric(R)
tmp <- 0

for (i in 1:R) {
  prop <- tmp + rnorm(1)
  u <- runif(1)
  A <- dnorm(tmp - prop, log = TRUE) - dnorm(prop - tmp, log = TRUE)
  A <- A + log(f(prop)) - log(f(tmp))
  if (log(u) < A) {
    tmp <- prop
  }
  y[i] <- tmp
}

lines(xx, f(xx), col = "blue", lty = 2)