Code
using Distributions
using Random
using Plots
using StatsBase
using SpecialFunctions
The Generalised Gamma (GG) distribution (Stacy 1962) 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 et al. (2021) and Jones and Noufaily (2015).
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.
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
.
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:
using Distributions
using Random
using Plots
using StatsBase
using SpecialFunctions
#=
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)
= log(gamma) .- nu * log(sigma) .- loggamma(nu/gamma) .+
val - 1) * log.(t) .- (t/sigma).^gamma
(nu 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)
= t.^gamma
tp = logccdf.(Gamma(nu/gamma, sigma^gamma), tp)
val 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)
= pdfGenGamma.(t, sigma, nu, gamma, true)
lpdf0 = ccdfGenGamma.(t, sigma, nu, gamma, true)
ls0 = lpdf0 .- ls0
val 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)
= -ccdfGenGamma.(t, sigma, nu, gamma, true)
H0 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)
= rand(n)
p = quantile.(Gamma(nu/gamma, sigma^gamma), p).^(1/gamma)
out return out
end
randGenGamma (generic function with 1 method)
#= Fix the seed =#
Random.seed!(123)
#= True values of the parameters =#
= 1
sigma0 = 3
nu0 = 2
gamma0 #= Simulation =#
= randGenGamma(1000, sigma0, nu0, gamma0); sim
#= Histogram and probability density function =#
histogram(sim, normalize=:pdf, color=:gray,
= range(0, 4, length=30), label = "")
bins plot!(t -> pdfGenGamma(t, sigma0, nu0, gamma0),
= "x", ylabel = "Density", title = "GenGamma distribution",
xlabel = (0,4), xticks = 0:1:4, label = "",
xlims = font(16, "Courier"), ytickfont = font(16, "Courier"),
xtickfont =18, yguidefontsize=18, linewidth=3,
xguidefontsize= "blue") linecolor
#= Empirical CDF and CDF =#
#= Empirical CDF=#
= ecdf(sim)
ecdfsim
#= ad hoc CDF =#
function cdfGenGamma(t, sigma, nu, gamma)
= 1 .- ccdfGenGamma.(t, sigma, nu, gamma)
val return val
end
plot(x -> ecdfsim(x), 0, 5, label = "ECDF", linecolor = "gray", linewidth=3)
plot!(t -> cdfGenGamma(t, sigma0, nu0, gamma0),
= "x", ylabel = "CDF vs. ECDF", title = "GenGamma distribution",
xlabel = (0,5), xticks = 0:1:5, label = "CDF",
xlims = font(16, "Courier"), ytickfont = font(16, "Courier"),
xtickfont =18, yguidefontsize=18, linewidth=3,
xguidefontsize= "blue") linecolor
#= Hazard function =#
plot(t -> hGenGamma(t, 0.5, 1.5, 0.75),
= "x", ylabel = "Hazard", title = "GenGamma distribution",
xlabel = (0,15), xticks = 0:1:15, label = "",
xlims = font(16, "Courier"), ytickfont = font(16, "Courier"),
xtickfont =18, yguidefontsize=18, linewidth=3,
xguidefontsize= "blue") linecolor