# STATS 4CI3 January 31, 2019

# Generating gamma random variables (shape>1) using
# exponential candidates in an accept-reject algorithm.

alpha <- 2.5
beta <- 4

# The ratio of densities
foverg <- function(x, alpha, beta, lambda) {
  dgamma(x, alpha, 1/beta)/dexp(x,lambda)
}

# Note that for lambda>1/\beta=0.25 in this case the
# ratio increases without bound

lambda <- 1/(alpha*beta)
# Calculate the upper bound on f(x)/g(x) when 
# lambda=1/(alpha*beta) which minimizes M.
M <- alpha^alpha*exp(1-alpha)/gamma(alpha)

#Generate the candidates
n <- 10000
Y <- -log(runif(n))/lambda
ratio <- foverg(Y, alpha,beta,lambda)
U <- runif(n)
accept <- U<=ratio/M
X <- Y[accept]