Monte-Carlo
  • Introduction
  • AIS
  • Generate R.V.
    • Special Distribution
    • Copulas
    • Minimum of Two Exponential
  • Gibbs
    • Comparing two groups
    • Linear Regression
    • Simulation of Exp-Abs-xy
  • Markov Chain
  • MC Approximation
  • MC Integration
    • Rao-Blackwellization
  • MC Optimization
  • MCMC
    • MCMC diagnostics
  • Metropolis-Hastings
    • Metropolis
    • Independent MH
    • Random Walk MH
    • ARMS MH
  • PBMCMC
  • RJMCMC
  • Diagnosing Convergence
  • SMCTC
  • Tempering
    • Parallel Tempering
  • Misc
    • R vs. Julia
  • References
Powered by GitBook
On this page
  1. Generate R.V.

Special Distribution

PreviousGenerate R.V.NextCopulas

Last updated 6 years ago

Gamma distribution

Firstly, let's introduce how to generate Ga(a,β){\mathcal G}a(a,\beta)Ga(a,β) where aaa is an integral ().

Attention! The parameter in (2.2.1) is scale parameter, thus the code to sample is

## sample from Ga(a, beta)
function rgamma_int(a::Int, beta)
    u = rand(a)
    return(-1.0 / beta * sum(log.(u)))
end
## sample from Ga(a, beta)
function rgamma_int(a::Int, beta)
    u = rand(a)
    return(-1.0 * beta * sum(log.(u)))
end

## density of Ga(alpha, beta)
include("func_gamma.jl")
function dgamma(x, alpha, beta)
    return(beta^alpha / lanczos_gamma(alpha) * x^(alpha-1) * exp(-1*beta*x))
end

## accept-reject algorithm
function rgamma(alpha = 1.5, beta = 2.1)
    a = Int(floor(alpha))
    b = a * beta / alpha
    M = exp(a * (log(a) - 1) - alpha * (log(alpha) - 1))
    while true
        x = rgamma_int(a, b)
        u = rand()
        cutpoint = dgamma(x, alpha, beta)/(M*dgamma(x, a, b))
        if u <= cutpoint
            return x
        end
    end
end

# example
println(rgamma())

If β=1\beta=1β=1, we can do some further calculations,

and would get more simpler expression,

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 :

where is to calculate gamma function via because there is not gamma function in Julia.

Robert and Casella (2013)
func_gamma.jl
Lanczos Approximation
Robert and Casella (2013)