#--------------------------------------------
# pi.mcmc.R: esempio minimalista di MCMC
# per campionare uniformemente il quadrato 
# {(-1<x<+1),(-1<y<+1}
#
# GdA, 26/2/09
#---------------------------------------------
 
#-- funzioni ----------------------------------
#-- (fare copia/incolla o mettere in uno script)
proponi.unif <- function(xy.old, d)
{
  phi <- runif(1)*2*pi
  r   <- runif(1)*d
  dx   <- r * cos(phi)
  dy   <- r * sin(phi)
  return ( c(xy.old[1]+dx, xy.old[2]+dy) )
}
 
campiona.unif <- function(n, xy0=c(0,0), raggio=1)
{
  x <- numeric()
  y <- numeric()
  x[1] <- xy0[1]
  y[1] <- xy0[2]
  for (i in c(2:n)) {
    New.xy <- proponi.unif( c(x[i-1], y[i-1]), raggio)
    if( (abs(New.xy[1]) <= 1) &&
       (abs(New.xy[2]) <= 1) ) {
      x[i] <- New.xy[1]
      y[i] <- New.xy[2]
    } else {
      x[i] <- x[i-1]
      y[i] <- y[i-1]
    }
  }
  return( cbind(x,y) )
}
 
#utilizzzazione (sessione interattiva) -----------------
xy <- campiona.unif(10000)
x <- xy[,1]
y <- xy[,2]
plot(x, y)
plot(x, ty='l')
#pi greco:
length(x[(x^2+y^2) <= 1]) /length(x) *4
 
# cambiamo raggio max di proposta:
# raggio grande:
xy.r10 <- campiona.unif(10000, r=10)
x.r10 <- xy.r10[,1]
y.r10 <- xy.r10[,2]
plot(x.r10, y.r10)
plot(x.r10, ty='l')
#pi greco:
length(x.r10[(x.r10^2+y.r10^2) <= 1]) /length(x.r10) *4
 
# raggio piccolo:
xy.r01 <- campiona.unif(10000, r=0.1)
x.r01 <- xy.r01[,1]
y.r01 <- xy.r01[,2]
plot(x.r01, y.r01)
plot(x.r01, ty='l')
#pi greco:
length(x.r01[(x.r01^2+y.r01^2) <= 1]) /length(x.r01) *4
 
 
# prove suggerite:
# - giocare con funzione di proposta, anche non simmetrica!
# - dipendenza da scelta punto iniziale (e come eventualmente
#   scartare la 'testa')
# - ottimizzazione efficienza e precisione su pi greco
# - graficare (con pochi passi) il moto nel piano xy