stdlib_stats_distribution_normal.fypp Source File


This file depends on

sourcefile~~stdlib_stats_distribution_normal.fypp~~EfferentGraph sourcefile~stdlib_stats_distribution_normal.fypp stdlib_stats_distribution_normal.fypp sourcefile~stdlib_kinds.fypp stdlib_kinds.fypp sourcefile~stdlib_stats_distribution_normal.fypp->sourcefile~stdlib_kinds.fypp sourcefile~stdlib_random.fypp stdlib_random.fypp sourcefile~stdlib_stats_distribution_normal.fypp->sourcefile~stdlib_random.fypp sourcefile~stdlib_error.f90 stdlib_error.f90 sourcefile~stdlib_stats_distribution_normal.fypp->sourcefile~stdlib_error.f90 sourcefile~stdlib_stats_distribution_uniform.fypp stdlib_stats_distribution_uniform.fypp sourcefile~stdlib_stats_distribution_normal.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_normal
    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

    real(dp), parameter  :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp
    integer :: kn(0:127)
    real(dp) :: wn(0:127), fn(0:127)
    logical  :: zig_norm_initialized = .false.

    public :: rvs_normal
    public :: pdf_normal
    public :: cdf_normal



    interface rvs_normal
    !! version: experimental
    !!
    !! Normal Distribution Random Variates
    !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html#
    !! rvs_normal-normal-distribution-random-variates))
    !!
        module procedure rvs_norm_0_rsp                  !0 dummy variable

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

        #:for k1, t1 in RC_KINDS_TYPES
        module procedure rvs_norm_array_${t1[0]}$${k1}$  !3 dummy variables
        #:endfor
    end interface rvs_normal



    interface pdf_normal
    !! version: experimental
    !!
    !! Normal Distribution Probability Density Function
    !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html#
    !! pdf_normal-normal-distribution-probability-density-function))
    !!
        #:for k1, t1 in RC_KINDS_TYPES
        module procedure pdf_norm_${t1[0]}$${k1}$
        #:endfor
    end interface pdf_normal



    interface cdf_normal
    !! version: experimental
    !!
    !! Normal Distribution Cumulative Distribution Function
    !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html#
    !! cdf_normal-normal-distribution-cumulative-distribution-function))
    !!
        #:for k1, t1 in RC_KINDS_TYPES
        module procedure cdf_norm_${t1[0]}$${k1}$
        #:endfor
    end interface cdf_normal





