## sample from Ga(a, beta)
function rgamma_int(a::Int, beta)
return(-1.0 * beta * sum(log.(u)))
## density of Ga(alpha, beta)
function dgamma(x, alpha, beta)
return(beta^alpha / lanczos_gamma(alpha) * x^(alpha-1) * exp(-1*beta*x))
## accept-reject algorithm
function rgamma(alpha = 1.5, beta = 2.1)
M = exp(a * (log(a) - 1) - alpha * (log(alpha) - 1))
cutpoint = dgamma(x, alpha, beta)/(M*dgamma(x, a, b))