#-------------------------------------------------------
# simula una funzione strana con x continua fra 0 e 1
# usando il metodo hit/miss
#
# GdA,  PED 25/11/2010
#-------------------------------------------------------

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


my.func <- function(x) {
   # att. non e' normalizzata, ma non ci interessa...
   return (abs(sin(x*10)*(x-0.5)^2)^0.2)
}

# controlliamo la funzione
x <- seq(0, 1, by=0.001)
plot(x, my.func(x), ty='l')
pausa()

# massimo (con fattore di tolleranza)
fmax = max( my.func(x) ) * 1.05

# funzione che estrae un numero random alla volta secondo my.func()
my.rfunc <- function() {
  ok = FALSE
  while (!ok) {
    x = runif(1)
    r = my.func(x) / fmax
        y = runif(1)
    if ( y <= r ) {
      ok = TRUE
    }
  }
  return(x)
}

# e ora generiamo n.evts eventi
n.evts = 100000
n = 10
p = 1/2

my.x = rep(0,n.evts)
for (i in 1:n.evts) {
  my.x[i] = my.rfunc()
}
hist(my.x,nc=200)


pausa()

#--------------------------------------------------------
# variazione sul tema: innanzitutto conserviamo i punti
# 'dentro' e quelli 'fuori'

n.evts = 30000 # cominciamo con pochi per vederli meglio

x.in  = rep(0, n.evts)
y.in  = rep(0, n.evts)
x.out  = rep(0, n.evts)
y.out  = rep(0, n.evts)

n.in  = 0
n.out = 0
for (i in 1:n.evts) {
    x = runif(1)
    y = runif(1)*fmax
    if ( y <= my.func(x) ) {
      n.in = n.in + 1
      x.in[n.in] = x
      y.in[n.in] = y 
    } else {
      n.out = n.out + 1
      x.out[n.out] = x
      y.out[n.out] = y
    }
}
x.in  = x.in[1:n.in]
y.in  = y.in[1:n.in]
x.out = x.out[1:n.out]
y.out = y.out[1:n.out]

x <- seq(0, 1, by=0.001)
plot(x, my.func(x), ty='l', xlim=c(0,1), ylim=c(0, fmax))
points(x.out, y.out, col='red', cex=2, pch='.')
points(x.in, y.in, col='black', cex=2, pch='.')

#-------  ------------------------------------------------------
save.plot = TRUE  # mettere FALSE se non si vuole salvare il plot
if (save.plot){
  nome.file <- "funzione_strana_hit-miss.png"
  png(filename=nome.file, width=1200, height=800)
  plot(x, my.func(x), ty='l', xlim=c(0,1), ylim=c(0, fmax))
  points(x.out, y.out, col='red', cex=2, pch='.')
  points(x.in, y.in, col='black', cex=2, pch='.')
  dev.off()
}
#---------------------------------------------------------------

# a questo punto, osservando lo 'scatter-plot', IDEA!
area.rettangolo = fmax*1
integrale.my.func = n.in / n.evts * area.rettangolo