Generalised Gamma Distribution

The Generalised Gamma Distribution

The Generalised Gamma (GG) distribution [1] is a three-parameter distribution with support on ${\mathbb R}_+$. The corresponding hazard function can accommodate bathtub, unimodal and monotone (increasing and decreasing) hazard shapes. The GG distribution has become popular in survival analysis due to its flexibility. Other flexible distributions that can account for these hazard shapes are discussed in @rubio:2021 and @jones:2015.

Probability Density Function

The pdf of the GG distribution is

\[f(t;\theta,\kappa,\delta) = \dfrac{\delta}{\Gamma\left(\frac{\kappa}{\delta}\right)\theta^\kappa} t^{{\kappa-1}}e^{{-\left(\frac{t}{\theta}\right)^{\delta}}},\]

where $\theta>0$ is a scale parameter, and $\kappa,\delta >0$ are shape parameters.

Cumulative Distribution Function

The CDF of the GG distribution is

\[F(t;\theta,\kappa,\delta) = {\frac {\gamma \left( \frac{\kappa}{\delta},\left(\frac{t}{\theta}\right)^{\delta}\right)}{\Gamma\left(\frac{\kappa}{\delta}\right)}},\]

where where $\gamma (\cdot )$ denotes the lower incomplete gamma function. The survival function can be obtained using the relationship $S(t;\theta,\kappa,\delta)=1-F(t;\theta,\kappa,\delta)$. An interesting relationship between the Gamma CDF ($G(t;\theta,\kappa)$, scale $\theta$ and shape $\kappa$) and the GG CDF is

\[F(t;\theta,\kappa,\delta) = G\left(t^\delta; \theta^\delta, \frac{\kappa}{\delta}\right).\]

This allows the implementation of the GG CDF using the Julia command Gamma.

Hazard Function

The hazard function of the GG distribution is

\[h(t;\theta,\kappa,\delta) = \dfrac{f(t;\theta,\kappa,\delta)}{1-F(t;\theta,\kappa,\delta)}.\]

The survival function can be obtained as $S(t;\theta,\kappa,\delta)=1-F(t;\theta,\kappa,\delta)$, and the cumulative hazard function as $H(t;\theta,\kappa,\delta) = -\log S(t;\theta,\kappa,\delta)$, as usual. The connection of the GG CDF with the Gamma distribution allows for writing these functions in terms of the Julia command Gamma as shown in the following code.

The following Julia code shows the implementation of the pdf, survival function, hazard function, cumulative hazard function, quantile function, and random number generation associated to the Generalised Gamma distribution using the Julia package HazReg.jl. Some illustrative examples are also presented.

See also:

• HazReg.jl Julia Package

• The Generalised Gamma Distribution

Required packages

using Distributions
using Random
using Plots
using StatsBase
using SpecialFunctions

Functions

#=
Generalised Gamma (GG) Distribution: Hazard, cumulative hazard,
probability density function, random number generation, and survival function
=#


#= Generalised Gamma (GG) probability density function.
t: positive argument

sigma: scale parameter

nu: shape parameter

gamma: shape parameter

logd: log scale (true or false)
=#

function pdfGenGamma(t, sigma, nu, gamma, logd::Bool = false)
    val = log(gamma) .- nu * log(sigma) .- loggamma(nu/gamma) .+
    (nu - 1) * log.(t) .- (t/sigma).^gamma
    if logd
        return val
    else
        return exp.(val)
    end
end

#= Generalised Gamma (GG) survival function.
t: positive argument

sigma: scale parameter

nu: shape parameter

gamma: shape parameter

logp: log scale (true or false)
=#

function ccdfGenGamma(t, sigma, nu, gamma, logp::Bool = false)
    tp = t.^gamma
    val = logccdf.(Gamma(nu/gamma, sigma^gamma), tp)
    if logp
        return val
    else
        return exp.(val)
    end
end


