# estrazioni random secondo una generica funzione
# non necessariamente normalizzata, in un certo intervallo,
# purché non negativa (il programma taglia la parte negativa!)

# definizione della funzione
# (per ora senza parametri passati dall'esterno: per altre funzioni
# è sufficiente definire altre funzioni, ad es. fun2, etc.)
fun1 <- function(x) {
  x^2 + x -abs(sin(x))^0.5 + 1.1
}

genera.hm <- function(n, fun, xmin=0, xmax=1) {
  # cerchiamo brutalmente il massimo di fun (con 1% di margine!)
  f.max <- max( fun(seq(xmin, xmax, len=10000)) ) * 1.01
  x <- runif(n, xmin, xmax)
  y <- runif(n) * f.max
  xr <- x[y <= fun(x)]
  return(xr)
}

# per visualizzare la funzione, scommentare
# x<- seq(-1,1,len=100); plot(x, fun1(x), ty='l')

# esempio di uso
n.evts <- 100000
xmin <- -1;
xmax <- 1

n <- n.evts
xr <- genera.hm(n, fun1, xmin, xmax)

# contiamo il nr di eventi generati
nr <- length(xr)
# se non sono tutti (efficienza dell'algoritmo!) andiamo avanti
eff <- nr/n.evts
# come sottoprodotto ci valutiamo l'integrale di max(0, fun(x))
f.max <- max( fun1(seq(xmin, xmax, len=10000)) ) * 1.01 # (*)
area.rett <- (xmax-xmin)*f.max
area.fun <- area.rett*eff

n1 <- (n.evts - nr) / eff  # stima di prove da effettuare (**)
if(n1 > 0) xr <- c(xr, genera.hm(n1, fun1, xmin, xmax) )

hist(xr, nc=100, main=sprintf("nr evts generati/richiesti: %d/%d\n [efficienza algoritmo: %.3f;  richiesti extra %.0f evts]\n [integrale = %.3f]",
                     length(xr), n.evts, eff, n1, area.fun) )

# NOTE:
#  (*) Ovviemanente non è una bella cosa che si usi il fattore
#      di "tolleranza" 1.01 usato all'interno di genera.hm()
#      
# (**) Il procedimento è approssimativo e al secondo giro
#      dà un numero di eventi molto prossimo a quello desiderato
#
#      => ciascuno si può sbizarrire a migliorare il codice