make.beta.log.density <- function(alpha=1,beta=1) {
  function(x) {
    log(dbeta(x,alpha,beta))
  }
}

unif.log.density <- function(x) { 0 }
unif.sampler <- function() { runif(1) }

a1 <- 2
a2 <- 5
mode <- (a1-1)/(a1+a2-2)  # if a1 > 1 and a2 > 1
beta.log.density <- make.beta.log.density(a1,a2)


x <- seq(0,1,by=0.01)
plot(x,rep(1,length(x)),type="l",ylim=c(0,2.5))
lines(x,dbeta(x,a1,a2),col="blue")
lines(x,rep(1,length(x))*exp(beta.log.density(mode)),col="red")




source("rejection-sampler.R")
d <- rejection.sampler(beta.log.density,unif.log.density,unif.sampler,1/exp(beta.log.density(mode)),10000)

hist(d,prob=TRUE,xlab="x",main=paste("Samples for a beta(",a1,",",a2,") distribution",sep=""))
lines(density(d),col="red")

x <- seq(0,1,by=0.01)
lines(x,dbeta(x,a1,a2),col="blue")