# STATS 4CI3/6CI3 February 4

# I meant to do this example of the Accept-Reject algorithm for
# generating Beta random variables provided both parameters are
# not less than 1.
#
# Please review this code and if there are questions let me know 
# in class on Wednesday.
#

rbeta.AR <- function(n, alpha, beta) {
  #
  # Accept-reject sampling for the beta distribution
  # with alpha>=1 and beta>=1
  #
  if ((alpha<1)||(beta<1))
    stop("Both parameters must be at least 1")
  if (alpha+beta>2)
    x0 <- (alpha-1)/(alpha+beta-2) 
  else x0 <- 0.5   
  M <- dbeta(x0, alpha, beta)
  Y <- runif(M*n)
  U <- runif(M*n)
  accept <- U<= dbeta(Y, alpha, beta)/M
  X <- Y[accept]
  n2 <- n-sum(accept)
  while (n2 > 0) {
    Y1 <- runif(M*n2)
    U1 <- runif(M*n2)
    accept1 <- U1<=dbeta(Y1, alpha, beta)/M
    X <- c(X, Y1[accept1])
    n2 <- n-length(X)
  }
  X[1:n]
}

# Two Monte Carlo estimators of pi considered in class
N <- 1000000

X <- runif(N)
Y <- runif(N)
p.hat <- mean(X^2+Y^2<1)
pi.hat1 <- 4*p.hat
se.1 <- 4*sqrt(p.hat*(1-p.hat)/N)

Eg <- mean(sqrt(1-X^2))
pi.hat2 <- 4*Eg
se.2 <- 4*sqrt((mean(1-X^2)-Eg^2)/N)