From edf8e54297139171d7a96925558ca3ccc4917709 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 18 Oct 2020 13:32:28 -0400 Subject: [PATCH 01/39] change in Makefile.manual changed the object name from stdlib_stats_distribution_implementation.o to stdlib_stats_distribution_imp.o --- src/Makefile.manual | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 1c731b9bb..129dbe553 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -12,7 +12,10 @@ SRC = f18estop.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ stdlib_stats_moment.f90 \ - stdlib_stats_var.f90 + stdlib_stats_var.f90 \ + stdlib_stats_distribution.f90 \ + stdlib_stats_distribution_rvs.f90 \ + stdlib_stats_distribution_implementation.f90 \ LIB = libstdlib.a @@ -61,7 +64,13 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o - +stdlib_stats_distribution_rvs.o: stdlib_kinds.o +stdlib_stats_distribution.o: \ + stdlib_error.o \ + stdlib_kinds.o \ + stdlib_stats_distribution_rvs.o \ +stdlib_stats_distribution_imp.o: \ + stdlib_stats_distribution.o # Fortran sources that are built from fypp templates stdlib_io.f90: stdlib_io.fypp stdlib_linalg.f90: stdlib_linalg.fypp @@ -71,3 +80,7 @@ stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp +stdlib_stats_distribution_rvs.f90: stdlib_stats_distribution_rvs.fypp +stdlib_stats_distribution.f90: stdlib_stats_distribution.fypp +stdlib_stats_distribution_implementation.f90: \ + stdlib_stats_distribution_implementation.fypp \ No newline at end of file From 2a495abc6144d8d8cec86d82c838578dbcf10f83 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 19:10:52 -0500 Subject: [PATCH 02/39] Add files via upload --- CMakeLists.txt | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6143b5ab5..d62689913 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,7 +22,9 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU) endif() elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel) add_compile_options(-warn declarations,general,usage,interfaces,unused) - add_compile_options(-standard-semantics) + if(NOT CMAKE_Fortran_COMPILER_VERSION VERSION_EQUAL 20.2.1.20200827) + add_compile_options(-standard-semantics) + endif() if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0) add_compile_options(-stand f15) else() @@ -35,7 +37,7 @@ endif() # --- compiler feature checks include(CheckFortranSourceCompiles) include(CheckFortranSourceRuns) -check_fortran_source_compiles("error stop i; end" f18errorstop SRC_EXT f90) +check_fortran_source_runs("i=0; error stop i; end" f18errorstop SRC_EXT f90) check_fortran_source_compiles("real, allocatable :: array(:, :, :, :, :, :, :, :, :, :); end" f03rank SRC_EXT f90) check_fortran_source_runs("use, intrinsic :: iso_fortran_env, only : real128; real(real128) :: x; x = x+1; end" f03real128) From 8e95e05633fbb8fc4da16f602250c764851c9cde Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 19:11:19 -0500 Subject: [PATCH 03/39] Add files via upload From 3d8e3e69ef4bce57b210e853362722c5058c12ce Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 19:25:20 -0500 Subject: [PATCH 04/39] Add files via upload --- src/CMakeLists.txt | 3 +++ src/Makefile.manual | 26 +++++++++++--------------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ea7403663..02604959e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,9 @@ # Create a list of the files to be preprocessed set(fppFiles + stdlib_bitsets.fypp + stdlib_bitsets_64.fypp + stdlib_bitsets_large.fypp stdlib_io.fypp stdlib_linalg.fypp stdlib_linalg_diag.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index 129dbe553..872f704c0 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -1,5 +1,8 @@ SRC = f18estop.f90 \ stdlib_ascii.f90 \ + stdlib_bitsets.f90 \ + stdlib_bitsets_64.f90 \ + stdlib_bitsets_large.f90 \ stdlib_error.f90 \ stdlib_io.f90 \ stdlib_kinds.f90 \ @@ -12,10 +15,7 @@ SRC = f18estop.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ stdlib_stats_moment.f90 \ - stdlib_stats_var.f90 \ - stdlib_stats_distribution.f90 \ - stdlib_stats_distribution_rvs.f90 \ - stdlib_stats_distribution_implementation.f90 \ + stdlib_stats_var.f90 LIB = libstdlib.a @@ -43,6 +43,9 @@ clean: # Fortran module dependencies f18estop.o: stdlib_error.o +stdlib_bitsets.o: stdlib_kinds.o +stdlib_bitsets_64.o: stdlib_bitsets.o +stdlib_bitsets_large.o: stdlib_bitsets.o stdlib_error.o: stdlib_optval.o stdlib_io.o: \ stdlib_error.o \ @@ -64,14 +67,11 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o -stdlib_stats_distribution_rvs.o: stdlib_kinds.o -stdlib_stats_distribution.o: \ - stdlib_error.o \ - stdlib_kinds.o \ - stdlib_stats_distribution_rvs.o \ -stdlib_stats_distribution_imp.o: \ - stdlib_stats_distribution.o + # Fortran sources that are built from fypp templates +stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp +stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp +stdlib_bitsets.f90: stdlib_bitsets.fypp stdlib_io.f90: stdlib_io.fypp stdlib_linalg.f90: stdlib_linalg.fypp stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp @@ -80,7 +80,3 @@ stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_rvs.f90: stdlib_stats_distribution_rvs.fypp -stdlib_stats_distribution.f90: stdlib_stats_distribution.fypp -stdlib_stats_distribution_implementation.f90: \ - stdlib_stats_distribution_implementation.fypp \ No newline at end of file From 7d4fc50a7013158c05ee7af5f20f4ee26d0dccf1 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:10:27 -0500 Subject: [PATCH 05/39] initial commit --- src/stdlib_stats_distribution_PRNG.fypp | 212 +++++++ src/stdlib_stats_distribution_beta.fypp | 241 ++++++++ src/stdlib_stats_distribution_gamma.fypp | 323 +++++++++++ src/stdlib_stats_distribution_normal.fypp | 366 ++++++++++++ src/stdlib_stats_distribution_special.fypp | 614 +++++++++++++++++++++ src/stdlib_stats_distribution_uniform.fypp | 482 ++++++++++++++++ 6 files changed, 2238 insertions(+) create mode 100644 src/stdlib_stats_distribution_PRNG.fypp create mode 100644 src/stdlib_stats_distribution_beta.fypp create mode 100644 src/stdlib_stats_distribution_gamma.fypp create mode 100644 src/stdlib_stats_distribution_normal.fypp create mode 100644 src/stdlib_stats_distribution_special.fypp create mode 100644 src/stdlib_stats_distribution_uniform.fypp diff --git a/src/stdlib_stats_distribution_PRNG.fypp b/src/stdlib_stats_distribution_PRNG.fypp new file mode 100644 index 000000000..d1bda107a --- /dev/null +++ b/src/stdlib_stats_distribution_PRNG.fypp @@ -0,0 +1,212 @@ +#:include "common.fypp" +module stdlib_stats_distribution_PRNG + use stdlib_kinds, only: int8, int16, int32, int64 + implicit none + private + integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64) + integer(int64), save :: st(4), si = 614872703977525537_int64 + logical, save :: seed_initialized = .false. + + public :: random_seed + public :: dist_rand + public :: jump + public :: long_jump + + + interface dist_rand + !! Version experimental + !! + !! Generation of random integers with different kinds + !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# + !! description)) + #:for k1, t1 in INT_KINDS_TYPES + module procedure dist_rand_${t1[0]}$${k1}$ + #:endfor + end interface dist_rand + + interface random_seed + !! Version experimental + !! + !! Set seed value for random number generator + !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# + !! description)) + !! + #:for k1, t1 in INT_KINDS_TYPES + module procedure random_distribution_seed_${t1[0]}$${k1}$ + #:endfor + end interface random_seed + + + contains + + #:for k1, t1 in INT_KINDS_TYPES + function dist_rand_${t1[0]}$${k1}$(n) result(res) + !! Random integer generation for various kinds + !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind + !! Result will be operated by bitwise operators to generate desired integer + !! and real pseudorandom numbers + !! + ${t1}$, intent(in) :: n + ${t1}$ :: res + integer :: k + + k = MAX_INT_BIT_SIZE - bit_size(n) + res = shiftr(xoshiro256ss( ), k) + end function dist_rand_${t1[0]}$${k1}$ + + #:endfor + + function xoshiro256ss( ) result (res) + ! Generate random 64-bit integers + ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! http://prng.di.unimi.it/xoshiro256starstar.c + ! + ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid + ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is + ! large enough for any parallel application, and it passes all tests we + ! are aware of. + ! + ! The state must be seeded so that it is not everywhere zero. If you have + ! a 64-bit seed, we suggest to seed a splitmix64 generator and use its + ! output to fill st. + ! + ! Fortran 90 version translated from C by Jim-215-Fisher + ! + integer(int64) :: res, t + + if(.not. seed_initialized) call random_distribution_seed_iint64(si,t) + res = rol64(st(2) * 5 , 7) * 9 + t = shiftl(st(2), 17) + st(3) = ieor(st(3), st(1)) + st(4) = ieor(st(4), st(2)) + st(2) = ieor(st(2), st(3)) + st(1) = ieor(st(1), st(4)) + st(3) = ieor(st(3), t) + st(4) = rol64(st(4), 45) + end function xoshiro256ss + + function rol64(x, k) result(res) + integer(int64), intent(in) :: x + integer, intent(in) :: k + integer(int64) :: t1, t2, res + + t1 = shiftr(x, (64 - k)) + t2 = shiftl(x, k) + res = ior(t1, t2) + end function rol64 + + + subroutine jump + ! This is the jump function for the xoshiro256ss generator. It is equivalent + ! to 2^128 calls to xoshiro256ss(); it can be used to generate 2^128 + ! non-overlapping subsequences for parallel computations. + ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! http://prng.di.unimi.it/xoshiro256starstar.c + ! + ! Fortran 90 version translated from C by Jim-215-Fisher + integer(int64) :: jp(4) = [1733541517147835066_int64, & + -3051731464161248980_int64, & + -6244198995065845334_int64, & + 4155657270789760540_int64] + integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 + integer :: i, j, k + + do i = 1, 4 + do j = 1, 64 + if(iand(jp(i), shiftl(c, j - 1)) /= 0) then + s1 = ieor(s1, st(1)) + s2 = ieor(s2, st(2)) + s3 = ieor(s3, st(3)) + s4 = ieor(s4, st(4)) + end if + k = xoshiro256ss( ) + end do + end do + st(1) = s1 + st(2) = s2 + st(3) = s3 + st(4) = s4 + end subroutine jump + + subroutine long_jump + ! This is the long-jump function for the xoshiro256ss generator. It is + ! equivalent to 2^192 calls to xoshiro256ss(); it can be used to generate + ! 2^64 starting points, from each of which jump() will generate 2^64 + ! non-overlapping subsequences for parallel distributed computations + ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! http://prng.di.unimi.it/xoshiro256starstar.c + ! + ! Fortran 90 version translated from C by Jim-215-Fisher + integer(int64) :: jp(4) = [8566230491382795199_int64, & + -4251311993797857357_int64, & + 8606660816089834049_int64, & + 4111957640723818037_int64] + integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 + integer(int32) :: i, j, k + + do i = 1, 4 + do j = 1, 64 + if(iand(jp(i), shiftl(c, j - 1)) /= 0) then + s1 = ieor(s1, st(1)) + s2 = ieor(s2, st(2)) + s3 = ieor(s3, st(3)) + s4 = ieor(s4, st(4)) + end if + k = xoshiro256ss() + end do + end do + st(1) = s1 + st(2) = s2 + st(3) = s3 + st(4) = s4 + end subroutine long_jump + + function splitmix64(s) result(res) + ! Written in 2015 by Sebastiano Vigna (vigna@acm.org) + ! This is a fixed-increment version of Java 8's SplittableRandom + ! generator. + ! See http://dx.doi.org/10.1145/2714064.2660195 and + ! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html + ! + ! It is a very fast generator passing BigCrush, and it can be useful if + ! for some reason you absolutely want 64 bits of state. + ! + ! Fortran 90 translated from C by Jim-215-Fisher + ! + integer(int64) :: res, int01, int02, int03 + integer(int64), intent(in), optional :: s + data int01, int02, int03/-7046029254386353131_int64, & + -4658895280553007687_int64, & + -7723592293110705685_int64/ + + if(present(s)) si = s + res = si + si = res + int01 + res = ieor(res, shiftr(res, 30)) * int02 + res = ieor(res, shiftr(res, 27)) * int03 + res = ieor(res, shiftr(res, 31)) + end function splitmix64 + + #:for k1, t1 in INT_KINDS_TYPES + subroutine random_distribution_seed_${t1[0]}$${k1}$(put, get) + !! Set seed value for random number generator + !! + ${t1}$, intent(in) :: put + ${t1}$, intent(out) :: get + integer(int64) :: tmp + integer :: i + + tmp = splitmix64(int(put, kind = int64)) + do i = 1, 10 + tmp = splitmix64( ) + end do + do i = 1, 4 + tmp = splitmix64( ) + st(i) = tmp + end do + get = int(tmp, kind = ${k1}$) + seed_initialized = .true. + end subroutine random_distribution_seed_${t1[0]}$${k1}$ + + #:endfor +end module stdlib_stats_distribution_PRNG \ No newline at end of file diff --git a/src/stdlib_stats_distribution_beta.fypp b/src/stdlib_stats_distribution_beta.fypp new file mode 100644 index 000000000..af1707660 --- /dev/null +++ b/src/stdlib_stats_distribution_beta.fypp @@ -0,0 +1,241 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +Module stdlib_stats_distribution_beta + use stdlib_kinds + use stdlib_error, only : error_stop + use stdlib_stats_distribution_uniform, only : uni=>uniform_distribution_rvs + use stdlib_stats_distribution_gamma, only : rgamma=>gamma_distribution_rvs + use stdlib_stats_distribution_special, only : beta, inbeta + + + implicit none + private + integer(int64), parameter :: INT_ONE = 1_int64 + + public :: beta_distribution_rvs + public :: beta_distribution_pdf + public :: beta_distribution_cdf + + interface beta_distribution_rvs + !! Version experimental + !! + !! Beta Distribution Random Variates + !! ([Specification](../page/specs/stdlib_stats_distribution_beta.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_rvs_${t1[0]}$${k1}$ ! 2 arguments + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_rvs_array_${t1[0]}$${k1}$ ! 3 arguments + #:endfor + end interface beta_distribution_rvs + + interface beta_distribution_pdf + !! Version experimental + !! + !! Beta Distribution Probability Density Function + !! ([Specification](../page/specs/stdlib_stats_distribution_beta.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_pdf_${t1[0]}$${k1}$ + #:endfor + end interface beta_distribution_pdf + + interface beta_distribution_cdf + !! Version experimental + !! + !! Beta Distribution Cumulative Distribution Function + !! ([Specification](../page/specs/stdlib_stats_distribution_beta.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_cdf_${t1[0]}$${k1}$ + #:endfor + end interface beta_distribution_cdf + + + contains + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function beta_dist_rvs_${t1[0]}$${k1}$(a, b) result(res) + ! Beta random variate + ! + ${t1}$, intent(in) :: a, b + ${t1}$ :: res, x, y, xx(2) + ${t1}$, parameter :: z = 0.0_${k1}$, one = 1.0_${k1}$ + + if(a <= z .or. b <= z) call error_stop("Error: Beta distribution" & + //" paramters a, b must be greater than zero") + if( a < one .or. b < one) then + do + xx = uni(z, one, 2) + x = xx(1) ** (one / a) + y = xx(2) ** (one / b) + y = x + y + if(y <= one .and. y /= z) exit + end do + else + do + x = rgamma(a); y = rgamma(b) + y = x + y + if( y /= z) exit + end do + endif + res = x / y + return + end function beta_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function beta_dist_rvs_${t1[0]}$${k1}$(a, b) result(res) + ! Beta distributed complex. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: a, b + ${t1}$ :: res + real(${k1}$), parameter :: z = 0.0_${k1}$ + real(${k1}$) :: tr, ti + + tr = beta_dist_rvs_r${k1}$(real(a), real(b)) + ti = beta_dist_rvs_r${k1}$(aimag(a), aimag(b)) + res = cmplx(tr, ti, kind = ${k1}$) + return + end function beta_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function beta_dist_rvs_array_${t1[0]}$${k1}$(a, b, array_size) result(res) + ${t1}$, intent(in) :: a, b + ${t1}$, allocatable :: res(:) + integer, intent(in) :: array_size + integer :: i + ${t1}$ :: x, y, xx(2) + real(${k1}$), parameter :: z = 0.0_${k1}$, one = 1.0_${k1}$ + + allocate(res(array_size)) + if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error: Beta" & + //" distribution paramters a, b must be greater than zero") + if( a < one .or. b < one) then + do i = 1, array_size + do + xx = uni(z, one, 2) + x = xx(1) ** (one / a) + y = xx(2) ** (one / b) + y = x + y + if(y <= one .and. y /= z) exit + end do + res(i) = x / y + end do + else + do i = 1, array_size + do + x = rgamma(a); y = rgamma(b) + y = x + y + if( y /= z) exit + end do + res(i) = x / y + end do + endif + return + end function beta_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + function beta_dist_rvs_array_${t1[0]}$${k1}$(a, b, array_size) result(res) + ${t1}$, intent(in) :: a, b + ${t1}$, allocatable :: res(:) + integer, intent(in) :: array_size + integer :: i + real(${k1}$), allocatable :: tr(:), ti(:) + + allocate(res(array_size), tr(array_size), ti(array_size)) + tr = beta_dist_rvs_array_r${k1}$(real(a), real(b), array_size) + ti = beta_dist_rvs_array_r${k1}$(aimag(a), aimag(b), array_size) + res = cmplx(tr, ti, kind = ${k1}$) + return + end function beta_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function beta_dist_pdf_${t1[0]}$${k1}$(x, a, b) result(res) + ! Beta distributed probability function + ! + ${t1}$, intent(in) :: x, a, b + real :: res + ${t1}$ :: z = 0.0_${k1}$, one = 1.0_${k1}$ + + if(a <= z .or. b <= z) call error_stop("Error: Beta distribution" & + //" parameters a, b must be greater than zero") + if(x == z) then + if(a <= one) then + res = huge(1.0) + else + res = z + endif + elseif(x == one) then + if(b <= one) then + res = huge(1.0) + else + res = z + endif + else + res = x ** (a - 1) * (1 - x) ** (b - 1) / beta(a, b) + endif + return + end function beta_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function beta_dist_pdf_${t1[0]}$${k1}$(x, a, b) result(res) + ${t1}$, intent(in) :: x, a, b + real :: res + + res = beta_dist_pdf_r${k1}$(real(x), real(a), real(b)) + res = res * beta_dist_pdf_r${k1}$(aimag(x), aimag(a), aimag(b)) + return + end function beta_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function beta_dist_cdf_${t1[0]}$${k1}$(x, a, b) result(res) + ! Beta cumulative distribution function + ! + ${t1}$, intent(in) :: x, a, b + real :: res + ${t1}$ :: z = 0.0_${k1}$, one = 1.0_${k1}$ + + if(a <= z .or. b <= z) call error_stop("Error: Beta distribution" & + //" parameters a, b must be greater than zero") + if(x == z) then + res = 0.0 + elseif(x == one) then + res = 1.0 + else + res = inbeta(x, a, b) + endif + return + end function beta_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function beta_dist_cdf_${t1[0]}$${k1}$(x, a, b) result(res) + ${t1}$, intent(in) :: x, a, b + real :: res + + res = beta_dist_cdf_r${k1}$(real(x), real(a), real(b)) + res = res * beta_dist_cdf_r${k1}$(aimag(x), aimag(a), aimag(b)) + end function beta_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + +end module stdlib_stats_distribution_beta \ No newline at end of file diff --git a/src/stdlib_stats_distribution_gamma.fypp b/src/stdlib_stats_distribution_gamma.fypp new file mode 100644 index 000000000..c25d47804 --- /dev/null +++ b/src/stdlib_stats_distribution_gamma.fypp @@ -0,0 +1,323 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +Module stdlib_stats_distribution_gamma + use stdlib_kinds + use stdlib_error, only : error_stop + use stdlib_stats_distribution_uniform, only : uni=>uniform_distribution_rvs + use stdlib_stats_distribution_normal, only : rnor=>normal_distribution_rvs + use stdlib_stats_distribution_special, only : ingamma=>ingamma_low, loggamma + + implicit none + private + integer(int64), parameter :: INT_ONE = 1_int64 + + public :: gamma_distribution_rvs + public :: gamma_distribution_pdf + public :: gamma_distribution_cdf + + interface gamma_distribution_rvs + !! Version experimental + !! + !! Gamma Distribution Random Variates + !! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure gamma_dist_rvs_1_${t1[0]}$${k1}$ ! 1 argument + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure gamma_dist_rvs_${t1[0]}$${k1}$ ! 2 arguments + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure gamma_dist_rvs_array_${t1[0]}$${k1}$ ! 3 arguments + #:endfor + end interface gamma_distribution_rvs + + interface gamma_distribution_pdf + !! Version experimental + !! + !! Gamma Distribution Probability Density Function + !! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure gamma_dist_pdf_${t1[0]}$${k1}$ + #:endfor + end interface gamma_distribution_pdf + + interface gamma_distribution_cdf + !! Version experimental + !! + !! Gamma Distribution Cumulative Distribution Function + !! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure gamma_dist_cdf_${t1[0]}$${k1}$ + #:endfor + end interface gamma_distribution_cdf + + + contains + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) result(res) + ! Gamma random variate + ! + ${t1}$, intent(in) :: shape + ${t1}$ :: res + ${t1}$ :: x, v, u, zz, tol = 1000 * epsilon(1.0_${k1}$) + ${t1}$, save :: alpha = 0._${k1}$, d, c, sq = 0.0331_${k1}$ + + if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" shape parameter must be greater than zero") + zz = shape + if(zz < 1._${k1}$) zz = 1._${k1}$ + zz + if(abs(zz - alpha) > tol) then + alpha = zz + d = alpha - 1._${k1}$ / 3._${k1}$ + c = 1._${k1}$ / (3._${k1}$ * sqrt(d)) + endif + do + do + x = rnor(0.0_${k1}$, 1.0_${k1}$) + v = 1._${k1}$ + c * x + v = v * v * v + if(v > 0._${k1}$) exit + end do + x = x * x + u = uni(1.0_${k1}$) + if(u < (1._${k1}$ - sq * x * x)) exit + if(log(u) < 0.5_${k1}$ * x + d * (1._${k1}$ - v + log(v))) exit + end do + res = d * v + if(shape < 1._${k1}$) then + u = uni(1.0_${k1}$) + res = res * u ** (1._${k1}$ / shape) + endif + return + end function gamma_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) result(res) + ! Gamma distributed complex. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: shape + ${t1}$ :: res + real(${k1}$) :: tr, ti + + tr = gamma_dist_rvs_1_r${k1}$(real(shape)) + ti = gamma_dist_rvs_1_r${k1}$(aimag(shape)) + res = cmplx(tr,ti, kind=${k1}$) + return + end function gamma_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate) & + result(res) + ${t1}$, intent(in) :: shape, rate + ${t1}$ :: res + ${t1}$ :: x, v, u, zz, tol = 1000 * epsilon(1.0_${k1}$) + ${t1}$, save :: alpha = 0._${k1}$, d, c, sq = 0.0331_${k1}$ + + if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" shape parameter must be greater than zero") + if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" rate parameter must be greater than zero") + zz = shape + if(zz < 1._${k1}$) zz = 1._${k1}$ + zz + if(abs(zz - alpha) > tol) then + alpha = zz + d = alpha - 1._${k1}$ / 3._${k1}$ + c = 1._${k1}$ / (3._${k1}$ * sqrt(d)) + endif + do + do + x = rnor(0.0_${k1}$, 1.0_${k1}$) + v = 1._${k1}$ + c * x + v = v * v * v + if(v > 0._${k1}$) exit + end do + x = x * x + u = uni(1.0_${k1}$) + if(u < (1._${k1}$ - sq * x * x)) exit + if(log(u) < 0.5_${k1}$ * x + d * (1._${k1}$ - v + log(v))) exit + end do + res = d * v + if(shape < 1._${k1}$) then + u = uni(1.0_${k1}$) + res = res * u ** (1._${k1}$ / shape) + endif + res = res / rate + return + end function gamma_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate) & + result(res) + ! Gamma distributed complex. The real part and imaginary part are & + ! independent of each other. + ! + ${t1}$, intent(in) :: shape, rate + ${t1}$ :: res + real(${k1}$) :: tr, ti + + tr = gamma_dist_rvs_r${k1}$(real(shape), real(rate)) + ti = gamma_dist_rvs_r${k1}$(aimag(shape), aimag(rate)) + res = cmplx(tr, ti, kind=${k1}$) + return + end function gamma_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function gamma_dist_rvs_array_${t1[0]}$${k1}$(shape, rate, array_size) & + result(res) + ${t1}$, intent(in) :: shape, rate + ${t1}$, allocatable :: res(:) + integer, intent(in) :: array_size + ${t1}$ :: x, v, u, zz, tol = 1000 * epsilon(1.0_${k1}$), re + ${t1}$, save :: alpha = 0._${k1}$, d, c, sq = 0.0331_${k1}$ + integer :: i + + if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" shape parameter must be greater than zero") + if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" rate parameter must be greater than zero") + allocate(res(array_size)) + zz = shape + if(zz < 1._${k1}$) zz = 1._${k1}$ + zz + if(abs(zz - alpha) > tol) then + alpha = zz + d = alpha - 1._${k1}$ / 3._${k1}$ + c = 1._${k1}$ / (3._${k1}$ * sqrt(d)) + endif + do i = 1, array_size + do + do + x = rnor(0.0_${k1}$, 1.0_${k1}$) + v = 1._${k1}$ + c * x + v = v * v * v + if(v > 0._${k1}$) exit + end do + x = x * x + u = uni(1.0_${k1}$) + if(u < (1._${k1}$ - sq * x * x)) exit + if(log(u) < 0.5_${k1}$ * x + d * (1._${k1}$ - v + log(v))) exit + end do + re = d * v + if(shape < 1._${k1}$) then + u = uni(1.0_${k1}$) + re = re * u ** (1._${k1}$ / shape) + endif + res(i) = re / rate + end do + return + end function gamma_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + function gamma_dist_rvs_array_${t1[0]}$${k1}$(shape, rate, array_size) & + result(res) + ${t1}$, intent(in) :: shape, rate + ${t1}$, allocatable :: res(:) + integer, intent(in) :: array_size + integer :: i + real(${k1}$) :: tr, ti + + allocate(res(array_size)) + do i = 1, array_size + tr = gamma_dist_rvs_r${k1}$(real(shape), real(rate)) + ti = gamma_dist_rvs_r${k1}$(aimag(shape), aimag(rate)) + res(i) = cmplx(tr, ti, kind=${k1}$) + end do + return + end function gamma_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function gamma_dist_pdf_${t1[0]}$${k1}$(x, shape, rate) & + result(res) + ! Gamma distributed probability function + ! + ${t1}$, intent(in) :: x, shape, rate + real :: res + + if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" rate parameter must be greaeter than zero") + if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" shape parameter must be greater than zero") + if(x <= 0.0_${k1}$) call error_stop("Error: Gamma distribution variate" & + //" x must be greater than zero") + if(x == 0.0_${k1}$) then + if(shape <= 1.0_${k1}$) then + res = huge(1.0) + 1.0 + else + res = 0.0_${k1}$ + endif + else + res = exp((shape - 1._${k1}$) * log(x) - x * rate + shape * & + log(rate) - loggamma(shape)) + endif + return + end function gamma_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function gamma_dist_pdf_${t1[0]}$${k1}$(x, shape, rate) & + result(res) + ${t1}$, intent(in) :: x, shape, rate + real :: res + + res = gamma_dist_pdf_r${k1}$(real(x), real(shape), real(rate)) + res = res * gamma_dist_pdf_r${k1}$(aimag(x), aimag(shape), aimag(rate)) + return + end function gamma_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate) & + result(res) + ! Gamma random cumulative distribution function + ! + ${t1}$, intent(in) :: x, shape, rate + real :: res + + if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" rate parameter must be greaeter than zero") + if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & + //" shape parameter must be greater than zero") + if(x <= 0.0_${k1}$) call error_stop("Error: Gamma distribution variate" & + //" x must be greater than zero") + res = ingamma(shape, rate * x) / gamma(shape) + return + end function gamma_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate) & + result(res) + ${t1}$, intent(in) :: x, shape, rate + real :: res + + res = gamma_dist_cdf_r${k1}$(real(x), real(shape), real(rate)) + res = res * gamma_dist_cdf_r${k1}$(aimag(x), aimag(shape), aimag(rate)) + end function gamma_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + +end module stdlib_stats_distribution_gamma \ No newline at end of file diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp new file mode 100644 index 000000000..8397ba08c --- /dev/null +++ b/src/stdlib_stats_distribution_normal.fypp @@ -0,0 +1,366 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +Module stdlib_stats_distribution_normal + use stdlib_kinds + use stdlib_error, only : error_stop + use stdlib_stats_distribution_PRNG, only : dist_rand + use stdlib_stats_distribution_uniform, only : uni=>uniform_distribution_rvs + + implicit none + private + + real(dp), parameter :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp + integer, save :: kn(0:127) + real(dp), save :: wn(0:127), fn(0:127) + logical, save :: zig_norm_initialized = .false. + + public :: normal_distribution_rvs + public :: normal_distribution_pdf + public :: normal_distribution_cdf + + interface normal_distribution_rvs + !! Version experimental + !! + !! Normal Distribution Random Variates + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + module procedure norm_dist_rvs_0_rsp !0 dummy variable + + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_rvs_${t1[0]}$${k1}$ !2 dummy variables + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_rvs_array_${t1[0]}$${k1}$ !3 dummy variables + #:endfor + end interface normal_distribution_rvs + + interface normal_distribution_pdf + !! Version experimental + !! + !! Normal Distribution Probability Density Function + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_pdf_${t1[0]}$${k1}$ + #:endfor + end interface normal_distribution_pdf + + interface normal_distribution_cdf + !! Version experimental + !! + !! Normal Distribution Cumulative Distribution Function + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_cdf_${t1[0]}$${k1}$ + #:endfor + end interface normal_distribution_cdf + + + 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 + ! + ! N.B. It is assumed that all integers are 32-bit. + ! + ! Latest version - 1 January 2001 + ! + real(dp), parameter :: M1 = 2147483648.0_dp + real(dp) :: dn = 3.442619855899_dp, tn, & + vn = 0.00991256303526217_dp, q + integer :: i + + 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. + return + end subroutine zigset + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_rvs_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 + ! original algorithm use 32bit + hz = dist_rand(1_int32) + + iz = iand( hz, 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 + + !original algorithm use 32bit + 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 + return + end function norm_dist_rvs_0_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Normal random variate (loc, scale) + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y + integer :: hz, iz + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + if( .not. zig_norm_initialized ) call zigset + iz = 0 + ! original algorithm use 32bit + hz = dist_rand(1_int32) + + iz = iand( hz, 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 + + !original algorithm use 32bit + 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 + res = res * scale + loc + return + end function norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Normal 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 = norm_dist_rvs_r${k1}$(real(loc), real(scale)) + ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + res = cmplx(tr, ti, kind=${k1}$) + return + end function norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + ${t1}$, allocatable :: res(:) + ${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: Normal distribution scale" & + //" parameter must be non-zero") + if( .not. zig_norm_initialized ) call zigset + allocate(res(array_size)) + do i = 1, array_size + iz = 0 + ! original algorithm use 32bit + 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 + + !original algorithm use 32bit + 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 + return + end function norm_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + integer :: i + ${t1}$, allocatable :: res(:) + real(${k1}$) :: tr, ti + + allocate(res(array_size)) + do i = 1, array_size + tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) + ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + res(i) = cmplx(tr, ti, kind=${k1}$) + end do + return + end function norm_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ! Normal distributed probability function + ! + ${t1}$, intent(in) :: x, loc, scale + real :: res + ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & + (sqrt_2_Pi * scale) + return + end function norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + res = norm_dist_pdf_r${k1}$(real(x), real(loc), real(scale)) + res = res * norm_dist_pdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) + return + end function norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ! Normal random cumulative distribution function + ! + ${t1}$, intent(in) :: x, loc, scale + real :: res + ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ + return + end function norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + res = norm_dist_cdf_r${k1}$(real(x), real(loc), real(scale)) + res = res * norm_dist_cdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) + return + end function norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + +end module stdlib_stats_distribution_normal \ No newline at end of file diff --git a/src/stdlib_stats_distribution_special.fypp b/src/stdlib_stats_distribution_special.fypp new file mode 100644 index 000000000..ca1992f9e --- /dev/null +++ b/src/stdlib_stats_distribution_special.fypp @@ -0,0 +1,614 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES +Module stdlib_stats_distribution_special + use stdlib_kinds + use stdlib_error, only : error_stop + + implicit none + private + real(qp), parameter :: D(0:10) = [2.48574089138753565546e-5_qp, & + 1.05142378581721974210_qp, & + -3.45687097222016235469_qp, & + 4.51227709466894823700_qp, & + -2.98285225323576655721_qp, & + 1.05639711577126713077_qp, & + -1.95428773191645869583e-1_qp, & + 1.70970543404441224307e-2_qp, & + -5.71926117404305781283e-4_qp, & + 4.63399473359905636708e-6_qp, & + -2.71994908488607703910e-9_qp] + real(qp), parameter :: R = 10.900511_qp, HALF = 0.5_qp, & + sqep = log(2.0_qp * sqrt(exp(1.0_qp) / acos(-1.0_qp))) + real(dp), parameter :: ep_machine = 2.2e-16_dp, dm = 1.0e-300_dp + + ! for stdlib_distribution internal use + + public :: loggamma, log_factorial + public :: ingamma_low, log_ingamma_low, ingamma_up, log_ingamma_up + public :: regamma_p, regamma_q + public :: beta, log_beta, inbeta + + interface loggamma + ! Logrithm of gamma function with real variable + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure l_gamma_${t1[0]}$${k1}$ + #:endfor + end interface loggamma + + interface log_factorial + ! Logrithm of factorial n!, integer variable + ! + #:for k1, t1 in INT_KINDS_TYPES + module procedure l_factorial_1_${t1[0]}$${k1}$ !1 dummy + #:endfor + + #: for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure l_factorial_${t1[0]}$${k1}$${k2}$ !2 dummy + #:endfor + #:endfor + end interface log_factorial + + + interface ingamma_low + ! Lower incomplete gamma function + ! + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + module procedure ingamma_low_${t1[0]}$${k1}$${k2}$ + #:endif + #:endfor + #:endfor + end interface ingamma_low + + interface log_ingamma_low + ! Logrithm of lower incomplete gamma function + ! + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$ + #:endif + #:endfor + #:endfor + end interface log_ingamma_low + + interface ingamma_up + ! Upper incomplete gamma function + ! + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + module procedure ingamma_up_${t1[0]}$${k1}$${k2}$ + #:endif + #:endfor + #:endfor + end interface ingamma_up + + interface log_ingamma_up + ! Logrithm of upper incomplete gamma function + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$ + #:endif + #:endfor + #:endfor + end interface log_ingamma_up + + interface regamma_p + ! Regularized (normalized) lower incomplete gamma function, P + ! + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + module procedure regamma_p_${t1[0]}$${k1}$${k2}$ + #:endif + #:endfor + #:endfor + end interface regamma_p + + interface regamma_q + ! Regularized (normalized) upper incomplete gamma function, Q + ! + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + module procedure regamma_q_${t1[0]}$${k1}$${k2}$ + #:endif + #:endfor + #:endfor + end interface regamma_q + + interface gpx + ! Evaluation of incomplete gamma function + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure gpx_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure gpx_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + end interface gpx + + interface beta + ! Beta function + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure beta_${t1[0]}$${k1}$ + #:endfor + end interface beta + + interface log_beta + ! Logrithm of beta function + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure l_beta_${t1[0]}$${k1}$ + #:endfor + end interface log_beta + + interface inbeta + ! Incomplete beta function + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure inbeta_${t1[0]}$${k1}$ + #:endfor + end interface inbeta + + + contains + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function l_gamma_${t1[0]}$${k1}$(x) result (res) + ! Log gamma function for any positive real number i,e, {R+} + ! + ${t1}$, intent(in) :: x + ${t1}$ :: res + real(qp) :: q, sum + integer :: i + + if(x <= 0._${k1}$) call error_stop("Error: Gamma function augument" & + //" must be greater than 0") + if(x == 1.0_${k1}$ .or. x == 2.0_${k1}$) then + res = 0.0_${k1}$ + else + q = real(x, qp) - HALF + sum = D(0) + do i=1, 10 + sum = sum + D(i) / (real(x, qp) - 1.0_qp + real(i, qp)) + end do + res = real(sqep + log(sum) - q + q * log(q + R), kind=${k1}$) + endif + return + end function l_gamma_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function l_factorial_1_${t1[0]}$${k1}$(n) result(res) + ! Log(n!) with single precision result, n is integer + ! + ${t1}$, intent(in) :: n + real :: res + + if(n < 0) call error_stop("Error: Factorial function augument must" & + //" be no less than 0") + select case(n) + case (0) + res = 0.0 + case (1) + res = 0.0 + case (2:) + res = loggamma(real(n+1, dp)) + end select + return + end function l_factorial_1_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function l_factorial_${t1[0]}$${k1}$${k2}$(n,x) result(res) + ! Log(n!) with required prescision for result, n is integer, x is a real & + ! for specified kind + ! + ${t1}$, intent(in) :: n + ${t2}$, intent(in) :: x + ${t2}$ :: res + + if(n < 0) call error_stop("Error: factorial function augument must" & + //" be no less than 0") + select case(n) + case (0) + res = 0.0_${k2}$ + case (1) + res = 0.0_${k2}$ + case (2:) + res = loggamma(real(n + 1, kind=${k2}$)) + end select + return + end function l_factorial_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function gpx_${t1[0]}$${k1}$(s, x) result(res) + ! Approximation of incomplete gamma G function + ! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and + ! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM + ! Transactions on Mathematical Software, March 2020. + ! + ! Fortran 90 program by Jim-215-Fisher + ! + ${t1}$, intent(in) :: x, s + real(dp) :: res + real(dp) :: a, b, g, c, d, y + integer :: n + + if(x < 0._${k1}$) then + call error_stop("Error: Incomplete gamma function with negative x" & + //" must come with integer of s") + elseif(s >= x) then + a = real(s, dp) + g = 1.0_dp / a + c = g + do + a = a + 1.0_dp + c = c * real(x, dp) / a + g = g + c + if(abs(c) < ep_machine) exit + end do + else + a = 1.0_dp + b = real(x + 1 - s, dp) + g = a / b + c = a / dm + d = 1.0_dp / b + n = 2 + do + a = -(n - 1) * real((n - 1 - s), dp) + b = real(x - s, dp) + 2 * n - 1.0_dp + d = d * a + b + if(d == 0.0_dp) d = dm + c = b + a / c + if(c == 0.0_dp) c = dm + d = 1.0_dp / d + y = c * d + g = g * y + n = n + 1 + if(abs(y - 1.0_dp) < ep_machine) exit + end do + endif + res = g + return + end function gpx_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function gpx_${t1[0]}$${k1}$${k2}$(s, x) result(res) + ! Approximation of incomplete gamma G function + ! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and + ! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM + ! Transactions on Mathematical Software, March 2020. + ! + ${t1}$, intent(in) :: s + ${t2}$, intent(in) :: x + real(dp) :: res + ${t2}$ :: p_lim + real(dp) :: a, b, g, c, d, y + integer :: n + + if(x < -9._${k2}$) then + p_lim = 5.0_${k2}$ * (sqrt(abs(x)) - 1.0_${k2}$) + elseif(x >= -9.0_${k2}$ .and. x <= 0.0_${k2}$) then + p_lim = 0.0_${k2}$ + else + p_lim = x + endif + if(real(s, ${k2}$) >= p_lim) then + a = real(s, dp) + g = 1.0_dp / a + c = g + do + a = a + 1.0_dp + c = c * x / a + g = g + c + if(abs(c) < ep_machine) exit + end do + elseif(x >= 0.0_${k2}$) then + a = 1.0_dp + b = real(x, dp) + (1 - s) + g = a / b + c = a / dm + d = 1.0_dp / b + n = 2 + do + a = -(n - 1) * real((n - s - 1), dp) + b = real(x - s, dp) + 2 * n - 1.0_dp + d = d * a + b + if(d == 0.0_dp) d = dm + c = b + a / c + if(c == 0.0_dp) c = dm + d = 1.0_dp / d + y = c * d + g = g * y + n = n + 1 + if(abs(y - 1.0_dp) < ep_machine) exit + end do + elseif(abs(x) > real(max(1_${k1}$, s - 1_${k1}$), ${k2}$)) then + a = real(-x, dp) + c = 1.0_dp / a + d = real(s - 1, dp) + b = c * (a - d) + n = 1 + do + c = d * (d - 1.0_dp) / (a * a) + d = d - 2.0_dp + y = c * ( a - d) + b = b + y + n = n + 1 + if(int(n, ${k1}$) > (s - 2_${k1}$) / 2_${k1}$ .or. y < b * & + ep_machine) exit + end do + if(y >= b * ep_machine .and. mod(s, 2_${k1}$) /= 0_${k1}$) & + b = b + d * c / a + g = ((-1) ** s * exp(-a + loggamma(real(s, dp)) - (s - 1) * & + log(a)) + b ) / a + endif + res = g + return + end function gpx_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(s, x) & + result(res) + ! Approximation of lower incomplete gamma function + ! + ${t1}$, intent(in) :: s + ${t2}$, intent(in) :: x + ${t2}$ :: res + real(dp) :: s1, y, xx, ss + + #:if t1[0] == "i" + if(s < 0_${k1}$) call error_stop("Error: Lower incomplete gamma" & + //" function input s value must be greater than 0") + #:else + if(s < 0._${k1}$) call error_stop("Error: Lower incomplete gamma" & + //" function input s value must be greater than 0") + #:endif + xx = real(x, dp); ss = real(s, dp) + if(x == 0.0_${k2}$) then + res = 0.0_${k2}$ + elseif(x > 0.0_${k2}$ .and. x <= real(s, ${k2}$)) then + s1 = -xx + ss * log(xx) + res = real(gpx(s,x) * exp(s1), kind=${k2}$) + elseif(x > real(s, ${k2}$)) then + s1 = loggamma(ss) + y = 1.0_dp - exp(-xx + ss * log(xx) - s1) * gpx(s,x) + res = real(y * exp(s1), kind=${k2}$) + else + s1 = -xx + ss * log(-xx) + res = real((-1)**s * gpx(s,x) * exp(s1), kind=${k2}$) + endif + return + end function ingamma_low_${t1[0]}$${k1}$${k2}$ + #:endif + #:endfor + #:endfor + + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(s, x) & + result(res) + ${t1}$, intent(in) :: s + ${t2}$, intent(in) :: x + ${t2}$ :: res + + res = log(ingamma_low(s,x)) + end function l_ingamma_low_${t1[0]}$${k1}$${k2}$ + + #:endif + #:endfor + #:endfor + + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(s, x) & + result(res) + ! Approximation of upper incomplete gamma function + ! + ${t1}$, intent(in) :: s + ${t2}$, intent(in) :: x + ${t2}$ :: res + + res = exp(loggamma(real(s, kind=${k2}$))) - ingamma_low(s,x) + return + end function ingamma_up_${t1[0]}$${k1}$${k2}$ + + #:endif + #:endfor + #:endfor + + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and ( k1 != k2)) + impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(s, x) & + result(res) + ${t1}$, intent(in) :: s + ${t2}$, intent(in) :: x + ${t2}$ :: res + + res = log(ingamma_up(s,x)) + end function l_ingamma_up_${t1[0]}$${k1}$${k2}$ + + #:endif + #:endfor + #:endfor + + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + impure elemental function regamma_p_${t1[0]}$${k1}$${k2}$(s, x) result(res) + ! Approximation of regulated incomplet gamma function P(s,x) + ! + ${t1}$, intent(in) :: s + ${t2}$, intent(in) :: x + ${t2}$ :: res + real(dp) :: s1, xx, ss + + #:if t1[0] == "i" + if(s < 0_${k1}$) call error_stop("Error: Lower incomplete gamma" & + //" function input s value must be greater than 0") + #:else + if(s < 0._${k1}$) call error_stop("Error: Lower incomplete gamma" & + //" function input s value must be greater than 0") + #:endif + xx = real(x, dp); ss = real(s, dp) + s1 = -xx + ss * log(abs(xx)) - loggamma(ss) + if(x == 0.0_${k2}$) then + res = 0.0_${k2}$ + elseif(x > 0.0_${k2}$ .and. x <= real(s, ${k2}$)) then + res = real(gpx(s,x) * exp(s1), kind=${k2}$) + elseif(x > real(s, ${k2}$)) then + res = 1.0_dp - exp(s1) * gpx(s,x) + else + res = real((-1)**s * gpx(s,x) * exp(s1), kind=${k2}$) + endif + return + end function regamma_p_${t1[0]}$${k1}$${k2}$ + + #:endif + #:endfor + #:endfor + + #:for k1, t1 in IR_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + #:if not ((t1[0] == "r") and (k1 != k2)) + impure elemental function regamma_q_${t1[0]}$${k1}$${k2}$(s, x) & + result(res) + ! Approximation of regulated incomplet gamma function Q(s,x) + ! + ${t1}$, intent(in) :: s + ${t2}$, intent(in) :: x + ${t2}$ :: res + + res = real(1.0_dp - regamma_p(s,x), kind=${k2}$) + return + end function regamma_q_${t1[0]}$${k1}$${k2}$ + + #:endif + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function beta_${t1[0]}$${k1}$(a, b) result(res) + ! Evaluation of beta function through gamma function + ! + ${t1}$, intent(in) :: a, b + ${t1}$ :: res + + if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error: Beta" & + //" function auguments a, b values must be greater than 0") + res = exp(loggamma(a) + loggamma(b) - loggamma(a+b)) + return + end function beta_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function l_beta_${t1[0]}$${k1}$(a, b) result(res) + ! Logrithm of beta function through log(gamma) + ! + ${t1}$, intent(in) :: a, b + ${t1}$ :: res + + if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error: Beta" & + //" function auguments a, b values must be greater than 0") + res = loggamma(a) + loggamma(b) - loggamma(a+b) + return + end function l_beta_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function inbeta_${t1[0]}$${k1}$(x, a, b) result(res) + ! Evaluation of incomplete beta function using continued fractions + ! "Computation of Special Functions" by S. Zhang and J. Jin, 1996 + ! + ${t1}$, intent(in) :: x, a, b + ${t1}$ :: res + integer :: n, k + real(dp) :: an, bn, g, c, d, y, s0, ak, ak2 + + if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error:" & + //" Incomplete beta function auguments a, b values must be" & + //" greater than 0") + s0 = (a + 1) / (a + b + 2) + an = 1.0_dp + bn = 1.0_dp + g = an / bn + c = an / dm + d = 1.0_dp / bn + n = 1 + if(x < real(s0, ${k1}$)) then + do + if(mod(n, 2) == 0) then + k = n / 2; ak = real(a + 2 * k, dp) + an = k * real(x, dp) * (b - k) / (ak * ak - ak) + else + k = (n - 1) / 2; ak = real(a + k, dp); ak2 = ak + k + an = - (ak + b) * ak * real(x, dp) / (ak2 * ak2 + ak2) + endif + d = d * an + bn + if(d == 0.0_dp) d = dm + c = bn + an / c + if(c == 0.0_dp) c = dm + d = 1.0_dp / d + y = c * d + g = g * y + n = n + 1 + if(abs(y - 1.0_dp) < ep_machine) exit + end do + g = x ** a * (1.0_${k1}$ - x) ** b * g / (a * beta(a, b)) + else + do + if(mod(n, 2) == 0) then + k = n / 2; ak = real(b + 2 * k, dp) + an = k * (1.0_dp - x) * (a - k) / (ak * ak - ak) + else + k = (n - 1) / 2; ak = b + k; ak2 = ak + k + an = - ak * (1.0_dp - x) * (a + ak) / (ak2 * ak2 + ak2) + endif + d = d * an + bn + if(d == 0.0_dp) d = dm + c = bn + an / c + if(c == 0.0_dp) c = dm + d = 1.0_dp / d + y = c * d + g = g * y + n = n + 1 + if(abs(y - 1.0_dp) < ep_machine) exit + end do + g = x ** a * (1.0_${k1}$ - x) ** b * g / (b * beta(a, b)) + g = 1.0_${k1}$ - g + endif + res = g + end function inbeta_${t1[0]}$${k1}$ + + #:endfor + +end module stdlib_stats_distribution_special diff --git a/src/stdlib_stats_distribution_uniform.fypp b/src/stdlib_stats_distribution_uniform.fypp new file mode 100644 index 000000000..572caaf37 --- /dev/null +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -0,0 +1,482 @@ +#: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_stats_distribution_PRNG, 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 :: uniform_distribution_rvs + public :: uniform_distribution_pdf + public :: uniform_distribution_cdf + public :: shuffle + + interface uniform_distribution_rvs + !! Version experimental + !! + !! Get uniformly distributed random variate for integer, real and complex + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! description)) + + module procedure unif_dist_rvs_0_rsp ! 0 dummy variable + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_rvs_1_${t1[0]}$${k1}$ ! 1 dummy variable + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_rvs_${t1[0]}$${k1}$ ! 2 dummy variables + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_rvs_array_${t1[0]}$${k1}$ ! 3 dummy variables + #:endfor + end interface uniform_distribution_rvs + + interface uniform_distribution_pdf + !! Version experiment + !! + !! Get uniform distribution probability density (pdf) for integer, real and complex variables + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! description)) + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_pdf_${t1[0]}$${k1}$ + #:endfor + end interface uniform_distribution_pdf + + interface uniform_distribution_cdf + !! Version experimental + !! + !! Get uniform distribution cumulative distribution function (cdf) for integer, real and complex variables + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! description)) + !! + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_cdf_${t1[0]}$${k1}$ + #:endfor + end interface uniform_distribution_cdf + + 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# + !! description)) + !! + #: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 unif_dist_rvs_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, n + integer :: zeros, bits_left, bits + + n = scale + if(n <= 0_${k1}$) call error_stop("Error: Uniform distribution scale" & + //" parameter must be positive") + zeros = leadz(n) + bits = bit_size(n) - zeros + mask = shiftr(not(0_${k1}$), zeros) + L1 : do + u = dist_rand(n) + res = iand(u, mask) + if(res <= n) exit L1 + bits_left = zeros + L2 : do + if(bits_left < bits) exit L2 + u = shiftr(u, bits) + res = iand(u, mask) + if(res <= n) exit L1 + bits_left = bits_left - bits + end do L2 + end do L1 + return + end function unif_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function unif_dist_rvs_${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: Uniform distribution" & + //" scale parameter must be non-zero") + res = loc + unif_dist_rvs_1_${t1[0]}$${k1}$(scale) + return + end function unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function unif_dist_rvs_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] + return + end function unif_dist_rvs_0_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function unif_dist_rvs_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: Uniform distribution" & + //" scale parameter must be non-zero") + res = scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) + return + end function unif_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function unif_dist_rvs_${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: Uniform distribution" & + //" scale parameter must be non-zero") + res = loc + scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) + return + end function unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function unif_dist_rvs_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, r2, tr, ti + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & + //" distribution scale parameter must be non-zero") + r1 = unif_dist_rvs_0_r${k1}$( ) + if(real(scale) == 0.0_${k1}$) then + ti = aimag(scale) * r1 + tr = 0.0_${k1}$ + elseif(aimag(scale) == 0.0_${k1}$) then + tr = real(scale) * r1 + ti = 0.0_${k1}$ + else + r2 = unif_dist_rvs_0_r${k1}$( ) + tr = real(scale) * r1 + ti = aimag(scale) * r2 + endif + res = cmplx(tr, ti, kind=${k1}$) + return + end function unif_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function unif_dist_rvs_${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, r2, tr, ti + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & + //" distribution scale parameter must be non-zero") + r1 = unif_dist_rvs_0_r${k1}$( ) + if(real(scale) == 0.0_${k1}$) then + tr = real(loc) + ti = aimag(loc) + aimag(scale) * r1 + elseif(aimag(scale) == 0.0_${k1}$) then + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + else + r2 = unif_dist_rvs_0_r${k1}$( ) + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + aimag(scale) * r2 + endif + res = cmplx(tr, ti, kind=${k1}$) + return + end function unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + function unif_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + ${t1}$, allocatable :: res(:) + ${t1}$ :: u, mask, n, nn + integer, intent(in) :: array_size + integer :: i, zeros, bits_left, bits + + n = scale + if(n == 0_${k1}$) call error_stop("Error: Uniform distribution" & + //" scale parameter must be non-zero") + allocate(res(array_size)) + zeros = leadz(n) + bits = bit_size(n) - zeros + mask = shiftr(not(0_${k1}$), zeros) + do i = 1, array_size + L1 : do + u = dist_rand(n) + nn = iand(u, mask) + if(nn <= n) exit L1 + bits_left = zeros + L2 : do + if(bits_left < bits) exit L2 + u = shiftr(u, bits) + nn = iand(u, mask) + if(nn <= n) exit L1 + bits_left = bits_left - bits + end do L2 + end do L1 + res(i) = loc + nn + end do + return + end function unif_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function unif_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + ${t1}$, allocatable :: res(:) + ${t1}$ :: t + integer, intent(in) :: array_size + integer(int64) :: tmp + integer :: i + + + if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & + //" scale parameter must be non-zero") + allocate(res(array_size)) + do i = 1, array_size + tmp = shiftr(dist_rand(INT_ONE), 11) + t = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + res(i) = loc + scale * t + enddo + return + end function unif_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + function unif_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + ${t1}$, allocatable :: res(:) + real(${k1}$) :: r1, r2, tr, ti + integer, intent(in) :: array_size + integer(int64) :: tmp + integer :: i + + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & + //" distribution scale parameter must be non-zero") + allocate(res(array_size)) + do i = 1, array_size + tmp = shiftr(dist_rand(INT_ONE), 11) + r1 = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + if(real(scale) == 0.0_${k1}$) then + tr = real(loc) + ti = aimag(loc) + aimag(scale) * r1 + elseif(aimag(scale) == 0.0_${k1}$) then + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + else + tmp = shiftr(dist_rand(INT_ONE), 11) + r2 = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + aimag(scale) * r2 + endif + res(i) = cmplx(tr, ti, kind=${k1}$) + enddo + return + end function unif_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + elemental function unif_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0) then + res = 0.0 + elseif(x < loc .or. x > (loc + scale)) then + res = 0.0 + else + res = 1. / (scale + 1_${k1}$) + end if + return + end function unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function unif_dist_pdf_${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 + elseif(x <= loc .or. x >= (loc + scale)) then + res = 0.0 + else + res = 1.0 / scale + end if + return + end function unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function unif_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + real(${k1}$) :: tr, ti + + tr = real(loc) + real(scale); ti = aimag(loc) + aimag(scale) + if(scale == (0.0_${k1}$,0.0_${k1}$)) then + res = 0.0 + elseif((real(x) >= real(loc) .and. real(x) <= tr) .and. & + (aimag(x) >= aimag(loc) .and. aimag(x) <= ti)) then + res = 1.0 / (real(scale) * aimag(scale)) + else + res = 0.0 + end if + return + end function unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + elemental function unif_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0) then + res = 0.0 + elseif(x < loc) then + res = 0.0 + elseif(x >= loc .and. x <= (loc + scale)) then + res = real((x - loc + 1_${k1}$)) / real((scale + 1_${k1}$)) + else + res = 1.0 + end if + return + end function unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function unif_dist_cdf_${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 + elseif(x < loc) then + res = 0.0 + elseif(x >= loc .and. x <= (loc + scale)) then + res = (x - loc) / scale + else + res = 1.0 + end if + return + end function unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function unif_dist_cdf_${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 + endif + r1 = real(x) < real(loc) + r2 = real(x) > (real(loc) + real(scale)) + i1 = aimag(x) < aimag(loc) + i2 = aimag(x) > (aimag(loc) + aimag(scale)) + if(r1 .or. i1) then + res = 0.0 + elseif((.not. r1) .and. (.not. r2) .and. i2) then + res = (real(x) - real(loc)) / real(scale) + elseif((.not. i1) .and. (.not. i2) .and. r2) then + res = (aimag(x) - aimag(loc)) / aimag(scale) + elseif((.not. r1) .and. (.not. r2) .and. (.not. i1) .and. (.not. i2)) & + then + res = (real(x) - real(loc)) * (aimag(x) - aimag(loc)) / & + (real(scale) * aimag(scale)) + elseif(r2 .and. i2)then + res = 1.0 + end if + return + end function unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + function shuffle_${t1[0]}$${k1}$( list ) result(res) + ${t1}$, intent(in) :: list(:) + ${t1}$, allocatable :: res(:) + ${t1}$ :: tmp + integer :: n, i, j + + n = size(list) + allocate(res(n), source=list) + do i = 1, n - 1 + j = uniform_distribution_rvs(n - i) + i + tmp = res(i) + res(i) = res(j) + res(j) = tmp + end do + return + end function shuffle_${t1[0]}$${k1}$ + + #:endfor +end module stdlib_stats_distribution_uniform \ No newline at end of file From 8ef351698afb9d93fa32cfcc5bdba6b648113a44 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:11:52 -0500 Subject: [PATCH 06/39] initial commit --- src/tests/stats/test_distribution_beta.f90 | 562 +++++++++++++++++++++ 1 file changed, 562 insertions(+) create mode 100644 src/tests/stats/test_distribution_beta.f90 diff --git a/src/tests/stats/test_distribution_beta.f90 b/src/tests/stats/test_distribution_beta.f90 new file mode 100644 index 000000000..f83f69bd4 --- /dev/null +++ b/src/tests/stats/test_distribution_beta.f90 @@ -0,0 +1,562 @@ +program test_distribution_beta + use stdlib_kinds + use stdlib_error, only : check + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_beta, beta_rvs => beta_distribution_rvs, & + beta_pdf => beta_distribution_pdf, & + beta_cdf => beta_distribution_cdf + + implicit none + real(sp), parameter :: sptol = 1000 * epsilon(1.0_sp) + real(dp), parameter :: dptol = 1000 * epsilon(1.0_dp) + real(qp), parameter :: qptol = 1000 * epsilon(1.0_qp) + logical :: warn = .true. + integer :: put, get + + real :: aa(2,3,4), bb(2,3,4), x(2,3,4) + complex :: a, b + + put = 1234567 + call random_seed(put, get) + print *, beta_cdf(0.4, 0.5,1.0) + ! a standard beta cumulative at 0.4 with a=0.5, b=1.0 + +! 0.842700839 + + print *, beta_cdf(0.8, 1.5,2.0) + ! a cumulative at 0.8 with a a=1.5, b=2.0 + +! 0.953988254 + + aa(:,:,:) = 2.0 + bb(:,:,:) = 3.0 + x = reshape(beta_rvs(2.0, 3.0, 24),[2,3,4]) + !beta random variates array with a a=2.0, b=3.0 + print *, beta_cdf(x,aa,bb) ! a rank 3 standard beta cumulative array + +! [0.710880339, 0.472411335, 0.578345954, 0.383050948, 0.870905757, +! 0.870430350, 0.170215249, 0.677347481, 0.620089889, 0.161825046, +! 4.17549349E-02, 0.510665894, 0.252201647, 0.911497891, 0.984424412, +! 0.635621786, 0.177783430, 0.414842933, 0.871342421, 0.338317066, +! 2.06879266E-02, 0.335232288, 0.907408893, 0.624871135] + + a = (.7, 2.1) + b = (0.5,1.0) + print *, beta_cdf((0.5,0.5),a,b) +stop + put = 1234567 + call random_seed(put, get) + + call test_beta_random_generator + + call test_beta_rvs_rsp + call test_beta_rvs_rdp + call test_beta_rvs_rqp + call test_beta_rvs_csp + call test_beta_rvs_cdp + call test_beta_rvs_cqp + + call test_beta_pdf_rsp + call test_beta_pdf_rdp + call test_beta_pdf_rqp + call test_beta_pdf_csp + call test_beta_pdf_cdp + call test_beta_pdf_cqp + + call test_beta_cdf_rsp + call test_beta_cdf_rdp + call test_beta_cdf_rqp + call test_beta_cdf_csp + call test_beta_cdf_cdp + call test_beta_cdf_cqp + stop + + contains + + subroutine test_beta_random_generator + integer :: i, j, freq(0:999), num=1000000 + real(dp) :: chisq, expct + + print *, "" + print *, "Test beta random generator with chi-squared" + freq = 0 + do i = 1, num + j = 1000 * beta_cdf(beta_rvs(2.0,1.5),2.0,1.5) + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + expct = num / 1000 + do i = 0, 999 + 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 beta random generator is : ", chisq + call check((chisq < 1143.9), & + msg="beta randomness failed chi-squared test", warn=warn) + end subroutine test_beta_random_generator + + subroutine test_beta_rvs_rsp + real(sp) :: res(10), a, b + integer :: i, n, k = 5 + integer :: put, get + real(sp) :: ans(10) = [0.744399697952416243334363353017662363_sp, & + 0.785582064888561104409969900319307115_sp, & + 0.290228167791285215460609614976244262_sp, & + 0.540957122824810300186112642749043673_sp, & + 0.866783498753591906081043187617603418_sp, & + 0.164859290722944895746841956125886140_sp, & + 0.752018475270934089892015814754187410_sp, & + 0.535463312713066631219371237676531884_sp, & + 0.438125081488452966935618841567308566_sp, & + 0.635255468090749026953184924665020348_sp] + + print *, "Test beta_distribution_rvs_rsp" + put = 639741825 + call random_seed(put, get) + a = 2.0_sp; b = 1.0_sp + do i = 1, 5 + res(i) = beta_rvs(a, b) + end do + res(6:10) = beta_rvs(a, b, k) + call check(all(abs(res - ans) < sptol), & + msg="beta_distribution_rvs_rsp failed", warn=warn) + end subroutine test_beta_rvs_rsp + + subroutine test_beta_rvs_rdp + real(dp) :: res(10), a, b + integer :: i, n, k = 5 + integer :: put, get + real(dp) :: ans(10) = [0.744399697952416243334363353017662363_dp, & + 0.785582064888561104409969900319307115_dp, & + 0.290228167791285215460609614976244262_dp, & + 0.540957122824810300186112642749043673_dp, & + 0.866783498753591906081043187617603418_dp, & + 0.164859290722944895746841956125886140_dp, & + 0.752018475270934089892015814754187410_dp, & + 0.535463312713066631219371237676531884_dp, & + 0.438125081488452966935618841567308566_dp, & + 0.635255468090749026953184924665020348_dp] + + print *, "Test beta_distribution_rvs_rdp" + put = 639741825 + call random_seed(put, get) + a = 2.0_dp; b = 1.0_dp + do i = 1, 5 + res(i) = beta_rvs(a, b) + end do + res(6:10) = beta_rvs(a, b, k) + call check(all(abs(res - ans) < dptol), & + msg="beta_distribution_rvs_rdp failed", warn=warn) + end subroutine test_beta_rvs_rdp + + subroutine test_beta_rvs_rqp + real(qp) :: res(10), a, b + integer :: i, n, k = 5 + integer :: put, get + real(qp) :: ans(10) = [0.744399697952416243334363353017662363_qp, & + 0.785582064888561104409969900319307115_qp, & + 0.290228167791285215460609614976244262_qp, & + 0.540957122824810300186112642749043673_qp, & + 0.866783498753591906081043187617603418_qp, & + 0.164859290722944895746841956125886140_qp, & + 0.752018475270934089892015814754187410_qp, & + 0.535463312713066631219371237676531884_qp, & + 0.438125081488452966935618841567308566_qp, & + 0.635255468090749026953184924665020348_qp] + + print *, "Test beta_distribution_rvs_rqp" + put = 639741825 + call random_seed(put, get) + a = 2.0_qp; b = 1.0_qp + do i = 1, 5 + res(i) = beta_rvs(a, b) + end do + res(6:10) = beta_rvs(a, b, k) + call check(all(abs(res - ans) < qptol), & + msg="beta_distribution_rvs_rqp failed", warn=warn) + end subroutine test_beta_rvs_rqp + + subroutine test_beta_rvs_csp + complex(sp) :: res(10), a, b + integer :: i, n, k = 5 + integer :: put, get + complex(sp) :: ans(10) = [(0.996593894945797558163416198727906880_sp, & + 0.132598000233262366166275960014792683_sp), & + (0.390100279193128267998817589162392551_sp, & + 0.594960539319605102054215589048597662_sp), & + (0.839296137442072004654073615963339625_sp, & + 0.811403350202500212373591232224849046_sp), & + (0.915863048173886740665234455471312086_sp, & + 0.791004162831226427993329892479385152_sp), & + (0.449544461366363638651136427537418482_sp, & + 6.648931970435189824360854746280901949E-0002_sp), & + (0.418599563122841869588123385796545717_sp, & + 2.682053757872834248817007196452821345E-0002_sp), & + (0.847048136644577210744803917816801455_sp, & + 0.728350780933130458744978339193529096_sp), & + (0.859055679164195250196327014741975834_sp, & + 0.677632443230547158536307254984227248_sp), & + (0.251814018668272502730175719600969280_sp, & + 6.410086007672458486925924399711432844E-0003_sp), & + (0.656763996944608477801486502669830627_sp, & + 0.870919913077248989108181221896598931_sp)] + + print *, "Test beta_distribution_rvs_csp" + put = 639741825 + call random_seed(put, get) + a = (2.0_sp, 0.7_sp); b = (0.8_sp, 1.2_sp) + do i = 1, 5 + res(i) = beta_rvs(a, b) + end do + res(6:10) = beta_rvs(a, b, k) + call check(all(abs(res - ans) < sptol), & + msg="beta_distribution_rvs_csp failed", warn=warn) + end subroutine test_beta_rvs_csp + + subroutine test_beta_rvs_cdp + complex(dp) :: res(10), a, b + integer :: i, n, k = 5 + integer :: put, get + complex(dp) :: ans(10) = [(0.996593894945797558163416198727906880_dp, & + 0.132598000233262366166275960014792683_dp), & + (0.390100279193128267998817589162392551_dp, & + 0.594960539319605102054215589048597662_dp), & + (0.839296137442072004654073615963339625_dp, & + 0.811403350202500212373591232224849046_dp), & + (0.915863048173886740665234455471312086_dp, & + 0.791004162831226427993329892479385152_dp), & + (0.449544461366363638651136427537418482_dp, & + 6.648931970435189824360854746280901949E-0002_dp), & + (0.418599563122841869588123385796545717_dp, & + 2.682053757872834248817007196452821345E-0002_dp), & + (0.847048136644577210744803917816801455_dp, & + 0.728350780933130458744978339193529096_dp), & + (0.859055679164195250196327014741975834_dp, & + 0.677632443230547158536307254984227248_dp), & + (0.251814018668272502730175719600969280_dp, & + 6.410086007672458486925924399711432844E-0003_dp), & + (0.656763996944608477801486502669830627_dp, & + 0.870919913077248989108181221896598931_dp)] + + print *, "Test beta_distribution_rvs_cdp" + put = 639741825 + call random_seed(put, get) + a = (2.0_dp, 0.7_dp); b = (0.8_dp, 1.2_dp) + do i = 1, 5 + res(i) = beta_rvs(a, b) + end do + res(6:10) = beta_rvs(a, b, k) + call check(all(abs(res - ans) < dptol), & + msg="beta_distribution_rvs_cdp failed", warn=warn) + end subroutine test_beta_rvs_cdp + + subroutine test_beta_rvs_cqp + complex(qp) :: res(10), a, b + integer :: i, n, k = 5 + integer :: put, get + complex(qp) :: ans(10) = [(0.996593894945797558163416198727906880_qp, & + 0.132598000233262366166275960014792683_qp), & + (0.390100279193128267998817589162392551_qp, & + 0.594960539319605102054215589048597662_qp), & + (0.839296137442072004654073615963339625_qp, & + 0.811403350202500212373591232224849046_qp), & + (0.915863048173886740665234455471312086_qp, & + 0.791004162831226427993329892479385152_qp), & + (0.449544461366363638651136427537418482_qp, & + 6.648931970435189824360854746280901949E-0002_qp), & + (0.418599563122841869588123385796545717_qp, & + 2.682053757872834248817007196452821345E-0002_qp), & + (0.847048136644577210744803917816801455_qp, & + 0.728350780933130458744978339193529096_qp), & + (0.859055679164195250196327014741975834_qp, & + 0.677632443230547158536307254984227248_qp), & + (0.251814018668272502730175719600969280_qp, & + 6.410086007672458486925924399711432844E-0003_qp), & + (0.656763996944608477801486502669830627_qp, & + 0.870919913077248989108181221896598931_qp)] + + print *, "Test beta_distribution_rvs_cqp" + put = 639741825 + call random_seed(put, get) + a = (2.0_qp, 0.7_qp); b = (0.8_qp, 1.2_qp) + do i = 1, 5 + res(i) = beta_rvs(a, b) + end do + res(6:10) = beta_rvs(a, b, k) + call check(all(abs(res - ans) < qptol), & + msg="beta_distribution_rvs_cqp failed", warn=warn) + end subroutine test_beta_rvs_cqp + + + + subroutine test_beta_pdf_rsp + real(sp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [1.97584832, 1.97584832, 1.97584832, 1.94792151, & + 1.05610514, 1.63085556, 1.44469929, 1.78598392, & + 1.03530371, 1.32163048, 1.95935822, 1.49064910, & + 1.22708333, 0.816426575, 1.93443334] + + print *, "Test beta_distribution_pdf_rsp" + put = 345987126 + call random_seed(put, get) + a = 2.0_sp; b = 1.0_sp + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_pdf(x1, a, b) + res(:, 2:5) = beta_pdf(x2, a, b) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="beta_distribution_pdf_rsp failed", warn=warn) + end subroutine test_beta_pdf_rsp + + subroutine test_beta_pdf_rdp + real(dp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [1.97584832, 1.97584832, 1.97584832, 1.94792151, & + 1.05610514, 1.63085556, 1.44469929, 1.78598392, & + 1.03530371, 1.32163048, 1.95935822, 1.49064910, & + 1.22708333, 0.816426575, 1.93443334] + + print *, "Test beta_distribution_pdf_rdp" + put = 345987126 + call random_seed(put, get) + a = 2.0_dp; b = 1.0_dp + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_pdf(x1, a, b) + res(:, 2:5) = beta_pdf(x2, a, b) + call check(all(abs(res - reshape(ans, [3,5])) < dptol), & + msg="beta_distribution_pdf_rdp failed", warn=warn) + end subroutine test_beta_pdf_rdp + + subroutine test_beta_pdf_rqp + real(qp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [1.97584832, 1.97584832, 1.97584832, 1.94792151, & + 1.05610514, 1.63085556, 1.44469929, 1.78598392, & + 1.03530371, 1.32163048, 1.95935822, 1.49064910, & + 1.22708333, 0.816426575, 1.93443334] + + print *, "Test beta_distribution_pdf_rqp" + put = 345987126 + call random_seed(put, get) + a = 2.0_qp; b = 1.0_qp + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_pdf(x1, a, b) + res(:, 2:5) = beta_pdf(x2, a, b) + call check(all(abs(res - reshape(ans, [3,5])) < qptol), & + msg="beta_distribution_pdf_rqp failed", warn=warn) + end subroutine test_beta_pdf_rqp + + subroutine test_beta_pdf_csp + complex(sp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [2.79032898, 2.79032898, 2.79032898, 0.831092536, & + 0.859560609, 0.833086491, 2.18498087, 0.611756265, & + 5.17011976, 0.805453956, 2.12247658, 0.969030082, & + 1.92922175, 0.700230777, 6.4548239] + + print *, "Test beta_distribution_pdf_csp" + put = 345987126 + call random_seed(put, get) + a = (2.0_sp, 0.7_sp); b = (0.8_sp, 1.2_sp) + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_pdf(x1, a, b) + res(:, 2:5) = beta_pdf(x2, a, b) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="beta_distribution_pdf_csp failed", warn=warn) + end subroutine test_beta_pdf_csp + + subroutine test_beta_pdf_cdp + complex(dp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [2.79032898, 2.79032898, 2.79032898, 0.831092536, & + 0.859560609, 0.833086491, 2.18498087, 0.611756265, & + 5.17011976, 0.805453956, 2.12247658, 0.969030082, & + 1.92922175, 0.700230777, 6.4548239] + + print *, "Test beta_distribution_pdf_cdp" + put = 345987126 + call random_seed(put, get) + a = (2.0_dp, 0.7_dp); b = (0.8_dp, 1.2_dp) + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_pdf(x1, a, b) + res(:, 2:5) = beta_pdf(x2, a, b) + call check(all(abs(res - reshape(ans, [3,5])) < dptol), & + msg="beta_distribution_pdf_cdp failed", warn=warn) + end subroutine test_beta_pdf_cdp + + subroutine test_beta_pdf_cqp + complex(qp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [2.79032898, 2.79032898, 2.79032898, 0.831092536, & + 0.859560609, 0.833086491, 2.18498087, 0.611756265, & + 5.17011976, 0.805453956, 2.12247658, 0.969030082, & + 1.92922175, 0.700230777, 6.4548239] + + print *, "Test beta_distribution_pdf_cqp" + put = 345987126 + call random_seed(put, get) + a = (2.0_qp, 0.7_qp); b = (0.8_qp, 1.2_qp) + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_pdf(x1, a, b) + res(:, 2:5) = beta_pdf(x2, a, b) + call check(all(abs(res - reshape(ans, [3,5])) < qptol), & + msg="beta_distribution_pdf_cqp failed", warn=warn) + end subroutine test_beta_pdf_cqp + + + subroutine test_beta_cdf_rsp + real(sp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [0.153344765, 0.153344765, 0.153344765, 0.639326215, & + 0.227737889, 0.832331538, 0.215463713, 0.609950244, & + 0.552298367, 0.936580479, 0.473157555, 0.375768840, & + 2.33022049E-02, 0.907276988, 0.230596066] + + print *, "Test beta_distribution_cdf_rsp" + put = 567985123 + call random_seed(put, get) + a = 2.0_sp; b = 2.0_sp + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_cdf(x1, a, b) + res(:, 2:5) = beta_cdf(x2, a, b) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg="beta_distribution_cdf_rsp failed", warn=warn) + end subroutine test_beta_cdf_rsp + + subroutine test_beta_cdf_rdp + real(dp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [0.153344765, 0.153344765, 0.153344765, 0.639326215, & + 0.227737889, 0.832331538, 0.215463713, 0.609950244, & + 0.552298367, 0.936580479, 0.473157555, 0.375768840, & + 2.33022049E-02, 0.907276988, 0.230596066] + + print *, "Test beta_distribution_cdf_rdp" + put = 567985123 + call random_seed(put, get) + a = 2.0_dp; b = 2.0_dp + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_cdf(x1, a, b) + res(:, 2:5) = beta_cdf(x2, a, b) + call check(all(abs(res - reshape(ans,[3,5])) < dptol), & + msg="beta_distribution_cdf_rdp failed", warn=warn) + end subroutine test_beta_cdf_rdp + + subroutine test_beta_cdf_rqp + real(qp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [0.153344765, 0.153344765, 0.153344765, 0.639326215, & + 0.227737889, 0.832331538, 0.215463713, 0.609950244, & + 0.552298367, 0.936580479, 0.473157555, 0.375768840, & + 2.33022049E-02, 0.907276988, 0.230596066] + + print *, "Test beta_distribution_cdf_rqp" + put = 567985123 + call random_seed(put, get) + a = 2.0_qp; b = 2.0_qp + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_cdf(x1, a, b) + res(:, 2:5) = beta_cdf(x2, a, b) + call check(all(abs(res - reshape(ans,[3,5])) < qptol), & + msg="beta_distribution_cdf_rqp failed", warn=warn) + end subroutine test_beta_cdf_rqp + + subroutine test_beta_cdf_csp + complex(sp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [7.06288144E-02, 7.06288144E-02, 7.06288144E-02, & + 3.14221904E-02, 0.229956210, 4.80622090E-02, & + 0.495213449, 0.541507423, 0.105880715, 0.194856107, & + 3.40392105E-02, 1.09316744E-02, 0.180127904, & + 0.654031873, 0.406583667] + + print *, "Test beta_distribution_cdf_csp" + put = 567985123 + call random_seed(put, get) + a = (2.0_sp, 0.7_sp); b = (0.8_sp, 1.2_sp) + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_cdf(x1, a, b) + res(:, 2:5) = beta_cdf(x2, a, b) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg="beta_distribution_cdf_csp failed", warn=warn) + end subroutine test_beta_cdf_csp + + subroutine test_beta_cdf_cdp + complex(dp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [7.06288144E-02, 7.06288144E-02, 7.06288144E-02, & + 3.14221904E-02, 0.229956210, 4.80622090E-02, & + 0.495213449, 0.541507423, 0.105880715, 0.194856107, & + 3.40392105E-02, 1.09316744E-02, 0.180127904, & + 0.654031873, 0.406583667] + + print *, "Test beta_distribution_cdf_cdp" + put = 567985123 + call random_seed(put, get) + a = (2.0_dp, 0.7_dp); b = (0.8_dp, 1.2_dp) + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_cdf(x1, a, b) + res(:, 2:5) = beta_cdf(x2, a, b) + call check(all(abs(res - reshape(ans,[3,5])) < dptol), & + msg="beta_distribution_cdf_cdp failed", warn=warn) + end subroutine test_beta_cdf_cdp + + subroutine test_beta_cdf_cqp + complex(qp) :: x1, x2(3,4), a, b + real :: res(3,5) + integer :: i, n + integer :: put, get + real :: ans(15) = [7.06288144E-02, 7.06288144E-02, 7.06288144E-02, & + 3.14221904E-02, 0.229956210, 4.80622090E-02, & + 0.495213449, 0.541507423, 0.105880715, 0.194856107, & + 3.40392105E-02, 1.09316744E-02, 0.180127904, & + 0.654031873, 0.406583667] + + print *, "Test beta_distribution_cdf_cqp" + put = 567985123 + call random_seed(put, get) + a = (2.0_qp, 0.7_qp); b = (0.8_qp, 1.2_qp) + x1 = beta_rvs(a, b) + x2 = reshape(beta_rvs(a, b, 12), [3,4]) + res(:,1) = beta_cdf(x1, a, b) + res(:, 2:5) = beta_cdf(x2, a, b) + call check(all(abs(res - reshape(ans,[3,5])) < qptol), & + msg="beta_distribution_cdf_cqp failed", warn=warn) + end subroutine test_beta_cdf_cqp + + +end program test_distribution_beta \ No newline at end of file From b7af23ce5aa7d9dfe5a94b48402c9a5c8e9fdd1d Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:12:27 -0500 Subject: [PATCH 07/39] Update CMakeLists.txt --- src/tests/stats/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 36ffc7aeb..89ce5d221 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -5,6 +5,7 @@ ADDTEST(moment) ADDTEST(rawmoment) ADDTEST(var) ADDTEST(varn) +ADDTEST(distribution_beta) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) From f5d37c887f5d5d5ea878f2e35f10f4309b257737 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:12:49 -0500 Subject: [PATCH 08/39] Update Makefile.manual --- src/tests/stats/Makefile.manual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index aacaf98ab..b8abb662b 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,3 +1,3 @@ -PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 +PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 test_distribution_beta.f90 include ../Makefile.manual.test.mk From d5c2e0b0624af1256daeb45df814fb8b8bb470e0 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:13:34 -0500 Subject: [PATCH 09/39] initial commit --- src/CMakeLists.txt | 75 +++++++--------------------------------- src/Makefile.manual | 83 ++------------------------------------------- 2 files changed, 15 insertions(+), 143 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 02604959e..3fbf70055 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,67 +1,16 @@ -### Pre-process: .fpp -> .f90 via Fypp +ADDTEST(corr) +ADDTEST(cov) +ADDTEST(mean) +ADDTEST(moment) +ADDTEST(rawmoment) +ADDTEST(var) +ADDTEST(varn) +ADDTEST(distribution_normal) -# Create a list of the files to be preprocessed -set(fppFiles - stdlib_bitsets.fypp - stdlib_bitsets_64.fypp - stdlib_bitsets_large.fypp - stdlib_io.fypp - stdlib_linalg.fypp - stdlib_linalg_diag.fypp - stdlib_optval.fypp - stdlib_stats.fypp - stdlib_stats_corr.fypp - stdlib_stats_cov.fypp - stdlib_stats_mean.fypp - stdlib_stats_moment.fypp - stdlib_stats_var.fypp - stdlib_quadrature.fypp - stdlib_quadrature_trapz.fypp - stdlib_quadrature_simps.fypp -) - - -# Custom preprocessor flags if(DEFINED CMAKE_MAXIMUM_RANK) - set(fyppFlags "-DMAXRANK=${CMAKE_MAXIMUM_RANK}") + if(${CMAKE_MAXIMUM_RANK} GREATER 7) + ADDTEST(mean_f03) + endif() elseif(f03rank) - set(fyppFlags) -else() - set(fyppFlags "-DVERSION90") -endif() - -fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) - -set(SRC - stdlib_ascii.f90 - stdlib_error.f90 - stdlib_kinds.f90 - stdlib_logger.f90 - stdlib_system.F90 - ${outFiles} -) - -add_library(fortran_stdlib ${SRC}) - -set(LIB_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/mod_files/) -set_target_properties(fortran_stdlib PROPERTIES - Fortran_MODULE_DIRECTORY ${LIB_MOD_DIR}) -target_include_directories(fortran_stdlib PUBLIC - $ - $ -) - -if(f18errorstop) - target_sources(fortran_stdlib PRIVATE f18estop.f90) -else() - target_sources(fortran_stdlib PRIVATE f08estop.f90) + ADDTEST(mean_f03) endif() - -add_subdirectory(tests) - -install(TARGETS fortran_stdlib - RUNTIME DESTINATION bin - ARCHIVE DESTINATION lib - LIBRARY DESTINATION lib - ) -install(DIRECTORY ${LIB_MOD_DIR} DESTINATION include) diff --git a/src/Makefile.manual b/src/Makefile.manual index 872f704c0..9a82079e0 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -1,82 +1,5 @@ -SRC = f18estop.f90 \ - stdlib_ascii.f90 \ - stdlib_bitsets.f90 \ - stdlib_bitsets_64.f90 \ - stdlib_bitsets_large.f90 \ - stdlib_error.f90 \ - stdlib_io.f90 \ - stdlib_kinds.f90 \ - stdlib_linalg.f90 \ - stdlib_linalg_diag.f90 \ - stdlib_logger.f90 \ - stdlib_optval.f90 \ - stdlib_quadrature.f90 \ - stdlib_quadrature_trapz.f90 \ - stdlib_stats.f90 \ - stdlib_stats_mean.f90 \ - stdlib_stats_moment.f90 \ - stdlib_stats_var.f90 -LIB = libstdlib.a +PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 \ + test_distribution_normal.f90 - - -OBJS = $(SRC:.f90=.o) -MODS = $(OBJS:.o=.mod) -SMODS = $(OBJS:.o=*.smod) - -.PHONY: all clean - -all: $(LIB) - -$(LIB): $(OBJS) - ar rcs $@ $(OBJS) - -clean: - $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) - -%.o: %.f90 - $(FC) $(FFLAGS) -c $< - -%.f90: %.fypp - fypp $(FYPPFLAGS) $< $@ - -# Fortran module dependencies -f18estop.o: stdlib_error.o -stdlib_bitsets.o: stdlib_kinds.o -stdlib_bitsets_64.o: stdlib_bitsets.o -stdlib_bitsets_large.o: stdlib_bitsets.o -stdlib_error.o: stdlib_optval.o -stdlib_io.o: \ - stdlib_error.o \ - stdlib_optval.o \ - stdlib_kinds.o -stdlib_linalg_diag.o: stdlib_kinds.o -stdlib_logger.o: stdlib_ascii.o stdlib_optval.o -stdlib_optval.o: stdlib_kinds.o -stdlib_quadrature.o: stdlib_kinds.o -stdlib_stats_mean.o: \ - stdlib_optval.o \ - stdlib_kinds.o \ - stdlib_stats.o -stdlib_stats_moment.o: \ - stdlib_optval.o \ - stdlib_kinds.o \ - stdlib_stats.o -stdlib_stats_var.o: \ - stdlib_optval.o \ - stdlib_kinds.o \ - stdlib_stats.o - -# Fortran sources that are built from fypp templates -stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp -stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp -stdlib_bitsets.f90: stdlib_bitsets.fypp -stdlib_io.f90: stdlib_io.fypp -stdlib_linalg.f90: stdlib_linalg.fypp -stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp -stdlib_quadrature.f90: stdlib_quadrature.fypp -stdlib_stats.f90: stdlib_stats.fypp -stdlib_stats_mean.f90: stdlib_stats_mean.fypp -stdlib_stats_moment.f90: stdlib_stats_moment.fypp -stdlib_stats_var.f90: stdlib_stats_var.fypp +include ../Makefile.manual.test.mk \ No newline at end of file From 734ab024c88d66daa21d03a2c1de3f17efaff321 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:14:13 -0500 Subject: [PATCH 10/39] initial commit --- src/Makefile.manual | 99 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 96 insertions(+), 3 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 9a82079e0..4617499ba 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -1,5 +1,98 @@ +noSRC = f18estop.f90 \ + stdlib_ascii.f90 \ + stdlib_bitsets.f90 \ + stdlib_bitsets_64.f90 \ + stdlib_bitsets_large.f90 \ + stdlib_error.f90 \ + stdlib_io.f90 \ + stdlib_kinds.f90 \ + stdlib_linalg.f90 \ + stdlib_linalg_diag.f90 \ + stdlib_logger.f90 \ + stdlib_optval.f90 \ + stdlib_quadrature.f90 \ + stdlib_quadrature_trapz.f90 \ + stdlib_stats.f90 \ + stdlib_stats_mean.f90 \ + stdlib_stats_moment.f90 \ + stdlib_stats_var.f90 \ + stdlib_stats_distribution_PRNG.f90 \ + stdlib_stats_distribution_uniform.f90 \ + stdlib_stats_distribution_normal.f90 + +LIB = libstdlib.a -PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 \ - test_distribution_normal.f90 -include ../Makefile.manual.test.mk \ No newline at end of file + +OBJS = $(SRC:.f90=.o) +MODS = $(OBJS:.o=.mod) +SMODS = $(OBJS:.o=*.smod) + +.PHONY: all clean + +all: $(LIB) + +$(LIB): $(OBJS) + ar rcs $@ $(OBJS) + +clean: + $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) + +%.o: %.f90 + $(FC) $(FFLAGS) -c $< + +%.f90: %.fypp + fypp $(FYPPFLAGS) $< $@ + +# Fortran module dependencies +f18estop.o: stdlib_error.o +stdlib_bitsets.o: stdlib_kinds.o +stdlib_bitsets_64.o: stdlib_bitsets.o +stdlib_bitsets_large.o: stdlib_bitsets.o +stdlib_error.o: stdlib_optval.o +stdlib_io.o: \ + stdlib_error.o \ + stdlib_optval.o \ + stdlib_kinds.o +stdlib_linalg_diag.o: stdlib_kinds.o +stdlib_logger.o: stdlib_ascii.o stdlib_optval.o +stdlib_optval.o: stdlib_kinds.o +stdlib_quadrature.o: stdlib_kinds.o +stdlib_stats_mean.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_moment.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_var.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_distribution_PRNG.o: stdlib_kinds.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o +stdlib_stats_distribution_normal.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution.PRNG.o \ + stdlib_stats_distribution.uniform.o + +# Fortran sources that are built from fypp templates +stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp +stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp +stdlib_bitsets.f90: stdlib_bitsets.fypp +stdlib_io.f90: stdlib_io.fypp +stdlib_linalg.f90: stdlib_linalg.fypp +stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp +stdlib_quadrature.f90: stdlib_quadrature.fypp +stdlib_stats.f90: stdlib_stats.fypp +stdlib_stats_mean.f90: stdlib_stats_mean.fypp +stdlib_stats_moment.f90: stdlib_stats_moment.fypp +stdlib_stats_var.f90: stdlib_stats_var.fypp +stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp +stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp +stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp \ No newline at end of file From 67a5093b41a5b9240c697deda6fb9c39af9f2ab5 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:14:38 -0500 Subject: [PATCH 11/39] initial commit --- src/CMakeLists.txt | 75 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 63 insertions(+), 12 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3fbf70055..02604959e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,16 +1,67 @@ -ADDTEST(corr) -ADDTEST(cov) -ADDTEST(mean) -ADDTEST(moment) -ADDTEST(rawmoment) -ADDTEST(var) -ADDTEST(varn) -ADDTEST(distribution_normal) +### Pre-process: .fpp -> .f90 via Fypp +# Create a list of the files to be preprocessed +set(fppFiles + stdlib_bitsets.fypp + stdlib_bitsets_64.fypp + stdlib_bitsets_large.fypp + stdlib_io.fypp + stdlib_linalg.fypp + stdlib_linalg_diag.fypp + stdlib_optval.fypp + stdlib_stats.fypp + stdlib_stats_corr.fypp + stdlib_stats_cov.fypp + stdlib_stats_mean.fypp + stdlib_stats_moment.fypp + stdlib_stats_var.fypp + stdlib_quadrature.fypp + stdlib_quadrature_trapz.fypp + stdlib_quadrature_simps.fypp +) + + +# Custom preprocessor flags if(DEFINED CMAKE_MAXIMUM_RANK) - if(${CMAKE_MAXIMUM_RANK} GREATER 7) - ADDTEST(mean_f03) - endif() + set(fyppFlags "-DMAXRANK=${CMAKE_MAXIMUM_RANK}") elseif(f03rank) - ADDTEST(mean_f03) + set(fyppFlags) +else() + set(fyppFlags "-DVERSION90") +endif() + +fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) + +set(SRC + stdlib_ascii.f90 + stdlib_error.f90 + stdlib_kinds.f90 + stdlib_logger.f90 + stdlib_system.F90 + ${outFiles} +) + +add_library(fortran_stdlib ${SRC}) + +set(LIB_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/mod_files/) +set_target_properties(fortran_stdlib PROPERTIES + Fortran_MODULE_DIRECTORY ${LIB_MOD_DIR}) +target_include_directories(fortran_stdlib PUBLIC + $ + $ +) + +if(f18errorstop) + target_sources(fortran_stdlib PRIVATE f18estop.f90) +else() + target_sources(fortran_stdlib PRIVATE f08estop.f90) endif() + +add_subdirectory(tests) + +install(TARGETS fortran_stdlib + RUNTIME DESTINATION bin + ARCHIVE DESTINATION lib + LIBRARY DESTINATION lib + ) +install(DIRECTORY ${LIB_MOD_DIR} DESTINATION include) From 229ec35ca3ff1343ea19b4004fbac55051d9696c Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:23:34 -0500 Subject: [PATCH 12/39] Update CMakeLists.txt --- src/CMakeLists.txt | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 02604959e..2af329d5a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,6 +18,12 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp + stdlib_stats_distribution_PRNG.fypp + stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp + stdlib_stats_distribution_special.fypp + stdlib_stats_distribution_gamma.fypp + stdlib_stats_distribution_beta.fypp ) From 25dce48e33ee21e31046666f82da914fbff69712 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:29:45 -0500 Subject: [PATCH 13/39] Update Makefile.manual --- src/Makefile.manual | 37 +++++++++++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 4617499ba..59af055b6 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -1,4 +1,4 @@ -noSRC = f18estop.f90 \ +SRC = f18estop.f90 \ stdlib_ascii.f90 \ stdlib_bitsets.f90 \ stdlib_bitsets_64.f90 \ @@ -18,7 +18,10 @@ noSRC = f18estop.f90 \ stdlib_stats_var.f90 \ stdlib_stats_distribution_PRNG.f90 \ stdlib_stats_distribution_uniform.f90 \ - stdlib_stats_distribution_normal.f90 + stdlib_stats_distribution_normal.f90 \ + stdlib_stats_distribution_special.f90 \ + stdlib_stats_distribution_gamma.f90 \ + stdlib_stats_distribution_beta.f90 LIB = libstdlib.a @@ -72,15 +75,30 @@ stdlib_stats_var.o: \ stdlib_stats.o stdlib_stats_distribution_PRNG.o: stdlib_kinds.o stdlib_stats_distribution_uniform.o: \ - stdlib_kinds.o \ + stdlib_kinds.o \ stdlib_error.o \ stdlib_stats_distribution_PRNG.o stdlib_stats_distribution_normal.o: \ - stdlib_kinds.o \ + stdlib_kinds.o \ stdlib_error.o \ - stdlib_stats_distribution.PRNG.o \ - stdlib_stats_distribution.uniform.o - + stdlib_stats_distribution_PRNG.o \ + stdlib_stats_distribution_uniform.o +stdlib_stats_distribution_special.o: \ + stdlib_kinds.o \ + stdlib_error.o +stdlib_stats_distribution_gamma.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_uniform.o \ + stdlib_stats_distribution_normal.o \ + stdlib_stats_distribution_special.o +stdlib_stats_distribution_beta.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_uniform.o \ + stdlib_stats_distribution_gamma.o \ + stdlib_stats_distribution_special.o + # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp @@ -95,4 +113,7 @@ stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp -stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp \ No newline at end of file +stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp +stdlib_stats_distribution_special.f90: stdlib_stats_distribution_special.fypp +stdlib_stats_distribution_gamma.f90: stdlib_stats_distribution_gamma.fypp +stdlib_stats_distribution_beta.f90: stdlib_stats_distribution_beta.fypp From e779c56ad16a0d5bf4083199331dc52cb0c1c3cf Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:31:55 -0500 Subject: [PATCH 14/39] initial commit --- doc/specs/stdlib_stats_distribution_beta.md | 240 ++++++++++++++++++++ 1 file changed, 240 insertions(+) create mode 100644 doc/specs/stdlib_stats_distribution_beta.md diff --git a/doc/specs/stdlib_stats_distribution_beta.md b/doc/specs/stdlib_stats_distribution_beta.md new file mode 100644 index 000000000..3d74c11d3 --- /dev/null +++ b/doc/specs/stdlib_stats_distribution_beta.md @@ -0,0 +1,240 @@ +--- + +title: stats_distribution +--- + +# Statistical Distributions -- Beta Distribution Module + +[TOC] + +## `beta_distribution_rvs` - beta distribution random variates + +### Status + +Experimental + +### Description + +With two auguments for shape parameters a>0, b>0, the function returns a beta distributed random variate Beta(a,b), also known as the beta distribution of the first kind. The function is elemental. For complex auguments, the real and imaginary parts are independent of each other. + +With three auguments, the function return a rank one array of beta distributed random variate Beta(a, b). + +### Syntax + +`result = [[stdlib_stats_distribution_beta(module):beta_distribution_rvs(interface)]](a, b [, array_size])` + +### Arguments + +`a` : has `intent(in)` ans is a scalar of type `real` or `complx`. + +`b`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. + +### Return value + +The result is a scalar or rank one array, with a size of `array_size`, of type `real` or `complx`. + +### Example + +```fortran +program demo_beta_rvs + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_beta, only: rbeta => beta_distribution_rvs + + implicit none + real :: aa(2,3,4), bb(2,3,4) + complex :: a, b + integer :: put, get + + put = 1234567 + call random_seed(put, get) + + print *, rbeta(0.5, 2.0) + !single standard beta random variate with shape a=0.5, b=2.0 + +! 3.38860527E-02 + + print *, rbeta(3.0,2.0) !beta random variate with a=3.0, b=2.0 + +! 0.570277154 + + aa(:,:,:) = 0.8; bb(:,:,:)=0.6 + print *, rbeta(aa, bb) + !a rank 3 array of 24 beta random variates with a=0.8, b=0.6 + +! [0.251766384, 0.578202426, 0.539138556, 0.210130826, 0.908130825, +! 0.880996943, 9.49194748E-03, 0.945992589, 0.290732056, 0.803920329, +! 7.64303207E-02, 0.943150401, 0.927998245, 0.831781328, 0.671169102, +! 0.983966410, 0.289062619, 0.801237404, 0.891931713, 0.897902310, +! 0.845606744, 1.50359496E-02, 0.913162351, 0.915781260] + + print *, rbeta(1.2,0.7,10) + ! an array of 10 random variates with a=1.2, b=0.7 + +! [3.69704030E-02, 0.748639643, 0.896924615, 0.350249082, 0.813054740, +! 0.448072791, 9.39368978E-02, 0.475665420, 0.567661405, 0.893835664] + + a = (3.0, 4.0) + b = (2.0, 0.7) + print *, rbeta(a, b) + !single complex beta random variate with real part of a = 3.0, b=2.0; imagainary part of a=4.0, b=0.7 + +! (0.730923295,0.878909111) + +end program demo_beta_rvs +``` + +## `beta_distribution_pdf` - beta probability density function + +### Status + +Experimental + +### Description + +The probability density function of the continuous beta distribution. + +$$f(x)=\frac{x^{a-1}(1-x)^{b-1}}{B(a,b)}; where\; B(a,b)=\frac{\Gamma (a)\Gamma(b)}{\Gamma(a+b)}$$ + +x is supported in [0, 1] + +### Syntax + +`result = [[stdlib_stats_distribution_beta(module):beta_distribution_pdf(interface)]](x, a, b)` + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`a` has `intent(in)` and is a scalar of type real` or `complx`. + +`b`: has `intent(in)` and is a scalar of type `real` or `complx`. + +The function is elemental, i.e., all auguments could be arrays conformable to each other. All arguments must have the same type. + +### Return value + +The result is a scalar or an array, with a shape conformable to auguments, of type `real`. + +### Example + +```fortran +program demo_beta_pdf + use stdlib_stats_distribution_PRNG, onyl : random_seed + use stdlib_stats_distribution_beta, only: rbeta => beta_distribution_rvs,& + beta_pdf => beta_distribution_pdf + + implicit none + real :: x(2,3,4),aa(2,3,4),bb(2,3,4) + complx :: a, b + integer :: put, get + + put = 1234567 + call random_seed(put, get) + + print *, beta_pdf(0.65, 1.0, 1.0) + !a probability density at 0.65 with a=1.0, b=1.0 + +! 1.00000000 + + aa(:,:,:) = 2.0 + bb(:,:,:) = 1.0 + x = reshape(rbeta(2.0, 1.0, 24),[2,3,4]) ! beta random variates array + print *, beta_pdf(x,aa,bb) ! a rank 3 beta probability density array + +! [1.59333873, 1.60591197, 1.26951313, 0.825298846, 1.84633350, 0.715566635, +! 0.589380264, 1.71169996, 1.20676148, 1.79188335, 0.853198409, 1.60287392, +! 0.818829894, 1.05774701, 0.608810604, 1.40428901, 1.45229220, 1.92566359, +! 1.81786251, 1.69419682, 1.60652530, 1.87064970, 1.78721440, 0.851465702] + + a = (1.0, 1.5) + b = (1.0, 2.) + print *, beta_pdf((0.5,0.3), a, b) + ! a complex expon probability density function at (1.5,1.0) with real part of a=1.0, b=1.0 and imaginary part of a=1.5, b=2.0 + +! 1.43777180 + +end program demo_beta_pdf +``` + +## `beta_distribution_cdf` - beta cumulative distribution function + +### Status + +Experimental + +### Description + +Cumulative distribution function of the beta continuous distribution + +$$F(x)=\frac{B(x; a, b)}{B(a, b)}; where \; B(x; a, b) = \int_{0}^{x}t^{a-1}(1-t)^{b-1}dt$$ + +x is supported in [0, 1] + +### Syntax + +`result = [[stdlib_stats_distribution_beta(module):beta_distribution_cdf(interface)]](x, shape, rate)` + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`a`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`b`: has `intent(in)` and is a scalar of type `real` or `complx`. + +The function is elemental, i.e., all auguments could be arrays conformable to each other. All arguments must have the same type. + +### Return value + +The result is a scalar of type `real` with a shape conformable to auguments. + +### Example + +```fortran +program demo_beta_cdf + use stdlib_stats_distribution_PRNG, onyl : random_seed + use stdlib_stats_distribution_beta, only: rbeta => beta_distribution_rvs,& + beta_cdf => beta_distribution_cdf + + implicit none + real :: x(2,3,4),aa(2,3,4),bb(2,3,4) + complx :: a, b + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, beta_cdf(0.4, 0.5,1.0) + ! a standard beta cumulative at 0.4 with a=0.5, b=1.0 + +! 0.632455528 + + print *, beta_cdf(0.8, 1.5,2.0) + ! a cumulative at 0.8 with a a=1.5, b=2.0 + +! 0.930204272 + + aa(:,:,:) = 2.0 + bb(:,:,:) = 3.0 + x = reshape(rbeta(2.0, 3.0, 24),[2,3,4]) + !beta random variates array with a a=2.0, b=3.0 + print *, beta_cdf(x,aa,bb) ! a rank 3 standard beta cumulative array + +! [0.671738267, 0.630592465, 0.557911158, 0.157377750, 0.808665335, +! 8.00214931E-02, 0.118469201, 0.868274391, 0.268321663, 0.850062668, +! 7.99631923E-02, 0.756867588, 0.201488778, 0.228500694, 0.100621507, +! 0.412083119, 0.480990469, 0.831927001, 0.791745305, 0.654913783, +! 0.528246164, 0.275093734, 0.757254481, 0.212538764] + + a = (.7, 2.1) + b = (0.5,1.0) + print *, beta_cdf((0.5,0.5),a,b) + !complex beta cumulative distribution at (0.5,0.5) with real part of a=0.7,b=0.5 and imaginary part of a=2.1,b=1.0 + +! 9.32183489E-02 + +end program demo_beta_cdf + +``` From fe4eb31c97670ecbf377271319d91e3152043ed0 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:33:20 -0500 Subject: [PATCH 15/39] Update index.md --- doc/specs/index.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/specs/index.md b/doc/specs/index.md index 91284c2df..0b75d8600 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -17,6 +17,8 @@ This is and index/directory of the specifications (specs) for each new module/fe - [optval](./stdlib_optval.html) - Fallback value for optional arguments - [quadrature](./stdlib_quadrature.html) - Numerical integration - [stats](./stdlib_stats.html) - Descriptive Statistics + - [stats_distribution_beta](./stdlib_stats_distribution_beta.html) - Beta Distribution + ## Missing specs From 8e3b71541be3c53c17df76821da91b0d63eb2c99 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:45:01 -0500 Subject: [PATCH 16/39] Update Makefile.manual --- src/Makefile.manual | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 59af055b6..eacc97dd8 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -14,7 +14,9 @@ SRC = f18estop.f90 \ stdlib_quadrature_trapz.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ - stdlib_stats_moment.f90 \ + stdlib_stats_moment_all.f90 \ + stdlib_stats_moment_mask.f90 \ + stdlib_stats_moment_scalar.f90 \ stdlib_stats_var.f90 \ stdlib_stats_distribution_PRNG.f90 \ stdlib_stats_distribution_uniform.f90 \ From 17013cd7888ff5f7ffacab84fec5747b2acff8d7 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:49:05 -0500 Subject: [PATCH 17/39] Update Makefile.manual --- src/Makefile.manual | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Makefile.manual b/src/Makefile.manual index eacc97dd8..c834c826e 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -14,6 +14,7 @@ SRC = f18estop.f90 \ stdlib_quadrature_trapz.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ + stdlib_stats_moment.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ From afb881e8bc1fd6d669e4141202e5ecb46a812f4b Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:52:40 -0500 Subject: [PATCH 18/39] Update Makefile.manual --- src/Makefile.manual | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/Makefile.manual b/src/Makefile.manual index 872f704c0..dd6f12708 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -15,6 +15,9 @@ SRC = f18estop.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ stdlib_stats_moment.f90 \ + stdlib_stats_moment_all.f90 \ + stdlib_stats_moment_mask.f90 \ + stdlib_stats_moment_scalar.f90 \ stdlib_stats_var.f90 LIB = libstdlib.a @@ -79,4 +82,7 @@ stdlib_quadrature.f90: stdlib_quadrature.fypp stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp +stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp +stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp +stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp From 8abb168c53e184b82b2b1976d732d145a3b80d21 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 18:00:43 -0500 Subject: [PATCH 19/39] Update CMakeLists.txt --- src/CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 02604959e..e2caa0bbc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,6 +14,9 @@ set(fppFiles stdlib_stats_cov.fypp stdlib_stats_mean.fypp stdlib_stats_moment.fypp + stdlib_stats_moment_all.fypp + stdlib_stats_moment_mask.fypp + stdlib_stats_moment_scalar.fypp stdlib_stats_var.fypp stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp From 157dba80e588a64f1e85505ddf11d9fe6c9979c0 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 18:19:44 -0500 Subject: [PATCH 20/39] Update CMakeLists.txt --- src/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e2caa0bbc..8cdccdf68 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,7 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp + stdlib_stats_distribution_PRNG.fypp ) From d7643cbc5d9ad4f25d9346e5cdc27a9508a49ab2 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 18:21:31 -0500 Subject: [PATCH 21/39] Update Makefile.manual --- src/Makefile.manual | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index dd6f12708..302d1e21b 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,7 +18,8 @@ SRC = f18estop.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 + stdlib_stats_var.f90 \ + stdlib_stats_distribution_PRNG.f90 LIB = libstdlib.a @@ -70,6 +71,7 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_distribution.PRNG.o: stdlib_kinds.o stdlib_error.o # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp @@ -86,3 +88,4 @@ stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp +stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp From 5651f5e62f92caa30326fe888f1b97616a9a90f6 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 18:22:09 -0500 Subject: [PATCH 22/39] Update CMakeLists.txt --- src/tests/stats/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 36ffc7aeb..38f9bb84b 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -5,6 +5,7 @@ ADDTEST(moment) ADDTEST(rawmoment) ADDTEST(var) ADDTEST(varn) +ADDTEST(distribution_PRNG) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) From 5afcba5801a34115ca638eb9cc95804af73d7f16 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 18:26:08 -0500 Subject: [PATCH 23/39] Update CMakeLists.txt --- src/tests/stats/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 38f9bb84b..36ffc7aeb 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -5,7 +5,6 @@ ADDTEST(moment) ADDTEST(rawmoment) ADDTEST(var) ADDTEST(varn) -ADDTEST(distribution_PRNG) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) From 4cb0041d95e536c2fc325b828ac083b086286d4f Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 18:26:30 -0500 Subject: [PATCH 24/39] Update CMakeLists.txt --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8cdccdf68..1704e12ab 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,7 +21,7 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - stdlib_stats_distribution_PRNG.fypp + ) From 91c1ad42b9abfbcc4544bb1b977467368c86846a Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 18:27:11 -0500 Subject: [PATCH 25/39] Update Makefile.manual --- src/Makefile.manual | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 302d1e21b..dd6f12708 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,8 +18,7 @@ SRC = f18estop.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 \ - stdlib_stats_distribution_PRNG.f90 + stdlib_stats_var.f90 LIB = libstdlib.a @@ -71,7 +70,6 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o -stdlib_stats_distribution.PRNG.o: stdlib_kinds.o stdlib_error.o # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp @@ -88,4 +86,3 @@ stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp From 04dbbd0ee24b721999eac42c6de5eca4f5c2e1ce Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 18:54:59 -0500 Subject: [PATCH 26/39] Update CMakeLists.txt --- src/CMakeLists.txt | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2af329d5a..f17389d56 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,12 +18,7 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - stdlib_stats_distribution_PRNG.fypp - stdlib_stats_distribution_uniform.fypp - stdlib_stats_distribution_normal.fypp - stdlib_stats_distribution_special.fypp - stdlib_stats_distribution_gamma.fypp - stdlib_stats_distribution_beta.fypp + ) From 3e42dc9a905a74235c6f99908316f2522340f248 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 18:56:42 -0500 Subject: [PATCH 27/39] Update Makefile.manual --- src/Makefile.manual | 43 +++---------------------------------------- 1 file changed, 3 insertions(+), 40 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index c834c826e..8b54c2144 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,14 +18,8 @@ SRC = f18estop.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 \ - stdlib_stats_distribution_PRNG.f90 \ - stdlib_stats_distribution_uniform.f90 \ - stdlib_stats_distribution_normal.f90 \ - stdlib_stats_distribution_special.f90 \ - stdlib_stats_distribution_gamma.f90 \ - stdlib_stats_distribution_beta.f90 - + stdlib_stats_var.f90 + LIB = libstdlib.a @@ -76,32 +70,7 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o -stdlib_stats_distribution_PRNG.o: stdlib_kinds.o -stdlib_stats_distribution_uniform.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o -stdlib_stats_distribution_normal.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o \ - stdlib_stats_distribution_uniform.o -stdlib_stats_distribution_special.o: \ - stdlib_kinds.o \ - stdlib_error.o -stdlib_stats_distribution_gamma.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_uniform.o \ - stdlib_stats_distribution_normal.o \ - stdlib_stats_distribution_special.o -stdlib_stats_distribution_beta.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_uniform.o \ - stdlib_stats_distribution_gamma.o \ - stdlib_stats_distribution_special.o - + # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp @@ -114,9 +83,3 @@ stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp -stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp -stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp -stdlib_stats_distribution_special.f90: stdlib_stats_distribution_special.fypp -stdlib_stats_distribution_gamma.f90: stdlib_stats_distribution_gamma.fypp -stdlib_stats_distribution_beta.f90: stdlib_stats_distribution_beta.fypp From 461908109f9ee4f7f25af6721a050cb106592727 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 18:59:01 -0500 Subject: [PATCH 28/39] Update Makefile.manual --- src/Makefile.manual | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index dd6f12708..73aabe79c 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,7 +18,13 @@ SRC = f18estop.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 + stdlib_stats_var.f90 \ + stdlib_stats_distribution_PRNG.f90 \ + stdlib_stats_distribution_uniform.f90 \ + stdlib_stats_distribution_normal.f90 \ + stdlib_stats_distribution_special.f90 \ + stdlib_stats_distribution_gamma.f90 \ + stdlib_stats_distribution_beta.f90 LIB = libstdlib.a @@ -70,6 +76,31 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_distribution_PRNG.o: stdlib_kinds.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o +stdlib_stats_distribution_normal.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o \ + stdlib_stats_distribution_uniform.o +stdlib_stats_distribution_special.o: \ + stdlib_kinds.o \ + stdlib_error.o +stdlib_stats_distribution_gamma.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_uniform.o \ + stdlib_stats_distribution_normal.o \ + stdlib_stats_distribution_special.o +stdlib_stats_distribution_beta.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_uniform.o \ + stdlib_stats_distribution_gamma.o \ + stdlib_stats_distribution_special.o # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp @@ -86,3 +117,9 @@ stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp +stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp +stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp +stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp +stdlib_stats_distribution_special.f90: stdlib_stats_distribution_special.fypp +stdlib_stats_distribution_gamma.f90: stdlib_stats_distribution_gamma.fypp +stdlib_stats_distribution_beta.f90: stdlib_stats_distribution_beta.fypp From e6c80b824b28618c74746d1caa04870c61ff146f Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 19:01:16 -0500 Subject: [PATCH 29/39] Update CMakeLists.txt --- src/CMakeLists.txt | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1704e12ab..a8dbc803f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,7 +21,12 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - + stdlib_stats_distribution_PRNG.fypp + stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp + stdlib_stats_distribution_special.fypp + stdlib_stats_distribution_gamma.fypp + stdlib_stats_distribution_beta.fypp ) From 3e89cdc6a7a24d1fa25a46ce27b5932f32b5561a Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 11:47:43 -0500 Subject: [PATCH 30/39] Update Makefile.manual From c7cf916ede4cce55059fb6d3a98e18260951d545 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 20:23:20 -0500 Subject: [PATCH 31/39] Update CMakeLists.txt --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 256bb6e71..2429f555a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,7 +21,7 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - + stdlib_stats_distribution_PRNG.fypp ) From ba70104445d27a714d414b5cdca3849b66308cd1 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 20:23:38 -0500 Subject: [PATCH 32/39] Update CMakeLists.txt --- src/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2429f555a..02413b415 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,7 +21,6 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - stdlib_stats_distribution_PRNG.fypp ) From 1a81a484195a138febeba2c8b15b29040cefb7a6 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 21:13:50 -0500 Subject: [PATCH 33/39] Update CMakeLists.txt --- src/CMakeLists.txt | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a8dbc803f..1704e12ab 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,12 +21,7 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - stdlib_stats_distribution_PRNG.fypp - stdlib_stats_distribution_uniform.fypp - stdlib_stats_distribution_normal.fypp - stdlib_stats_distribution_special.fypp - stdlib_stats_distribution_gamma.fypp - stdlib_stats_distribution_beta.fypp + ) From bc1e8a8d75115f45692dc5cdf7428be331a9c619 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 21:14:23 -0500 Subject: [PATCH 34/39] Update Makefile.manual --- src/Makefile.manual | 127 +++++++++++++++++++------------------------- 1 file changed, 54 insertions(+), 73 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 73aabe79c..9351a374a 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -1,35 +1,35 @@ +SRCFYPP =\ + stdlib_bitsets_64.fypp \ + stdlib_bitsets_large.fypp \ + stdlib_bitsets.fypp \ + stdlib_io.fypp \ + stdlib_linalg.fypp \ + stdlib_linalg_diag.fypp \ + stdlib_optval.fypp \ + stdlib_quadrature.fypp \ + stdlib_quadrature_trapz.fypp \ + stdlib_quadrature_simps.fypp \ + stdlib_stats.fypp \ + stdlib_stats_corr.fypp \ + stdlib_stats_cov.fypp \ + stdlib_stats_mean.fypp \ + stdlib_stats_moment.fypp \ + stdlib_stats_moment_all.fypp \ + stdlib_stats_moment_mask.fypp \ + stdlib_stats_moment_scalar.fypp \ + stdlib_stats_var.fypp + SRC = f18estop.f90 \ stdlib_ascii.f90 \ - stdlib_bitsets.f90 \ - stdlib_bitsets_64.f90 \ - stdlib_bitsets_large.f90 \ stdlib_error.f90 \ - stdlib_io.f90 \ stdlib_kinds.f90 \ - stdlib_linalg.f90 \ - stdlib_linalg_diag.f90 \ stdlib_logger.f90 \ - stdlib_optval.f90 \ - stdlib_quadrature.f90 \ - stdlib_quadrature_trapz.f90 \ - stdlib_stats.f90 \ - stdlib_stats_mean.f90 \ - stdlib_stats_moment.f90 \ - stdlib_stats_moment_all.f90 \ - stdlib_stats_moment_mask.f90 \ - stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 \ - stdlib_stats_distribution_PRNG.f90 \ - stdlib_stats_distribution_uniform.f90 \ - stdlib_stats_distribution_normal.f90 \ - stdlib_stats_distribution_special.f90 \ - stdlib_stats_distribution_gamma.f90 \ - stdlib_stats_distribution_beta.f90 + $(SRCGEN) LIB = libstdlib.a - +SRCGEN = $(SRCFYPP:.fypp=.f90) OBJS = $(SRC:.f90=.o) MODS = $(OBJS:.o=.mod) SMODS = $(OBJS:.o=*.smod) @@ -42,12 +42,12 @@ $(LIB): $(OBJS) ar rcs $@ $(OBJS) clean: - $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) + $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) $(SRCGEN) %.o: %.f90 $(FC) $(FFLAGS) -c $< -%.f90: %.fypp +$(SRCGEN): %.f90: %.fypp common.fypp fypp $(FYPPFLAGS) $< $@ # Fortran module dependencies @@ -60,10 +60,32 @@ stdlib_io.o: \ stdlib_error.o \ stdlib_optval.o \ stdlib_kinds.o -stdlib_linalg_diag.o: stdlib_kinds.o +stdlib_linalg.o: \ + stdlib_kinds.o +stdlib_linalg_diag.o: \ + stdlib_linalg.o \ + stdlib_kinds.o stdlib_logger.o: stdlib_ascii.o stdlib_optval.o stdlib_optval.o: stdlib_kinds.o stdlib_quadrature.o: stdlib_kinds.o +stdlib_quadrature_simps.o: \ + stdlib_quadrature.o \ + stdlib_error.o \ + stdlib_kinds.o +stdlib_quadrature_trapz.o: \ + stdlib_quadrature.o \ + stdlib_error.o \ + stdlib_kinds.o +stdlib_stats.o: \ + stdlib_kinds.o +stdlib_stats_corr.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_cov.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o stdlib_stats_mean.o: \ stdlib_optval.o \ stdlib_kinds.o \ @@ -72,54 +94,13 @@ stdlib_stats_moment.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_moment_all.o: \ + stdlib_stats_moment.o +stdlib_stats_moment_mask.o: \ + stdlib_stats_moment.o +stdlib_stats_moment_scalar.o: \ + stdlib_stats_moment.o stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o -stdlib_stats_distribution_PRNG.o: stdlib_kinds.o -stdlib_stats_distribution_uniform.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o -stdlib_stats_distribution_normal.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o \ - stdlib_stats_distribution_uniform.o -stdlib_stats_distribution_special.o: \ - stdlib_kinds.o \ - stdlib_error.o -stdlib_stats_distribution_gamma.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_uniform.o \ - stdlib_stats_distribution_normal.o \ - stdlib_stats_distribution_special.o -stdlib_stats_distribution_beta.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_uniform.o \ - stdlib_stats_distribution_gamma.o \ - stdlib_stats_distribution_special.o - -# Fortran sources that are built from fypp templates -stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp -stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp -stdlib_bitsets.f90: stdlib_bitsets.fypp -stdlib_io.f90: stdlib_io.fypp -stdlib_linalg.f90: stdlib_linalg.fypp -stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp -stdlib_quadrature.f90: stdlib_quadrature.fypp -stdlib_stats.f90: stdlib_stats.fypp -stdlib_stats_mean.f90: stdlib_stats_mean.fypp -stdlib_stats_moment.f90: stdlib_stats_moment.fypp -stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp -stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp -stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp -stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp -stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp -stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp -stdlib_stats_distribution_special.f90: stdlib_stats_distribution_special.fypp -stdlib_stats_distribution_gamma.f90: stdlib_stats_distribution_gamma.fypp -stdlib_stats_distribution_beta.f90: stdlib_stats_distribution_beta.fypp From f2c07cef3a1e96420d5a6e07559a09c0f7053030 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 21:19:44 -0500 Subject: [PATCH 35/39] Update Makefile.manual --- src/Makefile.manual | 37 ++++++++++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 9351a374a..ed1805f48 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -17,7 +17,13 @@ SRCFYPP =\ stdlib_stats_moment_all.fypp \ stdlib_stats_moment_mask.fypp \ stdlib_stats_moment_scalar.fypp \ - stdlib_stats_var.fypp + stdlib_stats_var.fypp \ + stdlib_stats_distribution_PRNG.fypp \ + stdlib_stats_distribution_uniform.fypp \ + stdlib_stats_distribution_normal.fypp \ + stdlib_stats_distribution_gamma.fypp \ + stdlib_stats_distribution_special.fypp \ + stdlib_stats_distribution_beta.fypp SRC = f18estop.f90 \ stdlib_ascii.f90 \ @@ -104,3 +110,32 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_distribution_PRNG.o: \ + stdlib_kinds.o \ + stdlib_error.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o +stdlib_stats_distribution_normal.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o \ + stdlib_stats_distribution_uniform.o +stdlib_stats_distribution_special.o: \ + stdlib_kinds.o \ + stdlib_error.o +stdlib_stats_distribution_gamma.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o \ + stdlib_stats_distribution_uniform.o \ + stdlib_stats_distribution_normal.o \ + stdlib_stats_distribution_special.o +stdlib_stats_distribution_beta.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o \ + stdlib_stats_distribution_uniform.o \ + stdlib_stats_distribution_special.o \ + stdlib_stats_distribution_gamma.o From 0741c17d2df8e48291c9ba678b64d6d18cd5ae20 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 21:20:15 -0500 Subject: [PATCH 36/39] Add files via upload --- doc/specs/stdlib_stats_distribution_beta.md | 63 +++++++++++---------- 1 file changed, 33 insertions(+), 30 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_beta.md b/doc/specs/stdlib_stats_distribution_beta.md index 3d74c11d3..78786cc9e 100644 --- a/doc/specs/stdlib_stats_distribution_beta.md +++ b/doc/specs/stdlib_stats_distribution_beta.md @@ -25,15 +25,15 @@ With three auguments, the function return a rank one array of beta distributed r ### Arguments -`a` : has `intent(in)` ans is a scalar of type `real` or `complx`. +`a` : has `intent(in)` ans is a scalar of type `real` or `complex`. -`b`: has `intent(in)` and is a scalar of type `real` or `complx`. +`b`: has `intent(in)` and is a scalar of type `real` or `complex`. `array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. ### Return value -The result is a scalar or rank one array, with a size of `array_size`, of type `real` or `complx`. +The result is a scalar or rank one array, with a size of `array_size`, of type `real` or `complex`. ### Example @@ -63,22 +63,23 @@ program demo_beta_rvs print *, rbeta(aa, bb) !a rank 3 array of 24 beta random variates with a=0.8, b=0.6 -! [0.251766384, 0.578202426, 0.539138556, 0.210130826, 0.908130825, -! 0.880996943, 9.49194748E-03, 0.945992589, 0.290732056, 0.803920329, -! 7.64303207E-02, 0.943150401, 0.927998245, 0.831781328, 0.671169102, -! 0.983966410, 0.289062619, 0.801237404, 0.891931713, 0.897902310, -! 0.845606744, 1.50359496E-02, 0.913162351, 0.915781260] +! 0.251766384 0.578202426 0.539138556 0.210130826 0.908130825 +! 0.880996943 9.49194748E-03 0.945992589 0.290732056 0.803920329 +! 7.64303207E-02 0.943150401 0.927998245 0.831781328 0.671169102 +! 0.983966410 0.289062619 0.801237404 0.891931713 0.897902310 +! 0.845606744 1.50359496E-02 0.913162351 0.915781260 print *, rbeta(1.2,0.7,10) ! an array of 10 random variates with a=1.2, b=0.7 -! [3.69704030E-02, 0.748639643, 0.896924615, 0.350249082, 0.813054740, -! 0.448072791, 9.39368978E-02, 0.475665420, 0.567661405, 0.893835664] +! 3.69704030E-02 0.748639643 0.896924615 0.350249082 0.813054740 +! 0.448072791 9.39368978E-02 0.475665420 0.567661405 0.893835664 a = (3.0, 4.0) b = (2.0, 0.7) print *, rbeta(a, b) - !single complex beta random variate with real part of a = 3.0, b=2.0; imagainary part of a=4.0, b=0.7 + !single complex beta random variate with real part of a = 3.0, b=2.0; + !imagainary part of a=4.0, b=0.7 ! (0.730923295,0.878909111) @@ -105,11 +106,11 @@ x is supported in [0, 1] ### Arguments -`x`: has `intent(in)` and is a scalar of type `real` or `complx`. +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. -`a` has `intent(in)` and is a scalar of type real` or `complx`. +`a` has `intent(in)` and is a scalar of type real` or `complex`. -`b`: has `intent(in)` and is a scalar of type `real` or `complx`. +`b`: has `intent(in)` and is a scalar of type `real` or `complex`. The function is elemental, i.e., all auguments could be arrays conformable to each other. All arguments must have the same type. @@ -127,7 +128,7 @@ program demo_beta_pdf implicit none real :: x(2,3,4),aa(2,3,4),bb(2,3,4) - complx :: a, b + complex :: a, b integer :: put, get put = 1234567 @@ -143,15 +144,16 @@ program demo_beta_pdf x = reshape(rbeta(2.0, 1.0, 24),[2,3,4]) ! beta random variates array print *, beta_pdf(x,aa,bb) ! a rank 3 beta probability density array -! [1.59333873, 1.60591197, 1.26951313, 0.825298846, 1.84633350, 0.715566635, -! 0.589380264, 1.71169996, 1.20676148, 1.79188335, 0.853198409, 1.60287392, -! 0.818829894, 1.05774701, 0.608810604, 1.40428901, 1.45229220, 1.92566359, -! 1.81786251, 1.69419682, 1.60652530, 1.87064970, 1.78721440, 0.851465702] +! 1.59333873 1.60591197 1.26951313 0.825298846 1.84633350 0.715566635 +! 0.589380264 1.71169996 1.20676148 1.79188335 0.853198409 1.60287392 +! 0.818829894 1.05774701 0.608810604 1.40428901 1.45229220 1.92566359 +! 1.81786251 1.69419682 1.60652530 1.87064970 1.78721440 0.851465702 a = (1.0, 1.5) b = (1.0, 2.) print *, beta_pdf((0.5,0.3), a, b) - ! a complex expon probability density function at (1.5,1.0) with real part of a=1.0, b=1.0 and imaginary part of a=1.5, b=2.0 + ! a complex expon probability density function at (1.5,1.0) with real part of + !a=1.0, b=1.0 and imaginary part of a=1.5, b=2.0 ! 1.43777180 @@ -178,11 +180,11 @@ x is supported in [0, 1] ### Arguments -`x`: has `intent(in)` and is a scalar of type `real` or `complx`. +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. -`a`: has `intent(in)` and is a scalar of type `real` or `complx`. +`a`: has `intent(in)` and is a scalar of type `real` or `complex`. -`b`: has `intent(in)` and is a scalar of type `real` or `complx`. +`b`: has `intent(in)` and is a scalar of type `real` or `complex`. The function is elemental, i.e., all auguments could be arrays conformable to each other. All arguments must have the same type. @@ -200,7 +202,7 @@ program demo_beta_cdf implicit none real :: x(2,3,4),aa(2,3,4),bb(2,3,4) - complx :: a, b + complex :: a, b integer :: seed_put, seed_get seed_put = 1234567 @@ -222,16 +224,17 @@ program demo_beta_cdf !beta random variates array with a a=2.0, b=3.0 print *, beta_cdf(x,aa,bb) ! a rank 3 standard beta cumulative array -! [0.671738267, 0.630592465, 0.557911158, 0.157377750, 0.808665335, -! 8.00214931E-02, 0.118469201, 0.868274391, 0.268321663, 0.850062668, -! 7.99631923E-02, 0.756867588, 0.201488778, 0.228500694, 0.100621507, -! 0.412083119, 0.480990469, 0.831927001, 0.791745305, 0.654913783, -! 0.528246164, 0.275093734, 0.757254481, 0.212538764] +! 0.671738267 0.630592465 0.557911158 0.157377750 0.808665335 +! 8.00214931E-02 0.118469201 0.868274391 0.268321663 0.850062668 +! 7.99631923E-02 0.756867588 0.201488778 0.228500694 0.100621507 +! 0.412083119 0.480990469 0.831927001 0.791745305 0.654913783 +! 0.528246164 0.275093734 0.757254481 0.212538764 a = (.7, 2.1) b = (0.5,1.0) print *, beta_cdf((0.5,0.5),a,b) - !complex beta cumulative distribution at (0.5,0.5) with real part of a=0.7,b=0.5 and imaginary part of a=2.1,b=1.0 + !complex beta cumulative distribution at (0.5,0.5) with real part of + !a=0.7,b=0.5 and imaginary part of a=2.1,b=1.0 ! 9.32183489E-02 From 2358f4eeb601760ac5e15c51d3d06badfa8c7e57 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 21:25:22 -0500 Subject: [PATCH 37/39] Add files via upload --- src/stdlib_stats_distribution_PRNG.fypp | 75 ++-------------------- src/stdlib_stats_distribution_beta.fypp | 18 +++--- src/stdlib_stats_distribution_gamma.fypp | 51 ++++++++------- src/stdlib_stats_distribution_normal.fypp | 16 ++--- src/stdlib_stats_distribution_special.fypp | 46 ++++++------- src/stdlib_stats_distribution_uniform.fypp | 46 +++++++------ 6 files changed, 101 insertions(+), 151 deletions(-) diff --git a/src/stdlib_stats_distribution_PRNG.fypp b/src/stdlib_stats_distribution_PRNG.fypp index d1bda107a..3fdbf0438 100644 --- a/src/stdlib_stats_distribution_PRNG.fypp +++ b/src/stdlib_stats_distribution_PRNG.fypp @@ -1,16 +1,16 @@ #:include "common.fypp" module stdlib_stats_distribution_PRNG use stdlib_kinds, only: int8, int16, int32, int64 + use stdlib_error implicit none private integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64) - integer(int64), save :: st(4), si = 614872703977525537_int64 + integer(int64), save :: st(4) ! internal states for xoshiro256ss function + integer(int64), save :: si = 614872703977525537_int64 ! default seed value logical, save :: seed_initialized = .false. public :: random_seed public :: dist_rand - public :: jump - public :: long_jump interface dist_rand @@ -51,6 +51,8 @@ module stdlib_stats_distribution_PRNG integer :: k k = MAX_INT_BIT_SIZE - bit_size(n) + if(k < 0) call error_stop("Error(dist_rand): Integer bit size is" & + //" greater than 64bit") res = shiftr(xoshiro256ss( ), k) end function dist_rand_${t1[0]}$${k1}$ @@ -96,71 +98,6 @@ module stdlib_stats_distribution_PRNG end function rol64 - subroutine jump - ! This is the jump function for the xoshiro256ss generator. It is equivalent - ! to 2^128 calls to xoshiro256ss(); it can be used to generate 2^128 - ! non-overlapping subsequences for parallel computations. - ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) - ! http://prng.di.unimi.it/xoshiro256starstar.c - ! - ! Fortran 90 version translated from C by Jim-215-Fisher - integer(int64) :: jp(4) = [1733541517147835066_int64, & - -3051731464161248980_int64, & - -6244198995065845334_int64, & - 4155657270789760540_int64] - integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 - integer :: i, j, k - - do i = 1, 4 - do j = 1, 64 - if(iand(jp(i), shiftl(c, j - 1)) /= 0) then - s1 = ieor(s1, st(1)) - s2 = ieor(s2, st(2)) - s3 = ieor(s3, st(3)) - s4 = ieor(s4, st(4)) - end if - k = xoshiro256ss( ) - end do - end do - st(1) = s1 - st(2) = s2 - st(3) = s3 - st(4) = s4 - end subroutine jump - - subroutine long_jump - ! This is the long-jump function for the xoshiro256ss generator. It is - ! equivalent to 2^192 calls to xoshiro256ss(); it can be used to generate - ! 2^64 starting points, from each of which jump() will generate 2^64 - ! non-overlapping subsequences for parallel distributed computations - ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) - ! http://prng.di.unimi.it/xoshiro256starstar.c - ! - ! Fortran 90 version translated from C by Jim-215-Fisher - integer(int64) :: jp(4) = [8566230491382795199_int64, & - -4251311993797857357_int64, & - 8606660816089834049_int64, & - 4111957640723818037_int64] - integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 - integer(int32) :: i, j, k - - do i = 1, 4 - do j = 1, 64 - if(iand(jp(i), shiftl(c, j - 1)) /= 0) then - s1 = ieor(s1, st(1)) - s2 = ieor(s2, st(2)) - s3 = ieor(s3, st(3)) - s4 = ieor(s4, st(4)) - end if - k = xoshiro256ss() - end do - end do - st(1) = s1 - st(2) = s2 - st(3) = s3 - st(4) = s4 - end subroutine long_jump - function splitmix64(s) result(res) ! Written in 2015 by Sebastiano Vigna (vigna@acm.org) ! This is a fixed-increment version of Java 8's SplittableRandom @@ -178,6 +115,8 @@ module stdlib_stats_distribution_PRNG data int01, int02, int03/-7046029254386353131_int64, & -4658895280553007687_int64, & -7723592293110705685_int64/ + ! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15, + ! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb if(present(s)) si = s res = si diff --git a/src/stdlib_stats_distribution_beta.fypp b/src/stdlib_stats_distribution_beta.fypp index af1707660..1ac6eb7c5 100644 --- a/src/stdlib_stats_distribution_beta.fypp +++ b/src/stdlib_stats_distribution_beta.fypp @@ -67,8 +67,8 @@ Module stdlib_stats_distribution_beta ${t1}$ :: res, x, y, xx(2) ${t1}$, parameter :: z = 0.0_${k1}$, one = 1.0_${k1}$ - if(a <= z .or. b <= z) call error_stop("Error: Beta distribution" & - //" paramters a, b must be greater than zero") + if(a <= z .or. b <= z) call error_stop("Error(beta_dist_rvs): Beta" & + //" distribution paramters a, b must be greater than zero") if( a < one .or. b < one) then do xx = uni(z, one, 2) @@ -117,9 +117,11 @@ Module stdlib_stats_distribution_beta ${t1}$ :: x, y, xx(2) real(${k1}$), parameter :: z = 0.0_${k1}$, one = 1.0_${k1}$ + if(a <= z .or. b <= z) call error_stop("Error(beta_dist_rvs_array):" & + //" Beta distribution paramters a, b must be greater than zero") + allocate(res(array_size)) - if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error: Beta" & - //" distribution paramters a, b must be greater than zero") + if( a < one .or. b < one) then do i = 1, array_size do @@ -171,8 +173,8 @@ Module stdlib_stats_distribution_beta real :: res ${t1}$ :: z = 0.0_${k1}$, one = 1.0_${k1}$ - if(a <= z .or. b <= z) call error_stop("Error: Beta distribution" & - //" parameters a, b must be greater than zero") + if(a <= z .or. b <= z) call error_stop("Error(beta_dist_pdf): Beta" & + //" distribution parameters a, b must be greater than zero") if(x == z) then if(a <= one) then res = huge(1.0) @@ -213,8 +215,8 @@ Module stdlib_stats_distribution_beta real :: res ${t1}$ :: z = 0.0_${k1}$, one = 1.0_${k1}$ - if(a <= z .or. b <= z) call error_stop("Error: Beta distribution" & - //" parameters a, b must be greater than zero") + if(a <= z .or. b <= z) call error_stop("Error(beta_dist_cdf): Beta" & + //" distribution parameters a, b must be greater than zero") if(x == z) then res = 0.0 elseif(x == one) then diff --git a/src/stdlib_stats_distribution_gamma.fypp b/src/stdlib_stats_distribution_gamma.fypp index c25d47804..488e37858 100644 --- a/src/stdlib_stats_distribution_gamma.fypp +++ b/src/stdlib_stats_distribution_gamma.fypp @@ -42,7 +42,7 @@ Module stdlib_stats_distribution_gamma !! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html# !! description)) !! - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES module procedure gamma_dist_pdf_${t1[0]}$${k1}$ #:endfor end interface gamma_distribution_pdf @@ -71,8 +71,9 @@ Module stdlib_stats_distribution_gamma ${t1}$ :: x, v, u, zz, tol = 1000 * epsilon(1.0_${k1}$) ${t1}$, save :: alpha = 0._${k1}$, d, c, sq = 0.0331_${k1}$ - if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" shape parameter must be greater than zero") + if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs): Gamma" & + //" distribution shape parameter must be greater than zero") + zz = shape if(zz < 1._${k1}$) zz = 1._${k1}$ + zz if(abs(zz - alpha) > tol) then @@ -127,10 +128,12 @@ Module stdlib_stats_distribution_gamma ${t1}$ :: x, v, u, zz, tol = 1000 * epsilon(1.0_${k1}$) ${t1}$, save :: alpha = 0._${k1}$, d, c, sq = 0.0331_${k1}$ - if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" shape parameter must be greater than zero") - if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" rate parameter must be greater than zero") + if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs): Gamma" & + //" distribution shape parameter must be greater than zero") + + if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs): Gamma" & + //" distribution rate parameter must be greater than zero") + zz = shape if(zz < 1._${k1}$) zz = 1._${k1}$ + zz if(abs(zz - alpha) > tol) then @@ -189,10 +192,12 @@ Module stdlib_stats_distribution_gamma ${t1}$, save :: alpha = 0._${k1}$, d, c, sq = 0.0331_${k1}$ integer :: i - if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" shape parameter must be greater than zero") - if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" rate parameter must be greater than zero") + if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs_array):" & + //" Gamma distribution shape parameter must be greater than zero") + + if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs_array):" & + //" Gamma distribution rate parameter must be greater than zero") + allocate(res(array_size)) zz = shape if(zz < 1._${k1}$) zz = 1._${k1}$ + zz @@ -254,12 +259,12 @@ Module stdlib_stats_distribution_gamma ${t1}$, intent(in) :: x, shape, rate real :: res - if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" rate parameter must be greaeter than zero") - if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" shape parameter must be greater than zero") - if(x <= 0.0_${k1}$) call error_stop("Error: Gamma distribution variate" & - //" x must be greater than zero") + if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" & + //" distribution rate parameter must be greaeter than zero") + if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" & + //" distribution shape parameter must be greater than zero") + if(x <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" & + //" distribution variate x must be greater than zero") if(x == 0.0_${k1}$) then if(shape <= 1.0_${k1}$) then res = huge(1.0) + 1.0 @@ -296,12 +301,12 @@ Module stdlib_stats_distribution_gamma ${t1}$, intent(in) :: x, shape, rate real :: res - if(rate <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" rate parameter must be greaeter than zero") - if(shape <= 0.0_${k1}$) call error_stop("Error: Gamma distribution" & - //" shape parameter must be greater than zero") - if(x <= 0.0_${k1}$) call error_stop("Error: Gamma distribution variate" & - //" x must be greater than zero") + if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" & + //" distribution rate parameter must be greaeter than zero") + if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" & + //" distribution shape parameter must be greater than zero") + if(x <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" & + //" distribution variate x must be greater than zero") res = ingamma(shape, rate * x) / gamma(shape) return end function gamma_dist_cdf_${t1[0]}$${k1}$ diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 8397ba08c..fecef90e8 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -163,8 +163,8 @@ Module stdlib_stats_distribution_normal ${t1}$ :: x, y integer :: hz, iz - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs): Normal" & + //" distribution scale parameter must be non-zero") if( .not. zig_norm_initialized ) call zigset iz = 0 ! original algorithm use 32bit @@ -235,8 +235,8 @@ Module stdlib_stats_distribution_normal ${t1}$ :: x, y, re integer :: hz, iz, i - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs_array):" & + //" Normal distribution scale parameter must be non-zero") if( .not. zig_norm_initialized ) call zigset allocate(res(array_size)) do i = 1, array_size @@ -311,8 +311,8 @@ Module stdlib_stats_distribution_normal real :: res ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_pdf):" & + //" Normal distribution scale parameter must be non-zero") res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & (sqrt_2_Pi * scale) return @@ -342,8 +342,8 @@ Module stdlib_stats_distribution_normal real :: res ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_cdf):" & + //" Normal distribution scale parameter must be non-zero") res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ return end function norm_dist_cdf_${t1[0]}$${k1}$ diff --git a/src/stdlib_stats_distribution_special.fypp b/src/stdlib_stats_distribution_special.fypp index ca1992f9e..1c71b0106 100644 --- a/src/stdlib_stats_distribution_special.fypp +++ b/src/stdlib_stats_distribution_special.fypp @@ -173,8 +173,8 @@ Module stdlib_stats_distribution_special real(qp) :: q, sum integer :: i - if(x <= 0._${k1}$) call error_stop("Error: Gamma function augument" & - //" must be greater than 0") + if(x <= 0._${k1}$) call error_stop("Error(l_gamma): Gamma function" & + //" augument must be greater than 0") if(x == 1.0_${k1}$ .or. x == 2.0_${k1}$) then res = 0.0_${k1}$ else @@ -197,8 +197,8 @@ Module stdlib_stats_distribution_special ${t1}$, intent(in) :: n real :: res - if(n < 0) call error_stop("Error: Factorial function augument must" & - //" be no less than 0") + if(n < 0) call error_stop("Error(l_factorial): Factorial function" & + //" augument must be no less than 0") select case(n) case (0) res = 0.0 @@ -221,8 +221,8 @@ Module stdlib_stats_distribution_special ${t2}$, intent(in) :: x ${t2}$ :: res - if(n < 0) call error_stop("Error: factorial function augument must" & - //" be no less than 0") + if(n < 0) call error_stop("Error(l_factorial): Factorial function" & + //" augument must be no less than 0") select case(n) case (0) res = 0.0_${k2}$ @@ -251,8 +251,8 @@ Module stdlib_stats_distribution_special integer :: n if(x < 0._${k1}$) then - call error_stop("Error: Incomplete gamma function with negative x" & - //" must come with integer of s") + call error_stop("Error(gpx): Incomplete gamma function with" & + //" negative x must come with integer of s") elseif(s >= x) then a = real(s, dp) g = 1.0_dp / a @@ -382,11 +382,11 @@ Module stdlib_stats_distribution_special real(dp) :: s1, y, xx, ss #:if t1[0] == "i" - if(s < 0_${k1}$) call error_stop("Error: Lower incomplete gamma" & - //" function input s value must be greater than 0") + if(s < 0_${k1}$) call error_stop("Error(ingamma_low): Lower" & + //" incomplete gamma function input s value must be greater than 0") #:else - if(s < 0._${k1}$) call error_stop("Error: Lower incomplete gamma" & - //" function input s value must be greater than 0") + if(s < 0._${k1}$) call error_stop("Error(ingamma_low): Lower" & + //" incomplete gamma function input s value must be greater than 0") #:endif xx = real(x, dp); ss = real(s, dp) if(x == 0.0_${k2}$) then @@ -471,11 +471,11 @@ Module stdlib_stats_distribution_special real(dp) :: s1, xx, ss #:if t1[0] == "i" - if(s < 0_${k1}$) call error_stop("Error: Lower incomplete gamma" & - //" function input s value must be greater than 0") + if(s < 0_${k1}$) call error_stop("Error(regamma_p): Lower incomplete" & + //" gamma function input s value must be greater than 0") #:else - if(s < 0._${k1}$) call error_stop("Error: Lower incomplete gamma" & - //" function input s value must be greater than 0") + if(s < 0._${k1}$) call error_stop("Error(regamma_p): Lower incomplete" & + //" gamma function input s value must be greater than 0") #:endif xx = real(x, dp); ss = real(s, dp) s1 = -xx + ss * log(abs(xx)) - loggamma(ss) @@ -521,8 +521,8 @@ Module stdlib_stats_distribution_special ${t1}$, intent(in) :: a, b ${t1}$ :: res - if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error: Beta" & - //" function auguments a, b values must be greater than 0") + if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error(beta):" & + //" Beta function auguments a, b values must be greater than 0") res = exp(loggamma(a) + loggamma(b) - loggamma(a+b)) return end function beta_${t1[0]}$${k1}$ @@ -536,8 +536,8 @@ Module stdlib_stats_distribution_special ${t1}$, intent(in) :: a, b ${t1}$ :: res - if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error: Beta" & - //" function auguments a, b values must be greater than 0") + if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error(l_beta):"& + //" Beta function auguments a, b values must be greater than 0") res = loggamma(a) + loggamma(b) - loggamma(a+b) return end function l_beta_${t1[0]}$${k1}$ @@ -554,9 +554,9 @@ Module stdlib_stats_distribution_special integer :: n, k real(dp) :: an, bn, g, c, d, y, s0, ak, ak2 - if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error:" & - //" Incomplete beta function auguments a, b values must be" & - //" greater than 0") + if(a <= 0._${k1}$ .or. b <= 0._${k1}$) call error_stop("Error(inbeta):"& + //" Incomplete beta function auguments a, b values must be greater" & + //" than 0") s0 = (a + 1) / (a + b + 2) an = 1.0_dp bn = 1.0_dp diff --git a/src/stdlib_stats_distribution_uniform.fypp b/src/stdlib_stats_distribution_uniform.fypp index 572caaf37..b538144f8 100644 --- a/src/stdlib_stats_distribution_uniform.fypp +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -42,7 +42,8 @@ Module stdlib_stats_distribution_uniform interface uniform_distribution_pdf !! Version experiment !! - !! Get uniform distribution probability density (pdf) for integer, real and complex variables + !! Get uniform distribution probability density (pdf) for integer, real and + !! complex variables !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# !! description)) @@ -54,7 +55,8 @@ Module stdlib_stats_distribution_uniform interface uniform_distribution_cdf !! Version experimental !! - !! Get uniform distribution cumulative distribution function (cdf) for integer, real and complex variables + !! Get uniform distribution cumulative distribution function (cdf) for + !! integer, real and complex variables !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# !! description)) !! @@ -66,7 +68,8 @@ Module stdlib_stats_distribution_uniform interface shuffle !! Version experimental !! - !! Fisher-Yates shuffle algorithm for a rank one array of integer, real and complex variables + !! Fisher-Yates shuffle algorithm for a rank one array of integer, real and + !! complex variables !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# !! description)) !! @@ -85,13 +88,14 @@ Module stdlib_stats_distribution_uniform ! 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, n integer :: zeros, bits_left, bits n = scale - if(n <= 0_${k1}$) call error_stop("Error: Uniform distribution scale" & - //" parameter must be positive") + if(n <= 0_${k1}$) call error_stop("Error(unif_dist_rvs_1): Uniform" & + //" distribution scale parameter must be positive") zeros = leadz(n) bits = bit_size(n) - zeros mask = shiftr(not(0_${k1}$), zeros) @@ -121,8 +125,8 @@ Module stdlib_stats_distribution_uniform ${t1}$, intent(in) :: loc, scale ${t1}$ :: res - if(scale == 0_${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale <= 0_${k1}$) call error_stop("Error(unif_dist_rvs): Uniform" & + //" distribution scale parameter must be positive") res = loc + unif_dist_rvs_1_${t1[0]}$${k1}$(scale) return end function unif_dist_rvs_${t1[0]}$${k1}$ @@ -153,8 +157,8 @@ Module stdlib_stats_distribution_uniform ${t1}$, intent(in) :: scale ${t1}$ :: res - if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs_1): " & + //"Uniform distribution scale parameter must be non-zero") res = scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) return end function unif_dist_rvs_1_${t1[0]}$${k1}$ @@ -169,8 +173,8 @@ Module stdlib_stats_distribution_uniform ${t1}$, intent(in) :: loc, scale ${t1}$ :: res - if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs): " & + //"Uniform distribution scale parameter must be non-zero") res = loc + scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) return end function unif_dist_rvs_${t1[0]}$${k1}$ @@ -187,8 +191,8 @@ Module stdlib_stats_distribution_uniform ${t1}$ :: res real(${k1}$) :: r1, r2, tr, ti - if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & - //" distribution scale parameter must be non-zero") + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(uni_dist_" & + //"rvs_1): Uniform distribution scale parameter must be non-zero") r1 = unif_dist_rvs_0_r${k1}$( ) if(real(scale) == 0.0_${k1}$) then ti = aimag(scale) * r1 @@ -219,8 +223,8 @@ Module stdlib_stats_distribution_uniform ${t1}$ :: res real(${k1}$) :: r1, r2, tr, ti - if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & - //" distribution scale parameter must be non-zero") + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(uni_dist_" & + //"rvs): Uniform distribution scale parameter must be non-zero") r1 = unif_dist_rvs_0_r${k1}$( ) if(real(scale) == 0.0_${k1}$) then tr = real(loc) @@ -249,8 +253,8 @@ Module stdlib_stats_distribution_uniform integer :: i, zeros, bits_left, bits n = scale - if(n == 0_${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(n == 0_${k1}$) call error_stop("Error(unif_dist_rvs_array): Uniform" & + //" distribution scale parameter must be non-zero") allocate(res(array_size)) zeros = leadz(n) bits = bit_size(n) - zeros @@ -287,8 +291,8 @@ Module stdlib_stats_distribution_uniform integer :: i - if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs_array):" & + //" Uniform distribution scale parameter must be non-zero") allocate(res(array_size)) do i = 1, array_size tmp = shiftr(dist_rand(INT_ONE), 11) @@ -311,8 +315,8 @@ Module stdlib_stats_distribution_uniform integer :: i - if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & - //" distribution scale parameter must be non-zero") + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(unif_dist_"& + //"rvs_array): Uniform distribution scale parameter must be non-zero") allocate(res(array_size)) do i = 1, array_size tmp = shiftr(dist_rand(INT_ONE), 11) From 7a5e9eae5c75438a6a94a9b5119f51a9ca6bd31b Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 21:28:24 -0500 Subject: [PATCH 38/39] Update CMakeLists.txt --- src/CMakeLists.txt | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 02413b415..e620a3a44 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,12 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp + stdlib_stats_distribution_PRNG.fypp + stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp + stdlib_stats_distribution_special.fypp + stdlib_stats_distribution_gamma.fypp + stdlib_stats_distribution_beta.fypp ) From c1f87bc86cc56001d18b411845c6bdab814a5a53 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 21:39:10 -0500 Subject: [PATCH 39/39] Update Makefile.manual --- src/Makefile.manual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index ed1805f48..a448317fa 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -122,7 +122,7 @@ stdlib_stats_distribution_normal.o: \ stdlib_error.o \ stdlib_stats_distribution_PRNG.o \ stdlib_stats_distribution_uniform.o -stdlib_stats_distribution_special.o: \ +stdlib_stats_distribution_special.o: \ stdlib_kinds.o \ stdlib_error.o stdlib_stats_distribution_gamma.o: \