diff --git a/doc/specs/index.md b/doc/specs/index.md index 6ea78b52e..fe93ec8c6 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -19,6 +19,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 diff --git a/doc/specs/stdlib_stats_distribution_beta.md b/doc/specs/stdlib_stats_distribution_beta.md new file mode 100644 index 000000000..78786cc9e --- /dev/null +++ b/doc/specs/stdlib_stats_distribution_beta.md @@ -0,0 +1,243 @@ +--- + +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 `complex`. + +`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 `complex`. + +### 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 `complex`. + +`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 `complex`. + +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) + complex :: 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 `complex`. + +`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 `complex`. + +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) + complex :: 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 + +``` 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 ) diff --git a/src/Makefile.manual b/src/Makefile.manual index 9351a374a..a448317fa 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 diff --git a/src/stdlib_stats_distribution_PRNG.fypp b/src/stdlib_stats_distribution_PRNG.fypp new file mode 100644 index 000000000..3fdbf0438 --- /dev/null +++ b/src/stdlib_stats_distribution_PRNG.fypp @@ -0,0 +1,151 @@ +#: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) ! 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 + + + 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) + 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}$ + + #: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 + + + 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/ + ! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15, + ! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb + + 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..1ac6eb7c5 --- /dev/null +++ b/src/stdlib_stats_distribution_beta.fypp @@ -0,0 +1,243 @@ +#: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_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) + 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}$ + + 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 < 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_dist_pdf): 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_dist_cdf): 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..488e37858 --- /dev/null +++ b/src/stdlib_stats_distribution_gamma.fypp @@ -0,0 +1,328 @@ +#: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_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 + 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_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 + 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_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 + 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_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 + 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_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}$ + + #: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..fecef90e8 --- /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(norm_dist_rvs): 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(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 + 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(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 + 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(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}$ + + #: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..1c71b0106 --- /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(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 + 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(l_factorial): 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(l_factorial): 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(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 + 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(ingamma_low): Lower" & + //" incomplete gamma function input s value must be greater than 0") + #:else + 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 + 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(regamma_p): Lower incomplete" & + //" gamma function input s value must be greater than 0") + #:else + 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) + 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):" & + //" 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(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}$ + + #: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(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 + 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..b538144f8 --- /dev/null +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -0,0 +1,486 @@ +#: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(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) + 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(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}$ + + #: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(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}$ + + #: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(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}$ + + #: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(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 + 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(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) + 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(unif_dist_rvs_array): 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(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) + 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(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) + 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 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) 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 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