my.rbinom <- function(n, p) { # MOLTO INEFFICIENTE!! fxm <- max(dbinom(0:n, n, p)) max.prove <- 100 prova = 1 while(prova < max.prove) { x <- sample(0:n)[1] fx <- dbinom(x, n, p) rf <- runif(1)*fxm if (rf < fx) break() } return(x) } # test test = TRUE if(test) { N = 10000 x <- rep(0, N) for (i in 1:N) x[i] <- my.rbinom(10, 0.5) plot(table(x)) }