diff --git a/doc/specs/index.md b/doc/specs/index.md index 153b769f7..d96c3acf4 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -29,6 +29,7 @@ This is and index/directory of the specifications (specs) for each new module/fe - [stats](./stdlib_stats.html) - Descriptive Statistics - [stats_distributions_uniform](./stdlib_stats_distribution_uniform.html) - Uniform Probability Distribution - [stats_distributions_normal](./stdlib_stats_distribution_normal.html) - Normal Probability Distribution + - [stats_distributions_exponential](./stdlib_stats_distribution_exponential.html) - Exponential Probability Distribution - [string\_type](./stdlib_string_type.html) - Basic string support - [strings](./stdlib_strings.html) - String handling and manipulation routines - [stringlist_type](./stdlib_stringlist_type.html) - 1-Dimensional list of strings diff --git a/doc/specs/stdlib_stats_distribution_exponential.md b/doc/specs/stdlib_stats_distribution_exponential.md new file mode 100644 index 000000000..9a53bee95 --- /dev/null +++ b/doc/specs/stdlib_stats_distribution_exponential.md @@ -0,0 +1,248 @@ +--- +title: stats_distribution_exponential +--- + +# Statistical Distributions -- Exponential Distribution Module + +[TOC] + +## `rvs_exp` - exponential distribution random variates + +### Status + +Experimental + +### Description + +An exponential distribution is the distribution of time between events in a Poisson point process. The inverse scale parameter `lambda` specifies the average time between events, also called the rate of events. + +Without argument the function returns a random sample from the standard exponential distribution `E(1)` with `lambda = 1`. + +With a single argument, the function returns a random sample from the exponential distribution `E(lambda)`. +For complex arguments, the real and imaginary parts are sampled independently of each other. + +With two arguments the function returns a rank one array of exponentially distributed random variates. + +Note: the algorithm used for generating exponetial random variates is fundamentally limited to double precision. Ref.: Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8). + +### Syntax + +`result = [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]]([lambda] [[, array_size]])` + +### Class + +Function + +### Arguments + +`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. The value of `lambda` has to be non-negative. + +`array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind. + +### Return value + +The result is a scalar or rank one array with a size of `array_size`, and has the same type of `lambda`. + +### Example + +```fortran +program demo_exponential_rvs + use stdlib_random, only : random_seed + use stdlib_stats_distribution_exponential, only: rexp => rvs_exp + + implicit none + real :: a(2,3,4) + complex :: scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, rexp( ) !single standard exponential random variate + +! 0.358690143 + + print *, rexp(2.0) !exponential random variate with lambda=2.0 + +! 0.816459715 + + print *, rexp(0.3, 10) !an array of 10 variates with lambda=0.3 + +! 1.84008647E-02 3.59742008E-02 0.136567295 0.262772143 3.62352766E-02 +! 0.547133625 0.213591918 4.10784185E-02 0.583882213 0.671128035 + + scale = (2.0, 0.7) + print *, rexp(scale) + !single complex exponential random variate with real part of lambda=2.0; + !imagainary part of lambda=0.7 + +! (1.41435969,4.081114382E-02) + +end program demo_exponential_rvs +``` + +## `pdf_exp` - exponential distribution probability density function + +### Status + +Experimental + +### Description + +The probability density function (pdf) of the single real variable exponential distribution: + +$$f(x)=\begin{cases} \lambda e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$ + +For a complex variable (x + y i) with independent real x and imaginary y parts, the joint probability density function +is the product of the corresponding marginal pdf of real and imaginary pdf (for more details, see +"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): + +$$f(x+\mathit{i}y)=f(x)f(y)=\begin{cases} \lambda_{x} \lambda_{y} e^{-(\lambda_{x} x + \lambda_{y} y)} &x\geqslant 0, y\geqslant 0 \\\\ 0 &otherwise\end{}$$ + +### Syntax + +`result = [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]](x, lambda)` + +### Class + +Elemental function + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. + +`lambda`: has `intent(in)` and is a scalar of type `real` or `complex`. + +All arguments must have the same type. + +### Return value + +The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments. + +### Example + +```fortran +program demo_exponential_pdf + use stdlib_random, only : random_seed + use stdlib_stats_distribution_exponential, only: exp_pdf => pdf_exp, & + rexp => rvs_exp + + implicit none + real :: x(2,3,4),a(2,3,4) + complex :: scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, exp_pdf(1.0,1.0) !a probability density at 1.0 in standard expon + +! 0.367879450 + + print *, exp_pdf(2.0,2.0) !a probability density at 2.0 with lambda=2.0 + +! 3.66312787E-02 + + x = reshape(rexp(0.5, 24),[2,3,4]) ! standard expon random variates array + a(:,:,:) = 0.5 + print *, exp_pdf(x, a) ! a rank 3 standard expon probability density + +! 0.457115263 0.451488823 0.492391467 0.485233188 0.446215510 +! 0.401670188 0.485127628 0.316924453 0.418474048 0.483173639 +! 0.307366133 0.285812140 0.448017836 0.426440030 0.403896868 +! 0.334653258 0.410376132 0.485370994 0.333617479 0.263791025 +! 0.249779820 0.457159877 0.495636940 0.482243657 + + scale = (1.0, 2.) + print *, exp_pdf((1.5,1.0), scale) + ! a complex expon probability density function at (1.5,1.0) with real part + !of lambda=1.0 and imaginary part of lambda=2.0 + +! 6.03947677E-02 + +end program demo_exponential_pdf +``` + +## `cdf_exp` - exponential cumulative distribution function + +### Status + +Experimental + +### Description + +Cumulative distribution function (cdf) of the single real variable exponential distribution: + +$$F(x)=\begin{cases}1 - e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$ + +For a complex variable (x + y i) with independent real x and imaginary y parts, the joint cumulative distribution +function is the product of corresponding marginal cdf of real and imaginary cdf (for more details, see +"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): + +$$F(x+\mathit{i}y)=F(x)F(y)=\begin{cases} (1 - e^{-\lambda_{x} x})(1 - e^{-\lambda_{y} y}) &x\geqslant 0, \;\; y\geqslant 0 \\\\ 0 &otherwise \end{}$$ + +### Syntax + +`result = [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]](x, lambda)` + +### Class + +Elemental function + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. + +`lambda`: has `intent(in)` and is a scalar of type `real` or `complex`. + +All arguments must have the same type. + +### Return value + +The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments. + +### Example + +```fortran +program demo_exponential_cdf + use stdlib_random, only : random_seed + use stdlib_stats_distribution_exponential, only : exp_cdf => cdf_exp, & + rexp => rvs_exp + + implicit none + real :: x(2,3,4),a(2,3,4) + complex :: scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, exp_cdf(1.0, 1.0) ! a standard exponential cumulative at 1.0 + +! 0.632120550 + + print *, exp_cdf(2.0, 2.0) ! a cumulative at 2.0 with lambda=2 + +! 0.981684387 + + x = reshape(rexp(0.5, 24),[2,3,4]) + ! standard exponential random variates array + a(:,:,:) = 0.5 + print *, exp_cdf(x, a) ! a rank 3 array of standard exponential cumulative + +! 8.57694745E-02 9.70223546E-02 1.52170658E-02 2.95336246E-02 +! 0.107568979 0.196659625 2.97447443E-02 0.366151094 0.163051903 +! 3.36527228E-02 0.385267735 0.428375721 0.103964329 0.147119939 +! 0.192206264 0.330693483 0.179247737 2.92580128E-02 0.332765043 +! 0.472417951 0.500440359 8.56802464E-02 8.72612000E-03 3.55126858E-02 + + scale = (0.5,1.0) + print *, exp_cdf((0.5,0.5),scale) + !complex exponential cumulative distribution at (0.5,0.5) with real part of + !lambda=0.5 and imaginary part of lambda=1.0 + +! 8.70351046E-02 + +end program demo_exponential_cdf + +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8338033f2..f1f9051dd 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -39,6 +39,7 @@ set(fppFiles stdlib_stats_moment_scalar.fypp stdlib_stats_distribution_uniform.fypp stdlib_stats_distribution_normal.fypp + stdlib_stats_distribution_exponential.fypp stdlib_stats_var.fypp stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index c271b6aa6..3cb304d89 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -41,6 +41,7 @@ SRCFYPP = \ stdlib_stats_moment_scalar.fypp \ stdlib_stats_distribution_uniform.fypp \ stdlib_stats_distribution_normal.fypp \ + stdlib_stats_distribution_exponential.fypp \ stdlib_stats_var.fypp \ stdlib_math.fypp \ stdlib_math_linspace.fypp \ @@ -203,6 +204,11 @@ stdlib_stats_distribution_normal.o: \ stdlib_error.o \ stdlib_random.o \ stdlib_stats_distribution_uniform.o +stdlib_stats_distribution_exponential.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_random.o \ + stdlib_stats_distribution_uniform.o stdlib_random.o: \ stdlib_kinds.o \ stdlib_error.o diff --git a/src/stdlib_stats_distribution_exponential.fypp b/src/stdlib_stats_distribution_exponential.fypp new file mode 100644 index 000000000..01d4e8eb8 --- /dev/null +++ b/src/stdlib_stats_distribution_exponential.fypp @@ -0,0 +1,322 @@ +#: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 diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index b3b70bd4d..651286708 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -7,6 +7,7 @@ set(fppFiles test_median.fypp test_distribution_uniform.fypp test_distribution_normal.fypp + test_distribution_exponential.fypp ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -22,6 +23,7 @@ ADDTEST(varn) ADDTEST(random) ADDTEST(distribution_uniform) ADDTEST(distribution_normal) +ADDTEST(distribution_exponential) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index 15797d509..23c027dbf 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -2,7 +2,8 @@ SRCFYPP = \ test_mean.fypp \ test_median.fypp \ test_distribution_uniform.fypp \ - test_distribution_normal.fypp + test_distribution_normal.fypp \ + test_distribution_exponential.fypp SRCGEN = $(SRCFYPP:.fypp=.f90) diff --git a/src/tests/stats/test_distribution_exponential.fypp b/src/tests/stats/test_distribution_exponential.fypp new file mode 100644 index 000000000..95dc998c1 --- /dev/null +++ b/src/tests/stats/test_distribution_exponential.fypp @@ -0,0 +1,276 @@ + +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +program test_distribution_expon + use stdlib_kinds, only : sp, dp, xdp, qp + use stdlib_error, only : check + use stdlib_random, only : random_seed + use stdlib_stats_distribution_exponential, only : expon_rvs => rvs_exp, & + expon_pdf => pdf_exp, expon_cdf => cdf_exp + + implicit none + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: ${k1}$tol = 1000 * epsilon(1.0_${k1}$) + #:endfor + logical :: warn = .true. + integer :: put, get + + put = 12345678 + call random_seed(put, get) + + + call test_exponential_random_generator + + + #:for k1, t1 in RC_KINDS_TYPES + call test_expon_rvs_${t1[0]}$${k1}$ + #:endfor + + + #:for k1, t1 in RC_KINDS_TYPES + call test_expon_pdf_${t1[0]}$${k1}$ + #:endfor + call test_expon_pdf_rsp + + + #:for k1, t1 in RC_KINDS_TYPES + call test_expon_cdf_${t1[0]}$${k1}$ + #:endfor + + + + + +contains + + subroutine test_exponential_random_generator + + integer, parameter :: num = 10000000, array_size = 1000 + integer :: i, j, freq(0:array_size) + real(dp) :: chisq, expct + + print *, "" + print *, "Test exponential random generator with chi-squared" + freq = 0 + do i = 1, num + j = 1000 * (1 - exp(- expon_rvs(1.0))) + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + expct = num / array_size + do i = 0, array_size - 1 + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & + //" 1143.92" + write(*,*) "Chi-squared for exponential random generator is : ", chisq + call check((chisq < 1143.9), & + msg="exponential randomness failed chi-squared test", warn=warn) + end subroutine test_exponential_random_generator + + + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_expon_rvs_${t1[0]}$${k1}$ + ${t1}$ :: res(10), scale + integer, parameter :: k = 5 + integer :: i, n + integer :: seed, get + #:if t1[0] == "r" + #! for real type + ${t1}$, parameter :: ans(10) = & + [0.609680481289574416337018192280083895_${k1}$, & + 0.137541023585612635452927558314210417_${k1}$, & + 0.134921508232253721063879462841820585_${k1}$, & + 1.33766060689229752493171569464417802_${k1}$, & + 0.111148487576340881943792737729381770_${k1}$, & + 0.533951653963536868966836361020492979_${k1}$, & + 1.96897428558727671799033487332053483_${k1}$, & + 0.371111977992924465160247867364281152_${k1}$, & + 0.811918715695663687862785688290993341_${k1}$, & + 0.404637854946697868759504975362991277_${k1}$] + #:else + #! for complex type + ${t1}$, parameter :: ans(10) = & + [(1.30645817419194517786503898345732266_${k1}$, & + 0.158701181060322271676454874977935106_${k1}$), & + (0.289117517640543687994027420375329869_${k1}$, & + 1.54345454641418945184428733997405138_${k1}$), & + (0.238175330520730461308127295134389521_${k1}$, & + 0.616098062265619464192503493485184250_${k1}$), & + (4.21923061197273582426500329997257485_${k1}$, & + 0.428206128453374382877209077728016710_${k1}$), & + (1.73982581934785075970596933205212874_${k1}$, & + 0.466889832630805233184044202341912994_${k1}$), & + (2.22649889847873832288931745486999202_${k1}$, & + 0.879109337848515628785697537331053851_${k1}$), & + (8.76802198822945553859296653951917464_${k1}$, & + 0.200128045239398311139211728004738688_${k1}$), & + (0.694821947760945587572020290930855262_${k1}$, & + 0.101964167346166995492113143812345625_${k1}$), & + (0.141476585024528208770330398432893829_${k1}$, & + 3.989655879458742013468417133900891716E-0002_${k1}$), & + (2.10676792861163792685325850990401309_${k1}$, & + 0.249356813451327473065187125310051027_${k1}$)] + #:endif + + print *, "Test exponential_distribution_rvs_${t1[0]}$${k1}$" + seed = 593742186 + call random_seed(seed, get) + + #:if t1[0] == "r" + #! for real type + scale = 1.5_${k1}$ + #:else + #! for complex type + scale = (0.7_${k1}$, 1.3_${k1}$) + #:endif + + do i = 1, k + res(i) = expon_rvs(scale) ! 1 dummy + end do + res(6:10) = expon_rvs(scale, k) ! 2 dummies + call check(all(abs(res - ans) < ${k1}$tol), & + msg="exponential_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_expon_rvs_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_expon_pdf_${t1[0]}$${k1}$ + + ${t1}$ :: x1, x2(3,4), scale + integer :: i, n + integer :: seed, get + real(${k1}$) :: res(3,5) + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [0.362692289054629718342313806171796533_${k1}$, & + 0.362692289054629718342313806171796533_${k1}$, & + 0.362692289054629718342313806171796533_${k1}$, & + 1.44877092399186122284289290051705535_${k1}$, & + 1.08871761038277651996081144393335589_${k1}$, & + 0.203258408490339213767867275283195750_${k1}$, & + 0.730004225568590859263284124264208147_${k1}$, & + 0.237394827760488451509080387833146683_${k1}$, & + 0.301732182586179598102005265289645959_${k1}$, & + 1.35079274124711914255014934401469271_${k1}$, & + 0.416578245043239337295928202660090263_${k1}$, & + 1.44039177901335374382803898226703593_${k1}$, & + 0.196044829271295768265275728683411055_${k1}$, & + 0.271373826917613661285112379170965958_${k1}$, & + 1.00108987409617105109732206933052664_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [0.112097715784191810518066563334849515_${k1}$, & + 0.112097715784191810518066563334849515_${k1}$, & + 0.112097715784191810518066563334849515_${k1}$, & + 4.72087485401191174735651518020251204E-0002_${k1}$, & + 3.69705018439006691768174449531170720E-0002_${k1}$, & + 8.69498969681198520061798177185735738E-0002_${k1}$, & + 0.128007654288233028296342302153338001_${k1}$, & + 1.13496395875758374774198906169957218E-0002_${k1}$, & + 0.294260498264128747413785056084385424_${k1}$, & + 4.66169813179250908948018478030960097E-0002_${k1}$, & + 2.84438693906889813143446828488861951E-0002_${k1}$, & + 0.161859307815385236742977105439660254_${k1}$, & + 4.22904796362406579112752522035325397E-0002_${k1}$, & + 0.176117981883470250164040199296778089_${k1}$, & + 0.107352342201327219885025541854724060_${k1}$] + #:endif + + print *, "Test exponential_distribution_pdf_${t1[0]}$${k1}$" + seed = 123987654 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + scale = 1.5_${k1}$ + #:else + #! for complex type + scale = (0.3_${k1}$, 1.6_${k1}$) + #:endif + + x1 = expon_rvs(scale) + x2 = reshape(expon_rvs(scale, 12), [3,4]) + res(:,1) = expon_pdf(x1, scale) + res(:, 2:5) = expon_pdf(x2, scale) + call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), & + msg="exponential_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_expon_pdf_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_expon_cdf_${t1[0]}$${k1}$ + + ${t1}$ :: x1, x2(3,4), scale + integer :: i, n + integer :: seed, get + real(${k1}$) :: res(3,5) + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [0.109257742653704886153815776449785051_${k1}$, & + 0.109257742653704886153815776449785051_${k1}$, & + 0.109257742653704886153815776449785051_${k1}$, & + 0.717506371795765265795319089684216215_${k1}$, & + 6.82471795435370961628021592837348251E-0002_${k1}$, & + 0.158022297254037860379992220663140220_${k1}$, & + 0.914579543576380160727189390750289231_${k1}$, & + 0.735445094339121647068624074363021598_${k1}$, & + 8.69845458684957375690771394578441361E-0002_${k1}$, & + 0.491195342629961409581199224477971938_${k1}$, & + 0.574283568793105916250099261345264380_${k1}$, & + 0.312823040527767907760475800138803955_${k1}$, & + 0.640029783598040153827956625977856239_${k1}$, & + 2.16202116731629451897815202649346917E-0002_${k1}$, & + 7.74788145547936974757767867581111655E-0002_${k1}$] + #:else + real(${k1}$), parameter :: ans(15) = & + [7.83931265220552191922145459533155073E-0002_${k1}$, & + 7.83931265220552191922145459533155073E-0002_${k1}$, & + 7.83931265220552191922145459533155073E-0002_${k1}$, & + 1.07845760925785109085652212151328215E-0002_${k1}$, & + 0.672623038706161724678635394849362256_${k1}$, & + 4.27264038113873579678831482902258168E-0002_${k1}$, & + 0.179649132114996961326498233168917293_${k1}$, & + 1.38375793985183014482681114776428612E-0002_${k1}$, & + 3.49246365297941076158369468479748612E-0002_${k1}$, & + 0.116869945417176368845403154176734792_${k1}$, & + 0.468462732010133566674397830557697485_${k1}$, & + 0.413506985517976634907329948218002431_${k1}$, & + 0.665679674838121942273909342901808398_${k1}$, & + 0.223748595107983772617787558595393205_${k1}$, & + 0.337722969540396286456937689606849800_${k1}$] + #:endif + + print *, "Test exponential_distribution_cdf_${t1[0]}$${k1}$" + seed = 621957438 + call random_seed(seed, get) + + #:if t1[0] == "r" + #! for real type + scale = 2.0_${k1}$ + #:else + scale = (1.3_${k1}$, 2.1_${k1}$) + #:endif + + x1 = expon_rvs(scale) + x2 = reshape(expon_rvs(scale, 12), [3,4]) + res(:,1) = expon_cdf(x1, scale) + res(:, 2:5) = expon_cdf(x2, scale) + call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), & + msg="exponential_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_expon_cdf_${t1[0]}$${k1}$ + + #:endfor + +end program test_distribution_expon