#
#  Example of rejection sampling (Section 9.5)
#
n <- 1000
alpha <- 2
beta <- 4
cc <- exp((alpha-1)*log(alpha-1)+(beta-1)*log(beta-1)-
          (alpha+beta-2)*log(alpha+beta-2))
theormean <- alpha/(alpha+beta)
theorvar <- alpha*beta/
            ((alpha+beta)*(alpha+beta)*(alpha+beta+2))
mean <- 0
ss <- 0
for (i in 1:n)
  {
     cont <- TRUE
     while (cont)
     {
        y <- runif(1)
        u <- runif(1)
        if (u <= exp((alpha-1)*log(y)+(beta-1)*log(1-y)))
        {
           x <- y
           mean <- mean + x/n
           ss <- ss + x*x
           cont <- FALSE
        }
     }
  }
var <- (ss-n*mean*mean)/(n-1)
cat("\n")
cat(" Alpha =",alpha,"Beta =",beta,"; Mean =",mean,"Variance =",var,"\n")
cat(" Theoretical values         ",theormean,"and       ",theorvar,"\n")
cat(" Ratios                     ",mean/theormean,"and       ",var/theorvar)
cat("\n\n")