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