#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
module stdlib_stats_distribution_normal
use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
use stdlib_kinds, only: sp, dp, xdp, qp, int32
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
impure 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
impure function rvs_norm_0_${t1[0]}$${k1}$ () result(res)
!
! Standard normal random variate (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
impure elemental &
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}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else
res = rvs_norm_0_${t1[0]}$${k1}$ ()
res = res*scale + loc
end if
end function rvs_norm_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
impure elemental 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
impure 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 (.not. zig_norm_initialized) call zigset
if (scale <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if
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
impure 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
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}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else
res = exp(-0.5_${k1}$*((x - loc)/scale)*(x - loc)/scale)/ &
(sqrt_2_Pi*scale)
end if
end function pdf_norm_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
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
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}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else
res = erfc(-(x - loc)/(scale*sqrt_2))/2.0_${k1}$
end if
end function cdf_norm_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
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
