# Task 1

sampleExpo <- function(lam, n) {
  u <- runif(n)
  x <- -log(u)/lam
  x
}

# Task 2

lam <- 2
x <- sampleExpo(lam, 1000)

hist(x, freq=FALSE, col="grey")
xi <- seq(0, 5, length.out=100)
lines(xi, lam*exp(-lam*xi), lwd=2)


# Task 3

sampleGammaInt <- function(k, lam) {
  sum(sampleExpo(lam, k))
}

# Task 4

alpha <- 5.7
k <- floor(alpha)
lam <- 2

M <- dgamma(alpha-k,alpha,lam)/dgamma(alpha-k,k,lam-1)

xi <- seq(0, 10, length.out=100)

plot(xi, M*dgamma(xi,k,lam-1), type="l", lty=2, col="red",
     xlab="x", ylab="Density")
lines(xi, dgamma(xi,alpha,lam), col="darkblue")
legend("topright", col=c("red","darkblue"),lty=c(2,1),
       legend=c("Mg(x) [scaled instrumental]","f(x) [target]"))

       
# Task 5

sampleGamma <- function(alpha, lam) {
  k <- floor(alpha)
  M <- dgamma(alpha-k,alpha,lam)/dgamma(alpha-k,k,lam-1)
  while (TRUE) {
    x <- sampleGammaInt(k,lam-1)
    prob.of.accept <- dgamma(x,alpha,lam)/(M*dgamma(x,k,lam-1))
    if (runif(1)<prob.of.accept)
      return(x)
  }
}

alpha <- 5.7
lam <- 2
x <- numeric(1000)
for (i in 1:1000)
  x[i] <- sampleGamma(alpha,lam)

# Task 6

hist(x, freq=FALSE, col="grey")
xi <- seq(0, 10, length.out=100)
lines(xi, dgamma(xi, alpha,lam), lwd=2)