#----------------------------------------------------------
#  Example of Metropolis algoritm in 1D
#  Variant with simulated annealing in order to find the mode
#
#  GdA 8/4/2016
#-----------------------------------------------------------

# 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 [now we pass d as parameter]
prop <- function(x, d=1) { x + runif(1, -d, d)}

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

# Metropolis algorithm
metropolis.annealing <- function(n, x0, fun, prop, d=1, tau=Inf) {
  # d is the starting width of the proposal
  # tau is the "time constant" of the Temperature, i.e T=1*exp[-t/tau]   
  x = rep(0,n)
  x[1] = x0
  for (t in 2:n) {            # we use now the index t (as time)
    x.p <- prop(x[t-1], d)
    T <- exp(-t/tau)                             #  Temperature
    A <- ( fun(x.p)/fun(x[t-1]) )^(1/T)          # cooling
    x[t] <- ifelse (runif(1) <= A, x.p, x[t-1])
  }
  return(x)
}

curve(fun, -5, 5, col='blue')

pause()

n <- 10000
x0 <- 0

x <-  metropolis.annealing(n, x0, fun, prop, d=0.1, tau=2000)

# -> try to play with d and tau and see how performances change.

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

pause()

plot(x, ty='l')