set.seed(325)
# Use the rejection sampling method to generate a random sample of size 10000 from the Gamma(2,2) distribution f(y).
# Draw a target function f(y) to obtain the idea of proposal distribution
f_y <- function(y){
 y^(2-1)*exp(-2*y)
}

y <- seq(0.01, 4, 0.01)
plot(y, f_y(y), type ="l", xlab = "", ylab = "")
# Use chi-sq (df = 3) distribution as a proposal distribution g(y).
y <- seq(0.01, 4, 0.01)
plot(y, f_y(y), type ="l", xlab = "", ylab = "")
lines(y, dchisq(y, df = 3), col = "red")
# Draw a f(y)/g(y) plot to choose c
y <- seq(0.01, 4, 0.01)
plot(y, f_y(y)/dchisq(y,df = 3), type ="l", xlab = "", ylab = "")
max(f_y(y)/dchisq(y, df = 3))
# choose c = 0.88

# Rejection Sampling
n <- 10000
k <- 0 # counter for accepted
j <- 0 # counter for iterations
x <- rep(0, n)

while(k < n){
  y <- rchisq(1, df = 3) #random variate from g(y) which is Chisq(df = 3)
  u <- runif(1)
  j <- j + 1
  if (u < f_y(y)/(0.88*dchisq(y,df = 3))){
    #we accept y
    k <- k + 1
    x[k] <- y
  }
}

mean(x)

# Compare empirical and theoretical percentiles
p <- seq(0.1, 0.9, 0.01)
Qhat <- quantile(x, p) #quantile of sample
Q <- qnorm(p) #theoretical quantiles
plot(Q, Qhat)

round(rbind(Q, Qhat),2)