# 
# Simulation of beta distribution by rejection sampling
# 
n <- 1000
theta <- rep(0,n)
alpha <- 3
beta <- 2
f <- function(x) { x^(alpha-1)*(1-x)^(beta-1) }
c <- (alpha-1)^(alpha-1)*(beta-1)^(beta-1)
h <- function() { runif(1) }
for (i in 1:n){
  defined <- FALSE
  while (!defined){
    Y <- h()
    U <- runif(1)
    if (U <= f(Y)/c*1){
      theta[i] <- Y
      defined <- TRUE
    }
  }
}
d <- density(theta, bw = "sj")
simname <- paste("Simulation of beta(",alpha,",",beta,")")
plot(d,xlab="",ylab="",xlim=c(0,1),ylim=c(0,2),main=simname)
par(new=T)
f <- function(x) dbeta(x,alpha,beta)
curve(f,0,1,xlim=c(0,1),ylim=c(0,2),xlab="",ylab="",main="")