#:include "common.fypp"
#:set GAMMA_REAL_KINDS_TYPES = [item for item in REAL_KINDS_TYPES if item[0] in ('sp', 'dp', 'xdp')]
#:set GAMMA_CMPLX_KINDS_TYPES = [item for item in CMPLX_KINDS_TYPES if item[0] in ('sp', 'dp', 'xdp')]
#:set RC_KINDS_TYPES = GAMMA_REAL_KINDS_TYPES + GAMMA_CMPLX_KINDS_TYPES
Module stdlib_stats_distribution_gamma
use ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_nan
use stdlib_kinds, only : sp, dp, xdp
use stdlib_error, only : error_stop
use stdlib_optval, only : optval
use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform
use stdlib_stats_distribution_normal, only : rnor=>rvs_normal
use stdlib_specialfunctions_gamma, only : regamma_p=>regularized_gamma_p
implicit none
intrinsic :: log_gamma
private
public :: rvs_gamma
public :: pdf_gamma
public :: cdf_gamma
interface rvs_gamma
!! Version experimental
!!
!! Gamma Distribution Random Variates
!! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html#
!! rvs_gamma-gamma-distribution-random-variates))
!!
#:for k1, t1 in RC_KINDS_TYPES
module procedure gamma_dist_rvs_1_${t1[0]}$${k1}$ ! 1 argument
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
module procedure gamma_dist_rvs_${t1[0]}$${k1}$ ! 2 arguments
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
module procedure gamma_dist_rvs_array_${t1[0]}$${k1}$ ! 3 arguments
#:endfor
end interface rvs_gamma
interface pdf_gamma
!! Version experimental
!!
!! Gamma Distribution Probability Density Function
!! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html#
!! pdf_gamma-gamma-distribution-probability-density-function))
!!
#:for k1, t1 in RC_KINDS_TYPES
module procedure gamma_dist_pdf_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
module procedure gamma_dist_pdf_impure_${t1[0]}$${k1}$
#:endfor
end interface pdf_gamma
interface cdf_gamma
!! Version experimental
!!
!! Gamma Distribution Cumulative Distribution Function
!! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html#
!! cdf_gamma-gamma-distribution-cumulative-distribution-function))
!!
#:for k1, t1 in RC_KINDS_TYPES
module procedure gamma_dist_cdf_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
module procedure gamma_dist_cdf_impure_${t1[0]}$${k1}$
#:endfor
end interface cdf_gamma
contains
#:for k1, t1 in GAMMA_REAL_KINDS_TYPES
impure elemental function gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) result(res)
! Gamma distribution random variate. "A Simple Method for Generating Gamma
! Variables", G. Marsaglia & W. W. Tsang, ACM Transactions on Mathematical
! Software, 26(3), 2000, p. 363
!
${t1}$, intent(in) :: shape
${t1}$ :: res
${t1}$ :: x, v, u, zz, d, c
${t1}$, parameter :: sq = 0.0331_${k1}$
if(shape <= 0.0_${k1}$) then
res = ieee_value(1.0_${k1}$, ieee_quiet_nan)
return
end if
zz = shape
if(zz < 1._${k1}$) zz = 1._${k1}$ + zz
!shift shape parameter > 1
d = zz - 1._${k1}$ / 3._${k1}$
c = 1._${k1}$ / (3._${k1}$ * sqrt(d))
do
do
x = rnor(0.0_${k1}$, 1.0_${k1}$)
v = 1._${k1}$ + c * x
v = v * v * v
if(v > 0._${k1}$) exit
end do
x = x * x
u = uni(1.0_${k1}$)
if(u < (1._${k1}$ - sq * x * x)) exit
if(log(u) < 0.5_${k1}$ * x + d * (1._${k1}$ - v + log(v))) exit
end do
res = d * v
if(shape < 1._${k1}$) then
!restore shape parameter < 1
u = uni(1.0_${k1}$)
res = res * u ** (1._${k1}$ / shape)
endif
end function gamma_dist_rvs_1_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_CMPLX_KINDS_TYPES
impure elemental function gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) result(res)
! Complex parameter gamma distributed. The real part and imaginary part are
! independent of each other.
!
${t1}$, intent(in) :: shape
${t1}$ :: res
res = cmplx(gamma_dist_rvs_1_r${k1}$(shape%re), &
gamma_dist_rvs_1_r${k1}$(shape%im), kind=${k1}$)
end function gamma_dist_rvs_1_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_REAL_KINDS_TYPES
impure elemental function gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate, loc) &
result(res)
!
${t1}$, intent(in) :: shape, rate
${t1}$, intent(in), optional :: loc
${t1}$ :: res
${t1}$ :: loc_
loc_ = optval(loc, 0.0_${k1}$)
if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs): Gamma" &
//" distribution rate parameter must be greater than zero")
res = gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) / rate + loc_
end function gamma_dist_rvs_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_CMPLX_KINDS_TYPES
impure elemental function gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate, loc) &
result(res)
! Complex parameter gamma distributed. The real part and imaginary part are &
! independent of each other.
!
${t1}$, intent(in) :: shape, rate
${t1}$, intent(in), optional :: loc
${t1}$ :: res
${t1}$ :: loc_
loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$))
res = cmplx(gamma_dist_rvs_r${k1}$(shape%re, rate%re) + loc_%re, &
gamma_dist_rvs_r${k1}$(shape%im, rate%im) + loc_%im, kind=${k1}$)
end function gamma_dist_rvs_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_REAL_KINDS_TYPES
function gamma_dist_rvs_array_${t1[0]}$${k1}$(shape, rate, array_size, loc) &
result(res)
!
${t1}$, intent(in) :: shape, rate
integer, intent(in) :: array_size
${t1}$, intent(in), optional :: loc
${t1}$ :: res(array_size)
integer :: i
${t1}$ :: loc_
loc_ = optval(loc, 0.0_${k1}$)
do i = 1, array_size
res(i) = gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate, loc_)
end do
end function gamma_dist_rvs_array_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_CMPLX_KINDS_TYPES
function gamma_dist_rvs_array_${t1[0]}$${k1}$(shape, rate, array_size, loc) &
result(res)
! Complex parameter gamma distributed. The real part and imaginary part are &
! independent of each other.
!
${t1}$, intent(in) :: shape, rate
integer, intent(in) :: array_size
${t1}$, intent(in), optional :: loc
${t1}$ :: res(array_size)
integer :: i
${t1}$ :: loc_
loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$))
do i = 1, array_size
res(i) = gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate, loc_)
end do
end function gamma_dist_rvs_array_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_REAL_KINDS_TYPES
elemental function gamma_dist_pdf_${t1[0]}$${k1}$(x, shape, rate, loc) &
result(res)
! Gamma distribution probability density function
!
${t1}$, intent(in) :: x, shape, rate
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res
${t1}$ :: xs
${t1}$ :: loc_
loc_ = optval(loc, 0.0_${k1}$)
xs = x - loc_
if(rate <= 0.0_${k1}$ .or. shape <= 0.0_${k1}$) then
res = ieee_value(1.0_${k1}$, ieee_quiet_nan)
return
end if
if(xs == 0.0_${k1}$) then
if(shape == 1.0_${k1}$) then
res = rate
else if(shape > 1.0_${k1}$) then
res = 0.0_${k1}$
else
res = huge(res)
endif
else if(xs < 0.0_${k1}$) then
res = 0.0_${k1}$
else
res = exp((shape - 1._${k1}$) * log(xs) - xs * rate + shape * &
log(rate) - log_gamma(shape))
endif
end function gamma_dist_pdf_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_REAL_KINDS_TYPES
impure elemental function gamma_dist_pdf_impure_${t1[0]}$${k1}$(x, shape, rate, err, loc) &
result(res)
! Gamma distribution probability density function (impure wrapper)
!
${t1}$, intent(in) :: x, shape, rate
integer, intent(out) :: err
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res
res = gamma_dist_pdf_${t1[0]}$${k1}$(x, shape, rate, loc)
err = 0
if(ieee_is_nan(res)) err = 1
end function gamma_dist_pdf_impure_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_CMPLX_KINDS_TYPES
elemental function gamma_dist_pdf_${t1[0]}$${k1}$(x, shape, rate, loc) &
result(res)
! Complex parameter gamma distributed. The real part and imaginary part are &
! independent of each other.
!
${t1}$, intent(in) :: x, shape, rate
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res
${t1}$ :: loc_
loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$))
res = gamma_dist_pdf_r${k1}$(x%re, shape%re, rate%re, loc_%re)
res = res * gamma_dist_pdf_r${k1}$(x%im, shape%im, rate%im, loc_%im)
end function gamma_dist_pdf_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_CMPLX_KINDS_TYPES
impure elemental function gamma_dist_pdf_impure_${t1[0]}$${k1}$(x, shape, rate, err, loc) &
result(res)
! Complex parameter gamma distributed. The real part and imaginary part are &
! independent of each other.
!
${t1}$, intent(in) :: x, shape, rate
integer, intent(out) :: err
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res
${t1}$ :: loc_
loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$))
res = gamma_dist_pdf_r${k1}$(x%re, shape%re, rate%re, loc_%re)
res = res * gamma_dist_pdf_r${k1}$(x%im, shape%im, rate%im, loc_%im)
err = 0
if(ieee_is_nan(res)) err = 1
end function gamma_dist_pdf_impure_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_REAL_KINDS_TYPES
impure elemental function gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate, loc) &
result(res)
! Gamma distribution cumulative distribution function
!
${t1}$, intent(in) :: x, shape, rate
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res, xs
${t1}$ :: loc_
loc_ = optval(loc, 0.0_${k1}$)
if(rate <= 0.0_${k1}$ .or. shape <= 0.0_${k1}$) then
res = ieee_value(1.0_${k1}$, ieee_quiet_nan)
return
end if
xs = x - loc_
if(xs <= 0.0_${k1}$) then
res = 0.0_${k1}$
return
end if
res = real(regamma_p(shape, xs * rate), kind=${k1}$)
end function gamma_dist_cdf_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_REAL_KINDS_TYPES
impure elemental function gamma_dist_cdf_impure_${t1[0]}$${k1}$(x, shape, rate, err, loc) &
result(res)
! Gamma distribution cumulative distribution function (impure wrapper)
!
${t1}$, intent(in) :: x, shape, rate
integer, intent(out) :: err
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res
res = gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate, loc)
err = 0
if(ieee_is_nan(res)) err = 1
end function gamma_dist_cdf_impure_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_CMPLX_KINDS_TYPES
impure elemental function gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate, loc) &
result(res)
! Complex parameter gamma distributed. The real part and imaginary part are &
! independent of each other.
!
${t1}$, intent(in) :: x, shape, rate
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res
${t1}$ :: loc_
loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$))
res = gamma_dist_cdf_r${k1}$(x%re, shape%re, rate%re, loc_%re)
res = res * gamma_dist_cdf_r${k1}$(x%im, shape%im, rate%im, loc_%im)
end function gamma_dist_cdf_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in GAMMA_CMPLX_KINDS_TYPES
impure elemental function gamma_dist_cdf_impure_${t1[0]}$${k1}$(x, shape, rate, err, loc) &
result(res)
! Complex parameter gamma distributed. The real part and imaginary part are &
! independent of each other.
!
${t1}$, intent(in) :: x, shape, rate
integer, intent(out) :: err
${t1}$, intent(in), optional :: loc
real(${k1}$) :: res
${t1}$ :: loc_
loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$))
res = gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate, loc_)
err = 0
if(ieee_is_nan(res)) err = 1
end function gamma_dist_cdf_impure_${t1[0]}$${k1}$
#:endfor
end module stdlib_stats_distribution_gamma
