#----------------------------------------------------------
#  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)
}


# plot the function
#postscript(file="metropolis_test_function.eps",onefile=FALSE,horizontal=FALSE,height=6,width=7,pointsize=11)
#old.mar<- par("mar")
#par(mar=c(4.1,4.1,0.1,0.1))
curve(fun, -5, 5, col='blue')
#par(mar=old.mar)
#dev.off()

pause()

# sample fun() with Metropolis
n <- 10000
x0 <- 0

# try to change d and tau
x <-  metropolis.annealing(n, x0, fun, prop, d=0.1, tau=2000)



#postscript(file="metropolis_test_function_hist_d1.eps",onefile=FALSE,horizontal=FALSE,height=6,width=7,pointsize=11)
#old.mar<- par("mar")
#par(mar=c(4.1,4.1,0.1,0.1))
hist(x, nc=100, col='cyan', main=NULL)
#par(mar=old.mar)
#dev.off()

pause()

plot(x, ty='l')