#----------------------------------------------------------
#  Example of Metropolis algoritm in 1D
#
#  GdA 5/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
prop <- function(x) {d=10; x + runif(1, -d, d)}

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

# Metropolis algorithm
metropolis <- function(n, x0, fun, prop) {
  x = rep(0,n)
  x[1] = x0
  for (i in 2:n) {
    x.p <- prop(x[i-1])
    A <- fun(x.p)/fun(x[i-1])
    x[i] <- ifelse (runif(1) <= A, x.p, x[i-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

x <-  metropolis(n, 0, fun, prop)

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

#postscript(file="metropolis_test_function_plot_d10.eps",onefile=FALSE,horizontal=FALSE,height=3,width=7,pointsize=11)
#old.mar<- par("mar")
#par(mar=c(4.1,4.1,0.1,0.1))
plot(x[1:400], col='blue', ty='b', pch=19, cex=0.4)
#par(mar=old.mar)
#dev.off()