#----------------------------------------------------------
#  Example of Metropolis algoritm in 1D
#
#  GdA 5/4/2016 (rev 10/5/2020)
#-----------------------------------------------------------

# function to sample 
fun <- function(x) (0.5*dnorm(x, -2.5, 0.8) + dnorm(x, 3.1, 1) ) * exp(-abs(x)/0.43)

# proposal function
prop <- function(x, d=1)  x + runif(1, -d, d) 

# pause
pause <- function() { cat ("\n >> Press Enter to continue\n"); scan() }

# Metropolis algorithm
metropolis <- function(n, x0, fun, prop, ...) {    # NOTE this R feature !
  x = rep(0,n)
  x[1] = x0
  for (i in 2:n) {
    x.p <- prop(x[i-1], d)
    A <- fun(x.p)/fun(x[i-1])
    x[i] <- ifelse (runif(1) <= A, x.p, x[i-1])
  }
  return(x)
}


# plot the function
curve(fun, -5, 5, col='blue')
pause()

# sample fun() with Metropolis
n <- 10000
x0 <- 0 #  -10
d = 1

x <-  metropolis(n, x0, fun, prop, d)

hist(x, nc=100, col='cyan', main=NULL)
pause()

plot(x[1:n], col='blue', ty='b', pch=19, cex=0.4)