# esempio introduttivo alle MCMC
#
# questo script può essere eseguito con il comando
# source("sassi_a_caso.R")
#

# funzione di comodo se si hanno molti grafici
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}


# "sassi lanciati a caso"

# vettori con le coordinate
n = 10000
x <- y <- rep(0, n)

# distanza massima del lancio
d.max <- 2

# limiti del quadrato (o, eventualmente, rettangolo)
xl <- c(-1,1)
yl <- c(-1,1)

# posizione iniziale arbitraria, ad esempio
x[1] <- 0
y[1] <- 0

for (i in 2:n) {
  # si lancia "a caso"
  phi <- runif(1, max=2*pi)
  d   <- runif(1, max=d.max)
  x.new <- x[i-1] + d*cos(phi)
  y.new <- y[i-1] + d*sin(phi)
  if((x.new >= xl[1]) & (x.new <= xl[2]) &
     (y.new >= yl[1]) & (y.new <= yl[2]) ) {
    x[i] <- x.new
    y[i] <- y.new
  } else {
    x[i] <- x[i-1]
    y[i] <- y[i-1] 
  } 
}
# plottiamo i punti, aggiungendo del rumore per dare l'idea
# dei mucchietti
plot(x+rnorm(n,0,0.005), y+rnorm(n,0,0.005), pch=19, cex=0.1, asp=1,
     xlab='x', ylab='y', main=sprintf("distanza max  %.1f",d.max))

# indichiamo i punti nel cerchio (attenzione ai mucchietti
# artificiosamente sparpagliati!)
symbols(0, 0, circles=1, bg=rgb(1,0.,0., 0.05), add=TRUE, inches=FALSE)

# condizione di punti nel cerchio
dentro <- ( x^2 + y^2 ) <= 1

# stima pi greco
pi.mc <- sum(dentro) / n * 4

cat(sprintf("pi greco stimato: %f\n",pi.mc))

# pausa
pausa()

# "storia" dei valori di x: -> si studi la dipendenza da d.max
# limitiamo il nr di valori da plottare per osservare meglio i dettagli
nmax=1000
plot(x[1:nmax], ty='l')

# "storia" dei valori di x: -> si studi la dipendenza da d.max
# limitiamo il nr di valori da plottare per osservare meglio i dettagli
pausa()
nmax=100
plot(x[1:nmax], y[1:nmax], ty='l', col='blue', asp=1,
     xlab='x', ylab='y', main=sprintf("distanza max  %.1f",d.max))

# altro modo di mostrare la storia della "catena"
pausa()
nmax=100
arrows(x[1:(nmax-1)], y[1:(nmax-1)], x[2:nmax], y[2:nmax],
       col='red', len=0.1, lwd=1.3)
points(x[1], y[1])
points(x[nmax], y[nmax])