Special Distribution

Search…
Gamma distribution
Firstly, let's introduce how to generate Ga(a,β){\mathcal G}a(a,\beta)Ga(a,β)
 where aaa
 is an integral (Robert and Casella (2013)).
Attention! The parameter in (2.2.1) is scale parameter, thus the code to sample is
1
## sample from Ga(a, beta)
2
function rgamma_int(a::Int, beta)
3
    u = rand(a)
4
    return(-1.0 / beta * sum(log.(u)))
5
end
Copied!
Then when considering the general Ga(α,β){\mathcal G}a(\alpha,\beta)Ga(α,β)
, we can use Ga(a,b){G}a(a, b)Ga(a,b)
 as instrumental distribution in Accept-Rejection algorithm as mentioned in Robert and Casella (2013):
1
## sample from Ga(a, beta)
2
function rgamma_int(a::Int, beta)
3
    u = rand(a)
4
    return(-1.0 * beta * sum(log.(u)))
5
end
6
​
7
## density of Ga(alpha, beta)
8
include("func_gamma.jl")
9
function dgamma(x, alpha, beta)
10
    return(beta^alpha / lanczos_gamma(alpha) * x^(alpha-1) * exp(-1*beta*x))
11
end
12
​
13
## accept-reject algorithm
14
function rgamma(alpha = 1.5, beta = 2.1)
15
    a = Int(floor(alpha))
16
    b = a * beta / alpha
17
    M = exp(a * (log(a) - 1) - alpha * (log(alpha) - 1))
18
    while true
19
        x = rgamma_int(a, b)
20
        u = rand()
21
        cutpoint = dgamma(x, alpha, beta)/(M*dgamma(x, a, b))
22
        if u <= cutpoint
23
            return x
24
        end
25
    end
26
end
27
​
28
# example
29
println(rgamma())
Copied!
where func_gamma.jl is to calculate gamma function via Lanczos Approximation because there is not gamma function in Julia.
If β=1\beta=1β=1
, we can do some further calculations,
and would get more simpler expression,
Generate R.V.
Copulas
Last modified 3yr ago
Copy link