#= Generalised Gamma (GG) hazard function.
t: positive argument

sigma: scale parameter

nu: shape parameter

gamma: shape parameter

logh: log scale (true or false)
=#

function hGenGamma(t, sigma, nu, gamma, logh::Bool = false)
    lpdf0 = pdfGenGamma.(t, sigma, nu, gamma, true)
    ls0 = ccdfGenGamma.(t, sigma, nu, gamma, true)
    val = lpdf0 .- ls0
    if logh
        return val
    else
        return exp.(val)
    end
end

#= Generalised Gamma (GG) cumulative hazard function.
t: positive argument

sigma: scale parameter

nu: shape parameter

gamma: shape parameter
=#

function chGenGamma(t, sigma, nu, gamma)
    H0 = -ccdfGenGamma.(t, sigma, nu, gamma, true)
    return H0
end


#= Generalised Gamma (GG) random number generation.
n: number of observations

sigma: scale parameter

nu: shape parameter

gamma: shape parameter
=#

function randGenGamma(n::Int, sigma, nu, gamma)
    p = rand(n)
    out = quantile.(Gamma(nu/gamma, sigma^gamma), p).^(1/gamma)
    return out
end

randGenGamma (generic function with 1 method)

Examples

Random number generation

#= Fix the seed =#
Random.seed!(123)
#= True values of the parameters =#
sigma0 = 1
nu0 = 3
gamma0 = 2
#= Simulation =#
sim = randGenGamma(1000, sigma0, nu0, gamma0);

1000-element Vector{Float64}:
 1.7888218704290257
 1.0191568555786799
 1.4259967883673976
 1.1024641124601198
 0.7837761620844778
 0.886450973482425
 0.9996803060951007
 1.674168073779725
 0.5387795557060041
 0.5888894291847084
 ⋮
 1.443963318480712
 1.3085724726017594
 1.3542030543904628
 1.3766989084919623
 1.0211053923759061
 1.0562031767094422
 0.39613846163829636
 1.2635696631089963
 0.6423956940997977

Some plots

#= Histogram and probability density function =#
histogram(sim, normalize=:pdf, color=:gray,
          bins = range(0, 4, length=30), label = "")
plot!(t -> pdfGenGamma(t, sigma0, nu0, gamma0),
      xlabel = "x", ylabel = "Density", title = "GenGamma distribution",
    xlims = (0,4),   xticks = 0:1:4, label = "",
    xtickfont = font(16, "Courier"),  ytickfont = font(16, "Courier"),
    xguidefontsize=18, yguidefontsize=18, linewidth=3,
    linecolor = "blue")

#= Empirical CDF and CDF =#

#= Empirical CDF=#
ecdfsim = ecdf(sim)

#= ad hoc CDF =#
function cdfGenGamma(t, sigma, nu, gamma)
val = 1 .- ccdfGenGamma.(t, sigma, nu, gamma)
return val
end

plot(x -> ecdfsim(x), 0, 5, label = "ECDF", linecolor = "gray", linewidth=3)
plot!(t -> cdfGenGamma(t, sigma0, nu0, gamma0),
      xlabel = "x", ylabel = "CDF vs. ECDF", title = "GenGamma distribution",
      xlims = (0,5),   xticks = 0:1:5, label = "CDF",
      xtickfont = font(16, "Courier"),  ytickfont = font(16, "Courier"),
      xguidefontsize=18, yguidefontsize=18, linewidth=3,
      linecolor = "blue")

#= Hazard function =#
plot(t -> hGenGamma(t, 0.5, 1.5, 0.75),
     xlabel = "x", ylabel = "Hazard", title = "GenGamma distribution",
     xlims = (0,15),   xticks = 0:1:15, label = "",
     xtickfont = font(16, "Courier"),  ytickfont = font(16, "Courier"),
     xguidefontsize=18, yguidefontsize=18, linewidth=3,
     linecolor = "blue")