Attention! The parameter in (2.2.1) is scale parameter, thus the code to sample is
## sample from Ga(a, beta)functionrgamma_int(a::Int, beta) u =rand(a)return(-1.0/ beta *sum(log.(u)))end
## sample from Ga(a, beta)functionrgamma_int(a::Int, beta) u =rand(a)return(-1.0* beta *sum(log.(u)))end## density of Ga(alpha, beta)include("func_gamma.jl")functiondgamma(x, alpha, beta)return(beta^alpha /lanczos_gamma(alpha) * x^(alpha-1) *exp(-1*beta*x))end## accept-reject algorithmfunctionrgamma(alpha =1.5, beta =2.1) a =Int(floor(alpha)) b = a * beta / alpha M =exp(a * (log(a) -1) - alpha * (log(alpha) -1))whiletrue x =rgamma_int(a, b) u =rand() cutpoint =dgamma(x, alpha, beta)/(M*dgamma(x, a, b))if u <= cutpointreturn xendendend# exampleprintln(rgamma())
Firstly, let's introduce how to generate Ga(a,β) where a is an integral (Robert and Casella (2013)).
Then when considering the general Ga(α,β), we can use Ga(a,b) as instrumental distribution in Accept-Rejection algorithm as mentioned in Robert and Casella (2013):