#---------------------------------------------
# simula binomiale a partire da random uniform
#
# GdA, PED 25/11/2010
#---------------------------------------------

my.rbinom <- function(n, p, fmax) {
  # genera numeri secondo binomiale (uno alla volta!)
  # usando solo runif()

  ok = FALSE

  while (!ok) {
    # genera x uniforme fra 0 e n
    x = floor( runif(1)*(n+1) )

    # ripesa gli eventi secondo la distribuzione binomiale
    r = dbinom(x, n, p) / fmax
    y = runif(1)
    if ( y <= r ) {
      ok = TRUE
    }
  }
  return(x)
}

n.evts = 10000
n = 10
p = 1/2

# calcola fmax
fmax = max( dbinom(c(0:n), n, p) )
n.binom = rep(0,n.evts)

for (i in 1:n.evts) {
  n.binom[i] = my.rbinom(n, p, fmax)
}

h <- hist(n.binom+0.1, nc=11) 

titolo = sprintf("binomiale (n=%d,p=%.1f) simulata con %d eventi", n, p, n.evts)
barplot( h$counts, names=h$breaks[1:(n+1)], main=titolo)
# provare a commentare il comando precedente per vedere cosa fa hist

# inoltre, per capire cosa viene immagazzinato in h
str(h)