stdlib_stats_distribution_exponential.fypp Source File


This file depends on

sourcefile~~stdlib_stats_distribution_exponential.fypp~~EfferentGraph sourcefile~stdlib_stats_distribution_exponential.fypp stdlib_stats_distribution_exponential.fypp sourcefile~stdlib_kinds.fypp stdlib_kinds.fypp sourcefile~stdlib_stats_distribution_exponential.fypp->sourcefile~stdlib_kinds.fypp sourcefile~stdlib_random.fypp stdlib_random.fypp sourcefile~stdlib_stats_distribution_exponential.fypp->sourcefile~stdlib_random.fypp sourcefile~stdlib_error.f90 stdlib_error.f90 sourcefile~stdlib_stats_distribution_exponential.fypp->sourcefile~stdlib_error.f90 sourcefile~stdlib_stats_distribution_uniform.fypp stdlib_stats_distribution_uniform.fypp sourcefile~stdlib_stats_distribution_exponential.fypp->sourcefile~stdlib_stats_distribution_uniform.fypp sourcefile~stdlib_random.fypp->sourcefile~stdlib_kinds.fypp sourcefile~stdlib_random.fypp->sourcefile~stdlib_error.f90 sourcefile~stdlib_optval.fypp stdlib_optval.fypp sourcefile~stdlib_random.fypp->sourcefile~stdlib_optval.fypp sourcefile~stdlib_error.f90->sourcefile~stdlib_optval.fypp sourcefile~stdlib_stats_distribution_uniform.fypp->sourcefile~stdlib_kinds.fypp sourcefile~stdlib_stats_distribution_uniform.fypp->sourcefile~stdlib_random.fypp sourcefile~stdlib_stats_distribution_uniform.fypp->sourcefile~stdlib_error.f90 sourcefile~stdlib_optval.fypp->sourcefile~stdlib_kinds.fypp

Contents


Source Code

#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
module stdlib_stats_distribution_exponential
    use stdlib_kinds, only : sp, dp, xdp, qp, int32
    use stdlib_error, only : error_stop
    use stdlib_random, only : dist_rand
    use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform

    implicit none
    private

    integer :: ke(0:255)
    real(dp) :: we(0:255), fe(0:255)
    logical  :: zig_exp_initialized = .false.

    public :: rvs_exp
    public :: pdf_exp
    public :: cdf_exp



    interface rvs_exp
    !! Version experimental
    !!
    !! Exponential Distribution Random Variates
    !! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
    !! rvs_exp-exponential-distribution-random-variates))
    !!
        module procedure rvs_exp_0_rsp                 !0 dummy variable

        #:for k1, t1 in RC_KINDS_TYPES
        module procedure rvs_exp_${t1[0]}$${k1}$       !1 dummy variable
        #:endfor

        #:for k1, t1 in RC_KINDS_TYPES
        module procedure rvs_exp_array_${t1[0]}$${k1}$ !2 dummy variables
        #:endfor
    end interface rvs_exp



    interface pdf_exp
    !! Version experimental
    !!
    !! Exponential Distribution Probability Density Function
    !! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
    !! pdf_exp-exponential-distribution-probability-density-function))
    !!
        #:for k1, t1 in RC_KINDS_TYPES
        module procedure pdf_exp_${t1[0]}$${k1}$
        #:endfor
    end interface pdf_exp



    interface cdf_exp
    !! Version experimental
    !!
    !! Exponential Cumulative Distribution Function
    !! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
    !! cdf_exp-exponential-distribution-cumulative-distribution-function))
    !!
        #:for k1, t1 in RC_KINDS_TYPES
        module procedure cdf_exp_${t1[0]}$${k1}$
        #:endfor
    end interface cdf_exp





