#-------------------------------------------------------
# test metropolis + simulated annealing
#
# GdA 16 marzo 2012
#-------------------------------------------------------

library(compiler)

# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

# funzione da riprodurre 
fun <- function(x) dnorm(x) +  dnorm(x, 1.0, 0.02) +  dnorm(x, -3.0, 0.1)

# implementazione di Metropolis, incluso annealing
# (funzione di proposta normale di sigma s.p)
s.metropolis <- function(n, y0, s.p, alpha) {
  # alpha -> infinito: no annealing
  y = rep(0,n)
  y[1] = y0
  for (i in 2:n) {
    prova = TRUE
    while(prova) {
      yt = y[i-1] + rnorm(1, 0, s.p)
      A = fun(yt) / fun(y[i-1])
      A = A^(1 + i/alpha)
      if (runif(1) <= A) {
          y[i] = yt
          break
      }
    }  
  }
  return(y)
}

# Compiliamo per velocizzare
metropolis  <- cmpfun(s.metropolis)

# mostra la funzione
x <- seq(-6,3,by=0.01)
plot(x,fun(x),ty='l')
pausa()

# chiama Metropolis senza ennealing
cat(sprintf("\n Metropolis semplice\n\n"))
y0 = 0
s.p = 1
n = 10000
alpha = Inf
y <-  metropolis(n, y0, s.p, alpha)

hist(y, nc=200, xlim=c(-6,3))

pausa()

# chiama Metropolis con alpha=10000
alpha = 10000
cat(sprintf("\nMetropolis + Annealing con alpha=%d \n", alpha))
y <-  metropolis(n, y0, s.p, alpha)

hist(y[1:1000], nc=100, xlim=c(-6,3))
pausa()
hist(y[9000:10000], nc=100, xlim=c(-6,3))
pausa()
plot(y,cex=0.2)

cat(sprintf("media ultimi 1000 eventi %f\n", mean(y[9001:10000])))
cat(sprintf("std    ultimi 1000 eventi %f\n", sd(y[9000:10000])))


pausa()
# chiama Metropolis con alpha=1000
alpha = 1000
cat(sprintf("\nMetropolis + Annealing con alpha=%d \n", alpha))
y <-  metropolis(n, y0, s.p, alpha)

hist(y[1:1000], nc=100,  xlim=c(-6,3))
pausa()
hist(y[9000:10000], nc=100, , xlim=c(-6,3))
pausa()
plot(y,cex=0.2)

cat(sprintf("media ultimi 1000 eventi %f\n", mean(y[9001:10000])))
cat(sprintf("std    ultimi 1000 eventi %f\n", sd(y[9000:10000])))

pausa()
# chiama Metropolis con alpha=1000
alpha = 1000
s.p = 0.5
cat(sprintf("\nMetropolis + Annealing con alpha=%d e s.p=%.1f\n", alpha, s.p))
y <-  metropolis(n, y0, s.p, alpha)

hist(y[1:1000], nc=100,  xlim=c(-6,3))
pausa()
hist(y[9000:10000], nc=100, , xlim=c(-6,3))
pausa()
plot(y,cex=0.2)

cat(sprintf("media ultimi 1000 eventi %f\n", mean(y[9001:10000])))
cat(sprintf("std    ultimi 1000 eventi %f\n", sd(y[9000:10000])))