contains

    subroutine zigset
    ! Marsaglia & Tsang generator for random normals & random exponentials.
    ! Translated from C by Alan Miller (amiller@bigpond.net.au), released as public
    ! domain (https://jblevins.org/mirror/amiller/)
    !
    ! 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 :: M1 = 2147483648.0_dp, vn = 0.00991256303526217_dp
        real(dp)            :: dn, tn, q
        integer :: i

        dn = 3.442619855899_dp
        tn = dn
    !tables for random normals
        q = vn * exp(HALF * dn * dn)
        kn(0) = int((dn / q) * M1, kind = int32)
        kn(1) = 0
        wn(0) = q / M1
        wn(127) = dn / M1
        fn(0) = ONE
        fn(127) = exp( -HALF * dn * dn )
        do  i = 126, 1, -1
            dn = sqrt( -TWO * log( vn / dn + exp( -HALF * dn * dn ) ) )
            kn(i+1) = int((dn / tn) * M1, kind = int32)
            tn = dn
            fn(i) = exp(-HALF * dn * dn)
            wn(i) = dn / M1
        end do
        zig_norm_initialized = .true.
    end subroutine zigset



    #:for k1, t1 in REAL_KINDS_TYPES
    function rvs_norm_0_${t1[0]}$${k1}$( ) result(res)
    !
    ! Standard normal random vairate (0,1)
    !
        ${t1}$ :: res
        ${t1}$, parameter  ::  r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r
        ${t1}$ ::  x, y
        integer :: hz, iz

        if(.not. zig_norm_initialized) call zigset
        iz = 0
        hz = dist_rand(1_int32)          !32bit random integer
        iz = iand( hz, 127 )             !random integer in [0, 127]
        if( abs( hz ) < kn(iz) ) then
            res = hz * wn(iz)
        else
            L1: do
                L2: if( iz == 0 ) then
                    do
                        x = -log( uni(1.0_${k1}$) ) * rr
                        y = -log( uni(1.0_${k1}$) )
                        if( y + y >= x * x ) exit
                    end do
                    res = r + x
                    if( hz <= 0 ) res = -res
                    exit L1
                end if L2
                x = hz * wn(iz)
                if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) <            &
                    exp(-HALF * x * x) ) then
                    res = x
                    exit L1
                end if
                hz = dist_rand(1_int32)
                iz = iand( hz, 127 )
                if( abs( hz ) < kn(iz) ) then
                    res = hz * wn(iz)
                    exit L1
                end if
            end do L1
        end if
    end function rvs_norm_0_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res)
    !
    ! Normal random variate (loc, scale)
    !
        ${t1}$, intent(in) :: loc, scale
        ${t1}$ :: res

        if(scale == 0._${k1}$) call error_stop("Error(rvs_norm): Normal"       &
            //" distribution scale parameter must be non-zero")
        res = rvs_norm_0_${t1[0]}$${k1}$(  )
        res = res * scale + loc
    end function rvs_norm_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in CMPLX_KINDS_TYPES
    function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res)
    !
    ! Normally distributed complex. The real part and imaginary part are       &
    ! independent of each other.
    !
        ${t1}$, intent(in) :: loc, scale
        ${t1}$ :: res
        real(${k1}$) :: tr, ti

        tr = rvs_norm_r${k1}$(loc % re, scale % re)
        ti = rvs_norm_r${k1}$(loc % im, scale % im)
        res = cmplx(tr, ti, kind=${k1}$)
    end function rvs_norm_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
        ${t1}$, intent(in) :: loc, scale
        integer, intent(in) :: array_size
        ${t1}$ :: res(array_size)
        ${t1}$, parameter  ::  r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r
        ${t1}$ ::  x, y, re
        integer :: hz, iz, i

        if(scale == 0._${k1}$) call error_stop("Error(rvs_norm_array): Normal" &
            //"distribution scale parameter must be non-zero")
        if(.not. zig_norm_initialized) call zigset
        do i = 1, array_size
            iz = 0
            hz = dist_rand(1_int32)
            iz = iand( hz, 127 )
            if( abs( hz ) < kn(iz) ) then
                re = hz * wn(iz)
            else
                L1: do
                    L2: if( iz == 0 ) then
                        do
                            x = -log( uni(1.0_${k1}$) ) * rr
                            y = -log( uni(1.0_${k1}$) )
                            if( y + y >= x * x ) exit
                        end do
                        re = r + x
                        if( hz <= 0 ) re = -re
                        exit L1
                    end if L2
                    x = hz * wn(iz)
                    if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) <        &
                        exp(-HALF * x * x) ) then
                        re = x
                        exit L1
                    end if

                    hz = dist_rand(1_int32)
                    iz = iand( hz, 127 )
                    if( abs( hz ) < kn(iz) ) then
                        re = hz * wn(iz)
                        exit L1
                    end if
                end do L1
            end if
            res(i) = re * scale + loc
        end do
    end function rvs_norm_array_${t1[0]}$${k1}$

    #:endfor



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

        do i = 1, array_size
            tr = rvs_norm_r${k1}$(loc % re, scale % re)
            ti = rvs_norm_r${k1}$(loc % im, scale % im)
            res(i) = cmplx(tr, ti, kind=${k1}$)
        end do
    end function rvs_norm_array_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
    !
    ! Normal distribution probability density function
    !
        ${t1}$, intent(in) :: x, loc, scale
        ${t1}$ :: res
        ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$))

        if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal"       &
            //"distribution scale parameter must be non-zero")
        res = exp(- 0.5_${k1}$ * ((x - loc) / scale) * (x - loc) / scale) /    &
              (sqrt_2_Pi * scale)
    end function pdf_norm_${t1[0]}$${k1}$

    #:endfor



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

        res = pdf_norm_r${k1}$(x % re, loc % re, scale % re)
        res = res * pdf_norm_r${k1}$(x % im, loc % im, scale % im)
    end function pdf_norm_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
    !
    ! Normal distribution cumulative distribution function
    !
        ${t1}$, intent(in) :: x, loc, scale
        ${t1}$ :: res
        ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$)

        if(scale == 0._${k1}$) call error_stop("Error(cdf_norm): Normal"       &
            //"distribution scale parameter must be non-zero")
        res = erfc(- (x - loc) / (scale * sqrt_2)) / 2.0_${k1}$
    end function cdf_norm_${t1[0]}$${k1}$

    #:endfor



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

        res = cdf_norm_r${k1}$(x % re, loc % re, scale % re)
        res = res * cdf_norm_r${k1}$(x % im, loc % im, scale % im)
    end function cdf_norm_${t1[0]}$${k1}$

    #:endfor

end module stdlib_stats_distribution_normal