diff --git a/doc/specs/index.md b/doc/specs/index.md index e4a5d8d76..a3b0a5def 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -24,6 +24,7 @@ This is and index/directory of the specifications (specs) for each new module/fe - [random](./stdlib_random.html) - Probability Distributions random number generator - [sorting](./stdlib_sorting.html) - Sorting of rank one arrays - [stats](./stdlib_stats.html) - Descriptive Statistics + - [stats_distributions_uniform](./stdlib_stats_distribution_uniform.html) - Uniform 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_uniform.md b/doc/specs/stdlib_stats_distribution_uniform.md new file mode 100644 index 000000000..485684b98 --- /dev/null +++ b/doc/specs/stdlib_stats_distribution_uniform.md @@ -0,0 +1,379 @@ +--- +title: stats_distribution_uniform +--- + +# Statistical Distributions -- Uniform Distribution Module + +[TOC] + +## `shuffle` - Using Fisher-Yates algorithm to generate a random permutation of a list + +### Status + +Experimental + +### Description + +Applying Fisher-Yates algorithm to generate an unbiased permutation for any list of intrinsic numerical data types. + +### Syntax + +`result = [[stdlib_stats_distribution_uniform(module):shuffle(interface)]]( list )` + +### Class + +Function. + +### Arguments + +`list`: argument has `intent(in)` and is a rank one array of `integer`, `real`, or `complex` type. + +### Return value + +Return a randomized rank one array of the input type. + +### Example + +```fortran +program demo_shuffle + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_uniform, only : shuffle + implicit none + integer :: seed_put, seed_get, i + real :: x(10) + integer :: n(10) + complex :: z(10) + + do i=1, 10 + n(i) = i + x(i) = real(i) + z(i) = cmplx(real(i), real(i)) + end do + seed_put = 32165498 + call random_seed(seed_put, seed_get) ! set and get current value of seed + print *, shuffle(n) ! get randomized n + +!10 6 9 2 8 1 3 5 7 4 + + print *, shuffle(x) ! get randomized x + +!5.0 10.0 9.0 4.0 3.0 8.0 2.0 1.0 7.0 6.0 + + print *, shuffle(z) ! get randomized z + +!(8.0, 8.0) (7.0, 7.0) (4.0, 4.0) (1.0, 1.0) (5.0, 5.0) +!(9.0, 9.0) (6.0, 6.0) (3.0, 3.0) (2.0, 2.0) (10.0, 10.0) + +end program demo_shuffle +``` + +## `rvs_uniform` - uniform distribution random variates + +### Status + +Experimental + +### Description + +Without argument the function returns a scalar standard uniformly distributed variate U(0,1) of `real` type with single precision on [0,1]. + +With single argument `scale` of `integer` type the function returns a scalar uniformly distributed variate of `integer` type on [0,scale]. This is the standard Rectangular distribution. + +With single argument `scale` of `real` or `complex` type the function returns a scalar uniformly distributed variate of `real` or `complex` type on [0, scale]. + +With double arguments `loc` and `scale` the function returns a scalar uniformly distributed random variates of `integer`, `real` or `complex` type on [loc, loc + scale] dependent of input type. + +With triple arguments `loc`, `scale` and `array_size` the function returns a rank one array of uniformly distributed variates of `integer`, `real` or `complex` type with an array size of `array_size`. + +For `complex` type, the real part and imaginary part are independent of each other. + +### Syntax + +`result = [[stdlib_stats_distribution_uniform(module):rvs_uniform(interface)]]([[loc,] scale] [[[,array_size]]])` + +### Class + +Elemental function (without the third argument). + +### Arguments + +`loc`: optional argument has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +`scale`: optional argument has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. + +`loc` and `scale` must have the same type and kind when both are present. + +### Return value + +The result is a scalar or a rank one array, with size of `array_size`, of type `integer`, `real` or `complex` depending on the input type. + +### Example + +```fortran +program demo_uniform_rvs + use stdlib_stats_distribution_PRNG, only:random_seed + use stdlib_stats_distribution_uniform, only:uni=> rvs_uniform + + implicit none + complex :: loc, scale + real :: a(3,4,5), b(3,4,5) + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, uni( ) !real standard uniform random variate in [0., 1.] +! 0.161520019 + + print *, uni(3.0) !an uniform random variate in [0., 3.] +! 1.65974522 + + print *, uni(-0.5, 1.0) !an uniform random variate in [-0.5, 0.5] +! 0.486900032 + + print *, uni(-1.0,2.0,10) + !an array of 10 uniform random variates in [-1., 1.] + +!0.884182811 -0.771520197 0.560377002 0.709313750 -7.12267756E-02 +!-0.431066573 0.497536063 -0.396331906 -0.325983286 0.137686729 + + print *, uni(20) !a random integer variate in [0, 20] +! 17 + + print *, uni(5,13) !a random integer variate in [5, 18] +! 15 + + print *, uni(3,19,10) !an array of 10 integer variates in [3,22] + +!7 16 16 12 9 21 19 4 3 19 + + loc = (-0.5, -0.5) + scale = (1.0, 1.0) + + print *, uni(scale) !a complex uniform random variate in unit square + +!(0.139202669, 0.361759573) + + print *, uni(loc,scale) + !a complex uniform random variate in [(-0.5, -0.5), (0.5, 0.5)] + +!(0.296536088,-0.143987954) + + print *, uni(loc, scale, 10) + !an array of 10 complex uniform random variate in [(-0.5, -0.5), (0.5, 0.5)] + +!(-0.302334785,-0.401923567) (0.281620383,9.534919262E-02) +! (-0.374348879,0.457528770) (0.442990601,-0.240510434) +! (-0.421572685,0.279313922) (-0.182090610,5.901372433E-02) +! (-7.864198089E-02,0.378484428) (-0.423258364,-0.201292425) +! (0.193327367,-0.353985727) (-0.397661150,0.355926156) + + a(:,:,:) = -0.5 + b(:,:,:) = 1.0 + + print *, uni(a,b) + !a rank 3 array of random variates in [-0.5,0.5] + +! -0.249188632 -0.199248433 -0.389813602 2.88307667E-03 0.238479793, +! 0.264856219 -0.205177426 -0.480921626 0.131218433 0.252170086, +! -0.303151041 -8.89462233E-02 -0.377370685 0.341802299 0.323204756, +! 0.358679056 -0.138909757 0.384329498 -0.109372199 0.132353067, +! 0.494320452 0.419343710 -0.103044361 0.461389005 0.403132677 +! 0.121850729 0.403839290 -0.349389791 0.490482628 0.156600773 +! 8.46788883E-02 -0.483680278 0.388107836 0.119698405 0.154214382 +! 0.153113484 0.236523747 0.155937552 -0.135760903 0.219589531 +! 0.394639254 6.30156994E-02 -0.342692465 -0.444846451 -0.215700030 +! 0.204189956 -0.208748132 0.355063021 8.98272395E-02 -0.237928331 +! 2.98077464E-02 -0.485149682 -8.06870461E-02 -0.372713923 +! -0.178335011 0.283877611 -2.13934183E-02 -9.21690464E-03 +! 4.56320047E-02 0.220112979 + +end program demo_uniform_rvs +``` + +## `pdf_uniform` - Uniform probability density function + +### Status + +Experimental + +### Description + +The probability density function of the uniform distribution. + +f(x) = 1 / (scale + 1); for discrete uniform distribution. + +f(x) = 1 / scale; for continuous uniform distribution. + +f(x) = 1 / (scale%re * scale%im); for complex uniform distribution. + +### Syntax + +`result = [[stdlib_stats_distribution_uniform(module):pdf_uniform(interface)]](x, loc, scale)` + +### Class + +Elemental function. + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +`loc`: has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +`scale`: has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +All three arguments must have the same type and kind. + +### Return value + +The result is a scalar or an array, with a shape conformable to arguments, of type `real`. + +### Example + +```fortran +program demo_uniform_pdf + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_uniform, only : uni_pdf => pdf_uniform, & + uni => rvs_uniform + + implicit none + complex :: loc, scale + real :: a(3,4,5), b(3,4,5), x(3,4,5) + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, uni_pdf(3, 2, 10) !probability density at 3 in range [2, 10] + +! 9.09090936E-02 + + print *, uni_pdf(0.5,0.0,1.0) !a probability density at 0.5 in [0., 1.] + +! 1.00000000 + + + print *, uni_pdf(0.7,-1.0,2.0) !a probability density at 0.7 in [-1., 1.] + +! 0.500000000 + + a(:,:,:) = 0.0 + b(:,:,:) = 2.0 + x = reshape(uni(0.,2.,60),[3,4,5])! uniform random variates array in [0., 2.] + print *, uni_pdf(x, a, b) ! probability density array in [0., 2.] + +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 +! 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 0.500000000 + + loc = (-0.5,-0.5) + scale = (1.0,1.0) + print *, uni_pdf((-0.1,0.2), loc, scale) + ! joint probability density at (-0.1,0.2) in [(-0.5, -0.5), (0.5, 0.5)] + +! 1.00000000 +end program demo_uniform_pdf + +``` + +## `cdf_uniform` - Uniform cumulative distribution function + +### Status + +Experimental + +### Description + +Cumulative distribution function of the uniform distribution + +F(x) = (x - loc + 1) / (scale + 1); for discrete uniform distribution. + +F(x) = (x - loc) / scale; for continuous uniform distribution. + +F(x) = (x%re - loc%re)(x%im - loc%im) / (scale%re * scale%im); for complex uniform distribution. + +### Syntax + +`result = [[stdlib_stats_distribution_uniform(module):cdf_uniform(interface)]](x, loc, scale)` + +### Class + +Elemental function. + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +`loc`: has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +`scale`: has `intent(in)` and is a scalar of type `integer`, `real` or `complex`. + +All three arguments must have the same type and kind. + +### Return value + +The result is a scalar or an array, with a shape conformable to arguments, of type `real`. + +### Example + +```fortran +program demo_uniform_cdf + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_uniform, only : uni_cdf => cdf_uniform, & + uni => rvs_uniform + + implicit none + real :: x(3,4,5), a(3,4,5), b(3,4,5) + complex :: loc, scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, uni_cdf(0.5,0.,1.) ! a cumulative at 0.5 in [0., 1.] + +!0.500000000 + + print *, uni_cdf(0.7,-1.0,2.0) ! a cumulative at 0.7 in [-1.0, 1.0] + +! 0.850000024 + + print *, uni_cdf(6, 2, 10) ! a cumulative at 6 in [2, 10] + +! 0.454545468 + + a(:,:,:) = -1.0 + b(:,:,:) = 2.0 + x = reshape(uni(-1.0,2.0,60),[3,4,5]) ! uniform random variates array + print *, uni_cdf(x,a,b) ! cumulative array in [-1.0, 1.0] + +!0.161520004 0.553248405 0.986900032 0.942091405 0.114239901 0.780188501 +! 0.854656875 0.464386612 0.284466714 0.748768032 0.301834047 0.337008357 +!0.568843365 0.596165061 0.180993259 0.614166319 0.214835495 7.98164606E-02 +!0.641274095 0.607101977 0.701139212 0.230517209 1.97925568E-02 0.857982159 +!0.712761045 0.139202654 0.361759573 0.796536088 0.356012046 0.197665215 +!9.80764329E-02 0.781620383 0.595349193 0.125651121 0.957528770 0.942990601 +!0.259489566 7.84273148E-02 0.779313922 0.317909390 0.559013724 0.421358019 +!0.878484428 7.67416358E-02 0.298707575 0.693327367 0.146014273 0.102338850 +!0.855926156 0.250811368 0.300751567 0.110186398 0.502883077 0.738479793 +!0.764856219 0.294822574 1.90783739E-02 0.631218433 0.752170086 0.196848959 + + loc = (0., 0.) + scale=(2., 1.) + print *, uni_cdf((1.2,0.5), loc, scale) + ! joint cumulative distribution at (1.2,0.5) in [(0.,0.), (2.,1.)] + +! 0.300000012 +end program demo_uniform_cdf + +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 520e105f6..bb9fb4fd8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -### Pre-process: .fpp -> .f90 via Fypp +#### Pre-process: .fpp -> .f90 via Fypp # Create a list of the files to be preprocessed set(fppFiles @@ -14,7 +14,7 @@ set(fppFiles stdlib_sorting.fypp stdlib_sorting_ord_sort.fypp stdlib_sorting_sort.fypp - stdlib_sorting_sort_index.fypp + stdlib_sorting_sort_index.fypp stdlib_stats.fypp stdlib_stats_corr.fypp stdlib_stats_cov.fypp @@ -24,6 +24,7 @@ set(fppFiles stdlib_stats_moment_all.fypp stdlib_stats_moment_mask.fypp stdlib_stats_moment_scalar.fypp + stdlib_stats_distribution_uniform.fypp stdlib_stats_var.fypp stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index 0bf5ef02d..179fc600f 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -26,6 +26,7 @@ SRCFYPP = \ stdlib_stats_moment_all.fypp \ stdlib_stats_moment_mask.fypp \ stdlib_stats_moment_scalar.fypp \ + stdlib_stats_distribution_uniform.fypp \ stdlib_stats_var.fypp \ stdlib_math.fypp \ stdlib_math_linspace.fypp \ @@ -145,6 +146,10 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_random.o stdlib_random.o: \ stdlib_kinds.o \ stdlib_error.o diff --git a/src/stdlib_stats_distribution_uniform.fypp b/src/stdlib_stats_distribution_uniform.fypp new file mode 100644 index 000000000..7c441a0b9 --- /dev/null +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -0,0 +1,518 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set ALL_KINDS_TYPES = INT_KINDS_TYPES + RC_KINDS_TYPES +module stdlib_stats_distribution_uniform + use stdlib_kinds + use stdlib_error, only : error_stop + use stdlib_random, only : dist_rand + + implicit none + private + + real(dp), parameter :: MESENNE_NUMBER = 1.0_dp / (2.0_dp ** 53 - 1.0_dp) + integer(int64), parameter :: INT_ONE = 1_int64 + + public :: rvs_uniform + public :: pdf_uniform + public :: cdf_uniform + public :: shuffle + + + interface rvs_uniform + !! version: experimental + !! + !! Get uniformly distributed random variate for integer, real and complex + !! variables. + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! rvs_uniform-uniform-distribution-random-variates)) + + module procedure rvs_unif_0_rsp ! 0 dummy variable + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure rvs_unif_1_${t1[0]}$${k1}$ ! 1 dummy variable + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure rvs_unif_${t1[0]}$${k1}$ ! 2 dummy variables + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure rvs_unif_array_${t1[0]}$${k1}$ ! 3 dummy variables + #:endfor + end interface rvs_uniform + + + interface pdf_uniform + !! version: experimental + !! + !! Get uniform distribution probability density (pdf) for integer, real and + !! complex variables. + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! pdf_uniform-uniform-probability-density-function)) + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure pdf_unif_${t1[0]}$${k1}$ + #:endfor + end interface pdf_uniform + + + interface cdf_uniform + !! version: experimental + !! + !! Get uniform distribution cumulative distribution function (cdf) for integer, + !! real and complex variables. + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! cdf_uniform-uniform-cumulative-distribution-function)) + !! + #:for k1, t1 in ALL_KINDS_TYPES + module procedure cdf_unif_${t1[0]}$${k1}$ + #:endfor + end interface cdf_uniform + + + interface shuffle + !! version: experimental + !! + !! Fisher-Yates shuffle algorithm for a rank one array of integer, real and + !! complex variables. + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! shuffle-using-fisher-yates-algorithm-to-generate-a-random-permutation-of-a-list)) + !! + #:for k1, t1 in ALL_KINDS_TYPES + module procedure shuffle_${t1[0]}$${k1}$ + #:endfor + end interface shuffle + + + + + + +contains + + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function rvs_unif_1_${t1[0]}$${k1}$(scale) result(res) + ! + ! Uniformly distributed integer in [0, scale] + ! Bitmask with rejection + ! https://www.pcg-random.org/posts/bounded-rands.html + ! + ! Fortran 90 translated from c by Jim-215-fisher + ! + ${t1}$, intent(in) :: scale + ${t1}$ :: res, u, mask + integer :: zeros, bits_left, bits + + if(scale <= 0_${k1}$) call error_stop("Error(rvs_unif_1): Uniform" & + //" distribution scale parameter must be positive") + zeros = leadz(scale) + bits = bit_size(scale) - zeros + mask = shiftr(not(0_${k1}$), zeros) + L1 : do + u = dist_rand(scale) + res = iand(u, mask) + if(res <= scale) exit L1 + bits_left = zeros + L2 : do + if(bits_left < bits) exit L2 + u = shiftr(u, bits) + res = iand(u, mask) + if(res <= scale) exit L1 + bits_left = bits_left - bits + end do L2 + end do L1 + end function rvs_unif_1_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function rvs_unif_${t1[0]}$${k1}$(loc, scale) result(res) + ! + ! Uniformly distributed integer in [loc, loc + scale] + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + + if(scale <= 0_${k1}$) call error_stop("Error(rvs_unif): Uniform" & + //" distribution scale parameter must be positive") + res = loc + rvs_unif_1_${t1[0]}$${k1}$(scale) + end function rvs_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function rvs_unif_0_${t1[0]}$${k1}$( ) result(res) + ! + ! Uniformly distributed float in [0,1] + ! Based on the paper by Frederic Goualard, "Generating Random Floating- + ! Point Numbers By Dividing Integers: a Case Study", Proceedings of + ! ICCS 2020, June 2020, Amsterdam, Netherlands + ! + ${t1}$ :: res + integer(int64) :: tmp + + tmp = shiftr(dist_rand(INT_ONE), 11) ! Get random from [0,2^53-1] + res = real(tmp * MESENNE_NUMBER, kind = ${k1}$) ! convert to [0,1] + end function rvs_unif_0_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function rvs_unif_1_${t1[0]}$${k1}$(scale) result(res) + ! + ! Uniformly distributed float in [0, scale] + ! + ${t1}$, intent(in) :: scale + ${t1}$ :: res + + if(scale == 0._${k1}$) call error_stop("Error(rvs_unif_1): " & + //"Uniform distribution scale parameter must be non-zero") + res = scale * rvs_unif_0_${t1[0]}$${k1}$( ) + end function rvs_unif_1_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function rvs_unif_${t1[0]}$${k1}$(loc, scale) result(res) + ! + ! Uniformly distributed float in [loc, loc + scale] + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + + if(scale == 0._${k1}$) call error_stop("Error(rvs_unif): " & + //"Uniform distribution scale parameter must be non-zero") + res = loc + scale * rvs_unif_0_${t1[0]}$${k1}$( ) + end function rvs_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function rvs_unif_1_${t1[0]}$${k1}$(scale) result(res) + ! + ! Uniformly distributed complex in [(0,0i), (scale, i(scale))] + ! The real part and imaginary part are independent of each other, so that + ! the joint distribution is on an unit square [(0,0i), (scale,i(scale))] + ! + ${t1}$, intent(in) :: scale + ${t1}$ :: res + real(${k1}$) :: r1, tr, ti + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(rvs_uni_" & + //"1): Uniform distribution scale parameter must be non-zero") + r1 = rvs_unif_0_r${k1}$( ) + if(scale % re == 0.0_${k1}$) then + ti = scale % im * r1 + tr = 0.0_${k1}$ + else if(scale % im == 0.0_${k1}$) then + tr = scale % re * r1 + ti = 0.0_${k1}$ + else + tr = scale % re * r1 + r1 = rvs_unif_0_r${k1}$( ) + ti = scale % im * r1 + end if + res = cmplx(tr, ti, kind=${k1}$) + end function rvs_unif_1_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function rvs_unif_${t1[0]}$${k1}$(loc, scale) result(res) + ! + ! Uniformly distributed complex in [(loc,iloc), (loc + scale, i(loc + + ! scale))]. + ! The real part and imaginary part are independent of each other, so that + ! the joint distribution is on an unit square [(loc,iloc), (loc + scale, + ! i(loc + scale))] + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + real(${k1}$) :: r1, tr, ti + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(rvs_uni_" & + //"): Uniform distribution scale parameter must be non-zero") + r1 = rvs_unif_0_r${k1}$( ) + if(scale % re == 0.0_${k1}$) then + tr = loc % re + ti = loc % im + scale % im * r1 + else if(scale % im == 0.0_${k1}$) then + tr = loc % re + scale % re * r1 + ti = loc % im + else + tr = loc % re + scale % re * r1 + r1 = rvs_unif_0_r${k1}$( ) + ti = loc % im + scale % im * r1 + end if + res = cmplx(tr, ti, kind=${k1}$) + end function rvs_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + function rvs_unif_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) + + integer, intent(in) :: array_size + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res(array_size) + ${t1}$ :: u, mask, nn + integer :: i, zeros, bits_left, bits + + if(scale == 0_${k1}$) call error_stop("Error(rvs_unif_array): " & + //"Uniform distribution scale parameter must be non-zero") + zeros = leadz(scale) + bits = bit_size(scale) - zeros + mask = shiftr(not(0_${k1}$), zeros) + do i = 1, array_size + L1 : do + u = dist_rand(scale) + nn = iand(u, mask) + if(nn <= scale) exit L1 + bits_left = zeros + L2 : do + if(bits_left < bits) exit L2 + u = shiftr(u, bits) + nn = iand(u, mask) + if(nn <= scale) exit L1 + bits_left = bits_left - bits + end do L2 + end do L1 + res(i) = loc + nn + end do + end function rvs_unif_array_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + function rvs_unif_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) + + integer, intent(in) :: array_size + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res(array_size) + ${t1}$ :: t + integer(int64) :: tmp + integer :: i + + + if(scale == 0._${k1}$) call error_stop("Error(rvs_unif_array):" & + //" Uniform distribution scale parameter must be non-zero") + do i = 1, array_size + tmp = shiftr(dist_rand(INT_ONE), 11) + t = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + res(i) = loc + scale * t + end do + end function rvs_unif_array_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + function rvs_unif_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) + + integer, intent(in) :: array_size + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res(array_size) + real(${k1}$) :: r1, tr, ti + integer(int64) :: tmp + integer :: i + + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(rvs_unif" & + //"_array): Uniform distribution scale parameter must be non-zero") + do i = 1, array_size + tmp = shiftr(dist_rand(INT_ONE), 11) + r1 = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + if(scale % re == 0.0_${k1}$) then + tr = loc % re + ti = loc % im + scale % im * r1 + else if(scale % im == 0.0_${k1}$) then + tr = loc % re + scale % re * r1 + ti = loc % im + else + tr = loc % re + scale % re * r1 + tmp = shiftr(dist_rand(INT_ONE), 11) + r1 = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + ti = loc % im + scale % im * r1 + end if + res(i) = cmplx(tr, ti, kind=${k1}$) + end do + end function rvs_unif_array_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res) + + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0_${k1}$) then + res = 0.0 + else if(x < loc .or. x > (loc + scale)) then + res = 0.0 + else + res = 1. / (scale + 1_${k1}$) + end if + end function pdf_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res) + + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0.0_${k1}$) then + res = 0.0 + else if(x < loc .or. x > (loc + scale)) then + res = 0.0 + else + res = 1.0 / scale + end if + end function pdf_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res) + + ${t1}$, intent(in) :: x, loc, scale + real :: res + real(${k1}$) :: tr, ti + + tr = loc % re + scale % re; ti = loc % im + scale % im + if(scale == (0.0_${k1}$,0.0_${k1}$)) then + res = 0.0 + else if((x % re >= loc % re .and. x % re <= tr) .and. & + (x % im >= loc % im .and. x % im <= ti)) then + res = 1.0 / (scale % re * scale % im) + else + res = 0.0 + end if + end function pdf_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res) + + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0_${k1}$) then + res = 0.0 + else if(x < loc) then + res = 0.0 + else if(x >= loc .and. x <= (loc + scale)) then + res = real((x - loc + 1_${k1}$)) / real((scale + 1_${k1}$)) + else + res = 1.0 + end if + end function cdf_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res) + + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0.0_${k1}$) then + res = 0.0 + else if(x < loc) then + res = 0.0 + else if(x >= loc .and. x <= (loc + scale)) then + res = (x - loc) / scale + else + res = 1.0 + end if + end function cdf_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res) + + ${t1}$, intent(in) :: x, loc, scale + real :: res + logical :: r1, r2, i1, i2 + + if(scale == (0.0_${k1}$,0.0_${k1}$)) then + res = 0.0 + return + end if + r1 = x % re < loc % re + r2 = x % re > (loc % re + scale % re) + i1 = x % im < loc % im + i2 = x % im > (loc % im + scale % im) + if(r1 .or. i1) then + res = 0.0 + else if((.not. r1) .and. (.not. r2) .and. i2) then + res = (x % re - loc % re) / scale % re + else if((.not. i1) .and. (.not. i2) .and. r2) then + res = (x % im - loc % im) / scale % im + else if((.not. r1) .and. (.not. r2) .and. (.not. i1) .and. (.not. i2)) & + then + res = (x % re - loc % re) * (x % im - loc % im) / & + (scale % re * scale % im) + else if(r2 .and. i2)then + res = 1.0 + end if + end function cdf_unif_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in ALL_KINDS_TYPES + function shuffle_${t1[0]}$${k1}$( list ) result(res) + + ${t1}$, intent(in) :: list(:) + ${t1}$ :: res(size(list)) + ${t1}$ :: tmp + integer :: n, i, j + + n = size(list) + res = list + do i = 1, n - 1 + j = rvs_uniform(n - i) + i + tmp = res(i) + res(i) = res(j) + res(j) = tmp + end do + end function shuffle_${t1[0]}$${k1}$ + + #:endfor +end module stdlib_stats_distribution_uniform diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 600e925e5..14fed4167 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -3,6 +3,7 @@ # Create a list of the files to be preprocessed set(fppFiles test_median.fypp + test_distribution_uniform.fypp ) # Custom preprocessor flags @@ -25,6 +26,7 @@ ADDTEST(rawmoment) ADDTEST(var) ADDTEST(varn) ADDTEST(distribution_PRNG) +ADDTEST(distribution_uniform) 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 55a21c5f4..515b6c113 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,5 +1,6 @@ SRCFYPP =\ - test_median.fypp + test_median.fypp \ + test_distribution_uniform.fypp SRCGEN = $(SRCFYPP:.fypp=.f90) diff --git a/src/tests/stats/common.fypp b/src/tests/stats/common.fypp index 9ebaac6b9..b8a51460f 100644 --- a/src/tests/stats/common.fypp +++ b/src/tests/stats/common.fypp @@ -18,6 +18,14 @@ #! Collected (kind, type) tuples for integer types #:set INT_KINDS_TYPES = list(zip(INT_KINDS, INT_TYPES)) +#! Complex kinds to be considered during templating +#:set CMPLX_KINDS = ["sp", "dp", "qp"] + +#! Complex types to be considered during templating +#:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS] + +#! Collected (kind, type) tuples for complex types +#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES)) #! Whether Fortran 90 compatible code should be generated #:set VERSION90 = defined('VERSION90') diff --git a/src/tests/stats/test_distribution_uniform.fypp b/src/tests/stats/test_distribution_uniform.fypp new file mode 100644 index 000000000..738f89bfe --- /dev/null +++ b/src/tests/stats/test_distribution_uniform.fypp @@ -0,0 +1,316 @@ +#:include "common.fypp" +#:set ALL_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +program test_distribution_uniform + use stdlib_error, only : check + use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, qp + use stdlib_random, only : random_seed, dist_rand + use stdlib_stats_distribution_uniform, uni_rvs => rvs_uniform, & + uni_pdf => pdf_uniform, & + uni_cdf => cdf_uniform + + implicit none + logical :: warn = .true. + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: ${k1}$tol = 1000 * epsilon(1.0_${k1}$) + #:endfor + integer :: put, get + + put = 135792468 + + call test_shuffle + + call test_uni_rvs_0 + + #:for k1, t1 in ALL_KINDS_TYPES + call test_uni_rvs_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + call test_uni_pdf_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + call test_uni_cdf_${t1[0]}$${k1}$ + #:endfor + + + + +contains + + subroutine test_shuffle + integer :: n(10) + integer, parameter :: na(10) = [10, 6, 9, 2, 8, 1, 3, 5, 7, 4] + real :: x(10) + real, parameter :: xa(10) = [5.0, 10.0, 9.0, 4.0, 3.0, 8.0, 2.0, 1.0, & + 7.0, 6.0] + complex :: z(10) + complex, parameter :: za(10) = [(8.0, 8.0), (7.0, 7.0), (4.0, 4.0), & + (1.0, 1.0), (5.0, 5.0), (9.0, 9.0), & + (6.0, 6.0), (3.0, 3.0), (2.0, 2.0), & + (10.0, 10.0)] + integer :: i, put, get + + do i = 1, 10 + n(i) = i + x(i) = real(i) + z(i) = cmplx(real(i),real(i)) + end do + put = 32165498 + call random_seed(put, get) + n(:) = shuffle(n) + x(:) = shuffle(x) + z(:) = shuffle(z) + call check(all(n == na), & + msg = "Integer shuffle failed test", warn=warn) + call check(all(x == xa), & + msg = "Real shuffle failed test", warn=warn) + call check(all(z == za), & + msg = "Complex shuffle failed test", warn=warn) + end subroutine test_shuffle + + + + subroutine test_uni_rvs_0 + integer, parameter :: num = 10000000 + integer, parameter :: array_size = 1000 + integer :: i, j, freq(0 : array_size - 1) + real(dp) :: chisq, expct + + print *,"" + print *, "Test uniform random generator with chi-squared" + freq = 0 + do i = 1, num + j = array_size * uni_rvs( ) + 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 value for chi-squared with 1000 d. of f. is" & + //" 1143.92" + write(*,*) "Chi-squared for uniform random generator is : ", chisq + call check((chisq < 1143.9) , & + msg = "uniform randomness failed chi-squared test", warn=warn) + end subroutine test_uni_rvs_0 + + + + #:for k1, t1 in ALL_KINDS_TYPES + subroutine test_uni_rvs_${t1[0]}$${k1}$ + ${t1}$ :: res(15), scale, loc + integer :: i, n, seed, get, k + + #:if k1 == "int8" + ${t1}$, parameter :: ans(15) = [47, 99, 43, 37, 48, 30, 27, 100, 30, & + 33, 21, 103, 55, 54, 110] + #:elif k1 == "int16" + ${t1}$, parameter :: ans(15) = [25, 4, 81, 98, 49, 34, 32, 62, 115, & + 112, 26, 20, 37, 100, 82] + #:elif k1 == "int32" + ${t1}$, parameter :: ans(15) = [19, 52, 56, 20, 59, 44, 34, 102, 19, & + 39, 60, 50, 97, 56, 67] + #:elif k1 == "int64" + ${t1}$, parameter :: ans(15) = [76, 45, 43, 75, 76, 15, 25, 24, 114, & + 113, 94, 29, 109, 93, 89] + #:elif t1[0] == "r" + #! for real type + ${t1}$, parameter :: ans(15) = & + [0.914826628538749186958511927514337003_${k1}$, & + 0.367330098664966409049981166390352882_${k1}$, & + 1.77591243057709280428468900936422870_${k1}$, & + 0.885921308987590139238932351872790605_${k1}$, & + 0.950735656542987861428173346212133765_${k1}$, & + -0.659562573857055134407545438079978339_${k1}$, & + -0.116661718506947176265953203255776316_${k1}$, & + 0.837391893893859151631886561517603695_${k1}$, & + -0.703954396598600540269075054311542772_${k1}$, & + 0.382592729851141566399519433616660535_${k1}$, & + -0.132472493978185168472805344208609313_${k1}$, & + -0.878723366294216184924081858298450243_${k1}$, & + -0.901660046141515819639877804547722917_${k1}$, & + -0.164090614147737401395943379611708224_${k1}$, & + -0.333886718190384290672056977200554684_${k1}$] + #:else + #! for complex type + ${t1}$, parameter :: ans(15) = & + [(0.457413314269374593479255963757168502_${k1}$, & + 0.183665049332483204524990583195176441_${k1}$), & + (0.887956215288546402142344504682114348_${k1}$, & + 0.442960654493795069619466175936395302_${k1}$), & + (0.475367828271493930714086673106066883_${k1}$, & + 0.170218713071472432796227280960010830_${k1}$), & + (0.441669140746526411867023398372111842_${k1}$, & + 0.918695946946929575815943280758801848_${k1}$), & + (0.148022801700699729865462472844228614_${k1}$, & + 0.691296364925570783199759716808330268_${k1}$), & + (-6.623624698909258423640267210430465639E-0002_${k1}$, & + 0.560638316852891907537959070850774879_${k1}$), & + (-0.450830023070757909819938902273861459_${k1}$, & + 0.917954692926131299302028310194145888_${k1}$), & + (-0.166943359095192145336028488600277342_${k1}$, & + 1.05997401970850635422038976685144007_${k1}$), & + (-0.429652190199228276035192664039641386_${k1}$, & + 0.523558274341032421628217008446881664_${k1}$), & + (0.427181091476823815433760955784237012_${k1}$, & + 1.34628934976074521312483511792379431_${k1}$), & + (-0.343281426018765739582860874179459643_${k1}$, & + 1.15357331316264255516301773241139017_${k1}$), & + (-0.127590074749816595467422075671493076_${k1}$, & + 1.06891199479835175001340985545539297_${k1}$), & + (0.262287586904722758163188700564205647_${k1}$, & + 1.29508919831907332032017166056903079_${k1}$), & + (-0.192677407376582732201342196276527829_${k1}$, & + 1.32794925614337933073016984053538181_${k1}$), & + (-0.264742129752461530234342035328154452_${k1}$, & + 1.01282963412172621886497836385387927_${k1}$)] + #:endif + + print *, "Test rvs_uniform_${t1[0]}$${k1}$" + seed = 258147369; k = 5 + call random_seed(seed, get) + #:if t1[0] == "i" + #! for integer type + loc = 15_${k1}$; scale = 100_${k1}$ + #:elif t1[0] == "r" + #! for real type + loc = -1.0_${k1}$; scale = 2.0_${k1}$ + #:else + #! for complex type + loc = (-0.5_${k1}$,0.5_${k1}$); scale = (1.0_${k1}$, 1.0_${k1}$) + #:endif + do i = 1, 5 + res(i) = uni_rvs(scale) ! 1 dummy + end do + do i = 6,10 + res(i) = uni_rvs(loc, scale) ! 2 dummies + end do + res(11:15) = uni_rvs(loc, scale, k) ! 3 dummies + #:if t1[0] == "i" + #! for integer type + call check(all(res == ans), & + msg="rvs_uniform_${t1[0]}$${k1}$ failed", warn=warn) + #:else + #! for real and complex types + call check(all(abs(res - ans) < ${k1}$tol), & + msg="rvs_uniform_${t1[0]}$${k1}$ failed", warn=warn) + #:endif + end subroutine test_uni_rvs_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in ALL_KINDS_TYPES + subroutine test_uni_pdf_${t1[0]}$${k1}$ + ${t1}$ :: x1, x2(3,4), loc, scale + integer :: seed, get, i + real :: res(3,5) + #:if t1[0] == "i" + #! for integer type + real, parameter :: ans(15) = [(1.96078438E-02, i=1,15)] + #:elif t1[0] == "r" + #! for real type + real, parameter :: ans(15) = [(0.5, i=1,15)] + #:else + #! for complex type + real, parameter :: ans(15) = [(1.0, i=1,15)] + #:endif + + print *, "Test pdf_uniform_${t1[0]}$${k1}$" + seed = 147258639 + call random_seed(seed, get) + #:if t1[0] == "i" + #! for integer type + loc = 0_${k1}$; scale = 50_${k1}$ + #:elif t1[0] == "r" + #! for real type + loc = 0.0_${k1}$; scale = 2.0_${k1}$ + #:else + #! for complex type + loc = (-0.5_${k1}$, 0.5_${k1}$); scale = (1.0_${k1}$, 1.0_${k1}$) + #:endif + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_pdf(x1, loc, scale) + res(:, 2:5) = uni_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg = "pdf_uniform_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_uni_pdf_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in ALL_KINDS_TYPES + subroutine test_uni_cdf_${t1[0]}$${k1}$ + ${t1}$ :: x1, x2(3,4), loc, scale + real :: res(3,5) + integer :: seed, get + #:if k1 == "int8" + real, parameter :: ans(15) = [0.435643554, 0.435643554, 0.435643554, & + 0.702970326, 0.653465331, 0.485148519, & + 0.386138618, 0.386138618, 0.336633652, & + 0.277227730, 0.237623766, 0.524752498, & + 0.732673287, 0.534653485, 0.415841579] + #:elif k1 == "int16" + real, parameter :: ans(15) = [0.178217828, 0.178217828, 0.178217828, & + 0.465346545, 0.673267305, 0.247524753, & + 0.158415839, 0.792079210, 0.742574275, & + 0.574257433, 0.881188095, 0.663366318, & + 0.524752498, 0.623762369, 0.178217828] + #:elif k1 == "int32" + real, parameter :: ans(15) = [0.732673287, 0.732673287, 0.732673287, & + 0.722772300, 0.792079210, 5.94059415E-02,& + 0.841584146, 0.405940592, 0.960396051, & + 0.534653485, 0.782178223, 0.861386120, & + 0.564356446, 0.613861382, 0.306930691] + #:elif k1 == "int64" + real, parameter :: ans(15) = [0.455445558, 0.455445558, 0.455445558, & + 0.277227730, 0.455445558, 0.930693090, & + 0.851485133, 0.623762369, 5.94059415E-02,& + 0.693069279, 0.544554472, 0.207920790, & + 0.306930691, 0.356435657, 0.128712878] + #:elif t1[0] == "r" + #! for real type + real, parameter :: ans(15) = [0.170192942, 0.170192942, 0.170192942, & + 0.276106149, 0.754930079, 0.406620681, & + 0.187742814, 0.651605546, 0.934733927, & + 0.151271492, 0.987674534, 0.130533904, & + 0.106271908, 9.27578658E-02, 0.203399420] + #:else + #! for complex type + real, parameter :: ans(15) = [4.69913185E-02, 4.69913185E-02, & + 4.69913185E-02, 0.306970179, 0.122334257,& + 0.141398594, 0.128925011, 9.85755492E-03,& + 8.16527531E-02, 0.163921610, & + 7.81712309E-02, 0.446415812, & + 5.31753292E-04, 0.101455867, 0.155276477] + #:endif + + print *, "Test cdf_uniform_${t1[0]}$${k1}$" + seed = 369147258 + call random_seed(seed, get) + #:if t1[0] == "i" + #! for integer type + loc = 14_${k1}$; scale = 100_${k1}$ + #:elif t1[0] == "r" + #! for real type + loc = 0.0_${k1}$; scale = 2.0_${k1}$ + #:else + #! for complex type + loc = (-0.5_${k1}$, -0.5_${k1}$); scale = (1.0_${k1}$, 1.0_${k1}$) + #:endif + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_cdf(x1, loc, scale) + res(:, 2:5) = uni_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg = "cdf_uniform_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_uni_cdf_${t1[0]}$${k1}$ + + #:endfor + +end program test_distribution_uniform