contains

    subroutine zigset
    ! Marsaglia & Tsang generator for random normals & random exponentials.
    ! Translated from C by Alan Miller (amiller@bigpond.net.au)
    !
    ! Marsaglia, G. & Tsang, W.W. (2000) 'The ziggurat method for generating
    ! random variables', J. Statist. Software, v5(8).
    !
    ! This is an electronic journal which can be downloaded from:
    ! http://www.jstatsoft.org/v05/i08
    !
    ! Latest version - 1 January 2001
    !
        real(dp), parameter :: M2 = 2147483648.0_dp, ve = 0.003949659822581572_dp
        real(dp), parameter :: ONE = 1.0_dp
        real(dp)            :: de, te, q
        integer :: i

        de = 7.697117470131487_dp
        te = de
    !tables for random exponetials
        q = ve * exp(de)
        ke(0) = int((de / q) * M2, kind = int32)
        ke(1) = 0
        we(0) = q / M2
        we(255) = de / M2
        fe(0) = ONE
        fe(255) = exp(- de)
        do  i = 254, 1, -1
            de = -log(ve / de + exp(- de))
            ke(i+1) = int(M2 * (de / te), kind = int32)
            te = de
            fe(i) = exp(- de)
            we(i) = de / M2
        end do
        zig_exp_initialized = .true.
    end subroutine zigset




    #:for k1, t1 in REAL_KINDS_TYPES
    function rvs_exp_0_${t1[0]}$${k1}$( ) result(res)
    !
    ! Standard exponential random variate (lambda=1)
    !
        ${t1}$ :: res, x
        ${t1}$, parameter :: r = 7.69711747013104972_${k1}$
        integer :: jz, iz

        if(.not. zig_exp_initialized ) call zigset
        iz = 0
        jz = dist_rand(1_int32)                !32bit random integer
        iz = iand( jz, 255 )                   !random integer in [0, 255]
        if( abs( jz ) < ke(iz) ) then
            res = abs(jz) * we(iz)
        else
            L1: do
                if( iz == 0 ) then
                    res = r - log( uni(1.0_${k1}$) )
                    exit L1
                end if
                x = abs( jz ) * we(iz)
                if(fe(iz) + uni(1.0_${k1}$) * (fe(iz-1) - fe(iz)) < exp(-x)) then
                    res = x
                    exit L1
                end if
                jz = dist_rand(1_int32)
                iz = iand( jz, 255 )
                if( abs( jz ) < ke(iz) ) then
                    res = abs( jz ) * we(iz)
                    exit L1
                end if
           end do L1
       endif
    end function rvs_exp_0_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in REAL_KINDS_TYPES
    function rvs_exp_${t1[0]}$${k1}$(lambda) result(res)
    !
    ! Exponential distributed random variate
    !
        ${t1}$, intent(in) :: lambda
        ${t1}$ :: res


        if(lambda <= 0.0_${k1}$) call error_stop("Error(rvs_exp): Exponen"   &
            //"tial distribution lambda parameter must be greater than zero")
        res = rvs_exp_0_${t1[0]}$${k1}$(  )
        res = res / lambda
    end function rvs_exp_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in CMPLX_KINDS_TYPES
    function rvs_exp_${t1[0]}$${k1}$(lambda) result(res)
        ${t1}$, intent(in) :: lambda
        ${t1}$ :: res
        real(${k1}$) :: tr, ti

        tr = rvs_exp_r${k1}$(lambda % re)
        ti = rvs_exp_r${k1}$(lambda % im)
        res = cmplx(tr, ti, kind=${k1}$)
    end function rvs_exp_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in REAL_KINDS_TYPES
    function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res)
        ${t1}$, intent(in) :: lambda
        integer, intent(in) :: array_size
        ${t1}$ :: res(array_size), x, re
        ${t1}$, parameter :: r = 7.69711747013104972_${k1}$
        integer :: jz, iz, i

        if(lambda <= 0.0_${k1}$) call error_stop("Error(rvs_exp_array): Exp" &
            //"onential distribution lambda parameter must be greater than zero")

        if(.not. zig_exp_initialized) call zigset
        do i = 1, array_size
            iz = 0
            jz = dist_rand(1_int32)
            iz = iand( jz, 255 )
            if( abs( jz ) < ke(iz) ) then
                re = abs(jz) * we(iz)
            else
                L1: do
                    if( iz == 0 ) then
                        re = r - log( uni(1.0_${k1}$) )
                        exit L1
                    end if
                    x = abs( jz ) * we(iz)
                    if(fe(iz) + uni(1.0_${k1}$)*(fe(iz-1)-fe(iz)) < exp(-x)) then
                        re = x
                        exit L1
                    end if
                    jz = dist_rand(1_int32)
                    iz = iand( jz, 255 )
                    if( abs( jz ) < ke(iz) ) then
                        re = abs( jz ) * we(iz)
                        exit L1
                    end if
               end do L1
            endif
            res(i) = re / lambda
        end do
    end function rvs_exp_array_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in CMPLX_KINDS_TYPES
    function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res)
        ${t1}$, intent(in) :: lambda
        integer, intent(in) :: array_size
        ${t1}$ :: res(array_size)
        integer :: i
        real(${k1}$) :: tr, ti

        do i = 1, array_size
            tr = rvs_exp_r${k1}$(lambda % re)
            ti = rvs_exp_r${k1}$(lambda % im)
            res(i) = cmplx(tr, ti, kind=${k1}$)
        end do
    end function rvs_exp_array_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
    !
    ! Exponential Distribution Probability Density Function
    !
        ${t1}$, intent(in) :: x, lambda
        real(${k1}$) :: res

        if(lambda <= 0.0_${k1}$) call error_stop("Error(pdf_exp): Expon"     &
            //"ential distribution lambda parameter must be greater than zero")
        if(x < 0.0_${k1}$) call error_stop("Error(pdf_exp): Exponential"     &
            //" distribution variate x must be non-negative")
        res = exp(- x * lambda) * lambda
    end function pdf_exp_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in CMPLX_KINDS_TYPES
    impure elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
        ${t1}$, intent(in) :: x, lambda
        real(${k1}$) :: res

        res = pdf_exp_r${k1}$(x % re, lambda % re)
        res = res * pdf_exp_r${k1}$(x % im, lambda % im)
    end function pdf_exp_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
    !
    ! Exponential Distribution Cumulative Distribution Function
    !
        ${t1}$, intent(in) :: x, lambda
        real(${k1}$) :: res

        if(lambda <= 0.0_${k1}$) call error_stop("Error(cdf_exp): Expon"     &
            //"ential distribution lambda parameter must be greater than zero")
        if(x < 0.0_${k1}$) call error_stop("Error(cdf_exp): Exponential"     &
            //" distribution variate x must be non-negative")
        res = 1.0_${k1}$ - exp(- x * lambda)
    end function cdf_exp_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in CMPLX_KINDS_TYPES
    impure elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
        ${t1}$, intent(in) :: x, lambda
        real(${k1}$) :: res

        res = cdf_exp_r${k1}$(x % re, lambda % re)
        res = res * cdf_exp_r${k1}$(x % im, lambda % im)
    end function cdf_exp_${t1[0]}$${k1}$

    #:endfor

end module stdlib_stats_distribution_exponential