diff --git a/doc/specs/index.md b/doc/specs/index.md index f709fb2ca..30e45c25e 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -35,6 +35,7 @@ This is an index/directory of the specifications (specs) for each new module/fea - [stats_distributions_uniform](./stdlib_stats_distribution_uniform.html) - Uniform Probability Distribution - [stats_distributions_normal](./stdlib_stats_distribution_normal.html) - Normal Probability Distribution - [stats_distributions_exponential](./stdlib_stats_distribution_exponential.html) - Exponential Probability Distribution + - [stats_distributions_beta](./stdlib_stats_distribution_beta.html) - Beta Probability Distribution - [stats_distributions_gamma](./stdlib_stats_distribution_gamma.html) - Gamma Probability Distribution - [string\_type](./stdlib_string_type.html) - Basic string support - [stringlist_type](./stdlib_stringlist_type.html) - 1-Dimensional list of strings diff --git a/doc/specs/stdlib_stats_distribution_beta.md b/doc/specs/stdlib_stats_distribution_beta.md new file mode 100644 index 000000000..d0f4ae7a1 --- /dev/null +++ b/doc/specs/stdlib_stats_distribution_beta.md @@ -0,0 +1,167 @@ +--- +title: stats_distribution_beta +--- + +# Statistical Distributions -- Beta Distribution Module + +[TOC] + +## `rvs_beta` - beta distribution random variates + +### Status + +Experimental + +### Description + +The beta distribution is a continuous probability distribution defined on the interval [0, 1], widely used for modeling random variables that represent proportions, probabilities, and other bounded quantities. It is defined by two shape parameters (\\(a\\) and \\(b\\)) that control the distribution's form. + +With two arguments (a, b), the function returns a random sample from the beta distribution \\(\\text{Beta}(a, b)\\). + +The optional `loc` parameter specifies the location (shift) of the distribution. + +With three or more arguments including `array_size`, the function returns a rank-1 array of beta distributed random variates. + +For complex shape parameters, the real and imaginary parts are sampled independently of each other. + +@note +For shape parameters less than 1, the function uses a uniform method. For parameters greater than or equal to 1, it uses the gamma ratio method[^1], where \\(X \\sim \\text{Beta}(a,b)\\) is generated as \\(X = \\frac{Y_1}{Y_1 + Y_2}\\) where \\(Y_1 \\sim \\Gamma(a,1)\\) and \\(Y_2 \\sim \\Gamma(b,1)\\). + +### Syntax + +`result = [[stdlib_stats_distribution_beta(module):rvs_beta(interface)]](a, b [[, loc]] [[, array_size]])` + +### Class + +Impure elemental function + +### Arguments + +`a`: has `intent(in)` and is a scalar of type `real` or `complex`. +If `a` is `real`, its value must be positive. If `a` is `complex`, both the real and imaginary components must be positive. This is the first shape parameter of the distribution. + +`b`: has `intent(in)` and is a scalar of type `real` or `complex`. +If `b` is `real`, its value must be positive. If `b` is `complex`, both the real and imaginary components must be positive. This is the second shape parameter of the distribution. + +`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. +Specifies the location (shift) of the distribution with default value 0.0. The distribution support is loc < x < loc + 1. + +`array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind. + +### Return value + +The result is a scalar or rank-1 array with a size of `array_size`, and the same type as `a`. If `a` or `b` is non-positive, the result is `NaN`. + +### Example + +```fortran +{!example/stats_distribution_beta/example_beta_rvs.f90!} +``` + +## `pdf_beta` - beta distribution probability density function + +### Status + +Experimental + +### Description + +The probability density function (pdf) of the single real variable beta distribution is: + +$$ f(x)= \\frac{x^{a-1}(1-x)^{b-1}}{B(a,b)} ,\\quad 00,\\ b>0 $$ + +where \\(a\\) and \\(b\\) are the shape parameters, and \\(B(a,b)\\) is the beta function. + +An optional `loc` parameter specifies the location (shift) of the distribution. + +For a complex variable \\(z=(x + y i)\\) with independent real \\(x\\) and imaginary \\(y\\) parts, the joint probability density function is the product of the corresponding real and imaginary marginal pdfs:[^2] + +$$f(x+\\mathit{i}y)=f(x)f(y)$$ + +### Syntax + +`result = [[stdlib_stats_distribution_beta(module):pdf_beta(interface)]](x, a, b [[, loc]])` + +### Class + +Impure elemental function + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. The point at which to evaluate the pdf. + +`a`: has `intent(in)` and is a scalar of type `real` or `complex`. The first shape parameter. +If `a` is `real`, its value must be positive. If `a` is `complex`, both the real and imaginary components must be positive. + +`b`: has `intent(in)` and is a scalar of type `real` or `complex`. The second shape parameter. +If `b` is `real`, its value must be positive. If `b` is `complex`, both the real and imaginary components must be positive. + +`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. The location (shift) parameter with default value 0.0. + +All arguments must have the same type. + +### Return value + +The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `a` or `b` is non-positive, the result is `NaN`. + +### Example + +```fortran +{!example/stats_distribution_beta/example_beta_pdf.f90!} +``` + +## `cdf_beta` - beta distribution cumulative distribution function + +### Status + +Experimental + +### Description + +Cumulative distribution function (cdf) of the single real variable beta distribution is: + +$$ F(x)= I_x(a, b),\\quad 00,\\ b>0 $$ + +where \\(I_x(a,b)\\) is the regularized incomplete beta function. + +An optional `loc` parameter specifies the location (shift) of the distribution. + +For a complex variable \\(z=(x + y i)\\) with independent real \\(x\\) and imaginary \\(y\\) parts, the joint cumulative distribution function is the product of the corresponding real and imaginary marginal cdfs:[^2] + +$$F(x+\\mathit{i}y)=F(x)F(y)$$ + +### Syntax + +`result = [[stdlib_stats_distribution_beta(module):cdf_beta(interface)]](x, a, b [[, loc]])` + +### Class + +Impure elemental function + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. The point at which to evaluate the cdf. + +`a`: has `intent(in)` and is a scalar of type `real` or `complex`. The first shape parameter. +If `a` is `real`, its value must be positive. If `a` is `complex`, both the real and imaginary components must be positive. + +`b`: has `intent(in)` and is a scalar of type `real` or `complex`. The second shape parameter. +If `b` is `real`, its value must be positive. If `b` is `complex`, both the real and imaginary components must be positive. + +`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. The location (shift) parameter with default value 0.0. + +All arguments must have the same type. + +### Return value + +The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `a` or `b` is non-positive, the result is `NaN`. + +### Example + +```fortran +{!example/stats_distribution_beta/example_beta_cdf.f90!} +``` + +[^1]: Devroye, Luc. _Non-Uniform Random Variate Generation_. Springer-Verlag, 1986 (Chapter IX, Section 3). + +[^2]: Miller, Scott, and Donald Childers. _Probability and random processes: With applications to signal processing and communications_. Academic Press, 2012 (p. 197). diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index efd1ae171..eac7c78b1 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -51,6 +51,7 @@ if (STDLIB_SPECIALMATRICES) endif() if (STDLIB_STATS) add_subdirectory(stats) + add_subdirectory(stats_distribution_beta) add_subdirectory(stats_distribution_exponential) add_subdirectory(stats_distribution_gamma) add_subdirectory(stats_distribution_normal) diff --git a/example/stats_distribution_beta/CMakeLists.txt b/example/stats_distribution_beta/CMakeLists.txt new file mode 100644 index 000000000..5e5b2663e --- /dev/null +++ b/example/stats_distribution_beta/CMakeLists.txt @@ -0,0 +1,4 @@ +ADD_EXAMPLE(beta_rvs) +ADD_EXAMPLE(beta_pdf) +ADD_EXAMPLE(beta_cdf) + diff --git a/example/stats_distribution_beta/example_beta_cdf.f90 b/example/stats_distribution_beta/example_beta_cdf.f90 new file mode 100644 index 000000000..19de5c4e2 --- /dev/null +++ b/example/stats_distribution_beta/example_beta_cdf.f90 @@ -0,0 +1,29 @@ +program example_beta_cdf + use stdlib_random, only: random_seed + use stdlib_stats_distribution_beta, only: rbeta => rvs_beta, & + beta_cdf => cdf_beta + + implicit none + real, parameter :: a = 2.0, b = 5.0 + real :: xarr(2, 5) + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + ! cumulative probability at x=0.3 for beta(2,5) distribution + print *, beta_cdf(0.3, a, b) + ! 0.579824865 + + ! cumulative probability at x=1.3 with loc=1.0 for beta(2,5) distribution + print *, beta_cdf(1.3, a, b, 1.0) + ! 0.579824746 + + ! generate random variates and compute their cdf + xarr = reshape(rbeta(a, b, 10), [2, 5]) + + print *, beta_cdf(xarr, a, b) + ! 0.686331749 0.625633657 0.625057578 0.158218294 0.786031485 + ! 7.17176571E-02 0.136123925 0.909627795 0.245356008 0.865481198 + +end program example_beta_cdf diff --git a/example/stats_distribution_beta/example_beta_pdf.f90 b/example/stats_distribution_beta/example_beta_pdf.f90 new file mode 100644 index 000000000..e7c34f974 --- /dev/null +++ b/example/stats_distribution_beta/example_beta_pdf.f90 @@ -0,0 +1,29 @@ +program example_beta_pdf + use stdlib_random, only: random_seed + use stdlib_stats_distribution_beta, only: rbeta => rvs_beta, & + beta_pdf => pdf_beta + + implicit none + real, parameter :: a = 2.0, b = 5.0 + real :: xarr(2, 5) + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + ! probability density at x=0.3 for beta(2,5) distribution + print *, beta_pdf(0.3, a, b) + ! 2.16089988 + + ! probability density at x=1.3 with loc=1.0 for beta(2,5) distribution + print *, beta_pdf(1.3, a, b, 1.0) + ! 2.16090035 + + ! generate random variates and compute their pdf + xarr = reshape(rbeta(a, b, 10), [2, 5]) + + print *, beta_pdf(xarr, a, b) + ! 1.85633695 2.04246974 2.04407597 2.16859674 1.47261345 + ! 1.67244565 2.07803488 0.819388986 2.38697886 1.08206940 + +end program example_beta_pdf diff --git a/example/stats_distribution_beta/example_beta_rvs.f90 b/example/stats_distribution_beta/example_beta_rvs.f90 new file mode 100644 index 000000000..ce739c29e --- /dev/null +++ b/example/stats_distribution_beta/example_beta_rvs.f90 @@ -0,0 +1,47 @@ +program example_beta_rvs + use stdlib_random, only: random_seed + use stdlib_stats_distribution_beta, only: rbeta => rvs_beta + + implicit none + real :: a_arr(2, 3, 4) + complex :: ca, cb + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + ! single beta random variate with a=2.0, b=5.0 (loc=0.0 by default) + print *, rbeta(2.0, 5.0) + ! 0.235164985 + + ! beta random variate with a=2.0, b=5.0, loc=1.0 + print *, rbeta(2.0, 5.0, 1.0) + ! 1.23516498 + + ! a rank-3 array of 24 beta random variates with a=0.5, b=0.5 + a_arr(:, :, :) = 0.5 + print *, rbeta(a_arr, a_arr) + ! 0.894186497 0.948506236 0.899142742 0.293822825 0.751733482 + ! 0.170928627 0.742042720 0.921871543 0.112629898 0.153393656 + ! 0.188625366 0.291826040 0.238829076 0.764039755 0.935611486 + ! 0.454867721 8.74810152E-03 0.258653969 0.963788986 + ! 0.202841997 0.689699173 0.537226677 0.721585333 0.891451001 + + ! an array of 10 random variates with a=2.0, b=5.0 (loc=0.0 by default) + print *, rbeta(2.0, 5.0, 10) + ! 2.59639323E-02 0.401881814 0.451093256 0.863215625 6.78956718E-03 + ! 0.316774905 0.141516894 0.199765816 0.616839588 0.555854380 + + ! an array of 10 random variates with a=2.0, b=5.0, loc=1.0 + print *, rbeta(2.0, 5.0, 10, 1.0) + ! 1.02596393 1.40188181 1.45109326 1.86321562 1.00678957 + ! 1.31677490 1.14151689 1.19976582 1.61683959 1.55585438 + + ca = (2.0, 3.0) + cb = (5.0, 4.0) + ! single complex beta random variate with real part a=2.0, b=5.0; + ! imaginary part a=3.0, b=4.0 (loc=(0,0) by default) + print *, rbeta(ca, cb) + ! (0.247691274,0.337867618) + +end program example_beta_rvs diff --git a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp index 988a3e73e..d60f6ae07 100644 --- a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp +++ b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp @@ -25,6 +25,7 @@ module stdlib_specialfunctions_gamma public :: lower_incomplete_gamma, log_lower_incomplete_gamma public :: upper_incomplete_gamma, log_upper_incomplete_gamma public :: regularized_gamma_p, regularized_gamma_q + public :: beta, log_beta, incomplete_beta @@ -184,6 +185,34 @@ module stdlib_specialfunctions_gamma + interface beta + !! Beta function + !! + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + module procedure beta_${t1[0]}$${k1}$ + #:endfor + end interface beta + + + + interface log_beta + !! Logarithm of beta function + !! + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + module procedure log_beta_${t1[0]}$${k1}$ + #:endfor + end interface log_beta + + + + interface incomplete_beta + !! Regularized incomplete beta function I_x(a,b) + !! + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + module procedure incomplete_beta_${t1[0]}$${k1}$ + #:endfor + end interface incomplete_beta + contains @@ -1258,4 +1287,145 @@ contains #:endfor #:endfor + + + ! + ! Beta function implementations + ! + + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + elemental function beta_${t1[0]}$${k1}$(a, b) result(res) + ! Beta function: beta(a,b) = gamma(a)*gamma(b)/gamma(a+b) + ! + ${t1}$, intent(in) :: a, b + ${t1}$ :: res + + if(a <= 0.0_${k1}$ .or. b <= 0.0_${k1}$) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + else + ! Use log-gamma for numerical stability + res = exp(log_gamma(a) + log_gamma(b) - log_gamma(a + b)) + end if + end function beta_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + elemental function log_beta_${t1[0]}$${k1}$(a, b) result(res) + ! Logarithm of beta function + ! + ${t1}$, intent(in) :: a, b + ${t1}$ :: res + + if(a <= 0.0_${k1}$ .or. b <= 0.0_${k1}$) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + else + res = log_gamma(a) + log_gamma(b) - log_gamma(a + b) + end if + end function log_beta_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + elemental function incomplete_beta_${t1[0]}$${k1}$(x, a, b) result(res) + ! Regularized incomplete beta function I_x(a,b) + ! Uses continued fraction expansion for numerical accuracy + ! + ${t1}$, intent(in) :: x, a, b + ${t1}$ :: res + ${t1}$ :: bt, xc + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + ${t1}$, parameter :: eps = epsilon(1.0_${k1}$) + integer, parameter :: maxit = 200 + + if((a <= zero .or. b <= zero) .or. (x < zero .or. x > one)) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + return + end if + + if(x == zero .or. x == one) then + res = x + return + end if + + ! Compute the logarithm of the beta function coefficient + bt = exp(log_gamma(a + b) - log_gamma(a) - log_gamma(b) & + + a * log(x) + b * log(one - x)) + + ! Use symmetry relation if necessary for better convergence + xc = x + if(x < (a + one) / (a + b + 2.0_${k1}$)) then + res = bt * betacf_${t1[0]}$${k1}$(x, a, b, eps, maxit) / a + else + ! Use symmetry relation I_x(a,b) = 1 - I_(1-x)(b,a) + xc = one - x + bt = exp(log_gamma(a + b) - log_gamma(a) - log_gamma(b) & + + b * log(xc) + a * log(one - xc)) + res = one - bt * betacf_${t1[0]}$${k1}$(xc, b, a, eps, maxit) / b + end if + end function incomplete_beta_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + pure function betacf_${t1[0]}$${k1}$(x, a, b, eps, maxit) result(res) + ! Continued fraction for incomplete beta function + ! Internal use only + ! + ${t1}$, intent(in) :: x, a, b, eps + integer, intent(in) :: maxit + ${t1}$ :: res + ${t1}$ :: aa, c, d, del, h, qab, qam, qap, fpmin + integer :: m, m2 + logical :: converged + ${t1}$, parameter :: one = 1.0_${k1}$ + + fpmin = tiny(1.0_${k1}$) / eps + + qab = a + b + qap = a + one + qam = a - one + c = one + d = one - qab * x / qap + if(abs(d) < fpmin) d = fpmin + d = one / d + h = d + converged = .false. + + do m = 1, maxit + m2 = 2 * m + aa = m * (b - m) * x / ((qam + m2) * (a + m2)) + d = one + aa * d + if(abs(d) < fpmin) d = fpmin + c = one + aa / c + if(abs(c) < fpmin) c = fpmin + d = one / d + h = h * d * c + aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)) + d = one + aa * d + if(abs(d) < fpmin) d = fpmin + c = one + aa / c + if(abs(c) < fpmin) c = fpmin + d = one / d + del = d * c + h = h * del + if(abs(del - one) < eps) then + converged = .true. + exit + end if + end do + + if (converged) then + res = h + else + res = ieee_value(one, ieee_quiet_nan) + end if + end function betacf_${t1[0]}$${k1}$ + + #:endfor + + end module stdlib_specialfunctions_gamma diff --git a/src/stats/CMakeLists.txt b/src/stats/CMakeLists.txt index 548138ee9..2c7be18b1 100644 --- a/src/stats/CMakeLists.txt +++ b/src/stats/CMakeLists.txt @@ -5,6 +5,7 @@ set(stats_fppFiles stdlib_random.fypp stdlib_stats_corr.fypp stdlib_stats_cov.fypp + stdlib_stats_distribution_beta.fypp stdlib_stats_distribution_exponential.fypp stdlib_stats_distribution_gamma.fypp stdlib_stats_distribution_normal.fypp diff --git a/src/stats/stdlib_stats_distribution_beta.fypp b/src/stats/stdlib_stats_distribution_beta.fypp new file mode 100644 index 000000000..0ad5bac43 --- /dev/null +++ b/src/stats/stdlib_stats_distribution_beta.fypp @@ -0,0 +1,370 @@ +#:include "common.fypp" +#:set BETA_REAL_KINDS_TYPES = [item for item in REAL_KINDS_TYPES if item[0] in ('sp', 'dp', 'xdp')] +#:set BETA_CMPLX_KINDS_TYPES = [item for item in CMPLX_KINDS_TYPES if item[0] in ('sp', 'dp', 'xdp')] +#:set RC_KINDS_TYPES = BETA_REAL_KINDS_TYPES + BETA_CMPLX_KINDS_TYPES +module stdlib_stats_distribution_beta + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_nan + use stdlib_kinds, only : sp, dp, xdp + use stdlib_error, only : error_stop + use stdlib_optval, only : optval + use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform + use stdlib_stats_distribution_gamma, only : rgamma=>rvs_gamma + use stdlib_specialfunctions_gamma, only : beta, incomplete_beta, log_beta + + implicit none + private + + public :: rvs_beta + public :: pdf_beta + public :: cdf_beta + + + interface rvs_beta + !! Version experimental + !! + !! Beta Distribution Random Variates + !! ([Specification](../page/specs/stdlib_stats_distribution_beta.html# + !! rvs_beta-beta-distribution-random-variates)) + !! + #: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 rvs_beta + + + interface pdf_beta + !! Version experimental + !! + !! Beta Distribution Probability Density Function + !! ([Specification](../page/specs/stdlib_stats_distribution_beta.html# + !! pdf_beta-beta-distribution-probability-density-function)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_pdf_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_pdf_impure_${t1[0]}$${k1}$ + #:endfor + end interface pdf_beta + + + interface cdf_beta + !! Version experimental + !! + !! Beta Distribution Cumulative Distribution Function + !! ([Specification](../page/specs/stdlib_stats_distribution_beta.html# + !! cdf_beta-beta-distribution-cumulative-distribution-function)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_cdf_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure beta_dist_cdf_impure_${t1[0]}$${k1}$ + #:endfor + end interface cdf_beta + + + + +contains + + #:for k1, t1 in BETA_REAL_KINDS_TYPES + impure elemental function beta_dist_rvs_${t1[0]}$${k1}$(a, b, loc) & + result(res) + ! Beta random variate using gamma variates + ! If a < 1 or b < 1, uses uniform method, otherwise uses gamma method + ! + ${t1}$, intent(in) :: a, b + ${t1}$, intent(in), optional :: loc + ${t1}$ :: res, x, y, xx(2) + ${t1}$ :: loc_ + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + loc_ = optval(loc, 0.0_${k1}$) + + if(a <= zero .or. b <= zero) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + return + end if + + if(a < one .or. b < one) then + ! Use uniform method for small parameters + do + xx = uni(zero, one, 2) + x = xx(1) ** (one / a) + y = xx(2) ** (one / b) + y = x + y + if(y <= one .and. y /= zero) exit + end do + else + ! Use gamma method for larger parameters + do + x = rgamma(a, one) + y = rgamma(b, one) + y = x + y + if(y /= zero) exit + end do + end if + + res = x / y + loc_ + end function beta_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_CMPLX_KINDS_TYPES + impure elemental function beta_dist_rvs_${t1[0]}$${k1}$(a, b, loc) & + result(res) + ! Complex parameter beta distributed. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: a, b + ${t1}$, intent(in), optional :: loc + ${t1}$ :: res + ${t1}$ :: loc_ + + loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$)) + + res = cmplx(beta_dist_rvs_r${k1}$(a%re, b%re) + loc_%re, & + beta_dist_rvs_r${k1}$(a%im, b%im) + loc_%im, kind=${k1}$) + end function beta_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_REAL_KINDS_TYPES + function beta_dist_rvs_array_${t1[0]}$${k1}$(a, b, array_size, loc) & + result(res) + ! + ${t1}$, intent(in) :: a, b + integer, intent(in) :: array_size + ${t1}$, intent(in), optional :: loc + ${t1}$ :: res(array_size) + integer :: i + ${t1}$ :: loc_ + + loc_ = optval(loc, 0.0_${k1}$) + + do i = 1, array_size + res(i) = beta_dist_rvs_${t1[0]}$${k1}$(a, b, loc_) + end do + end function beta_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_CMPLX_KINDS_TYPES + function beta_dist_rvs_array_${t1[0]}$${k1}$(a, b, array_size, loc) & + result(res) + ! Complex parameter beta distributed. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: a, b + integer, intent(in) :: array_size + ${t1}$, intent(in), optional :: loc + ${t1}$ :: res(array_size) + integer :: i + ${t1}$ :: loc_ + + loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$)) + + do i = 1, array_size + res(i) = beta_dist_rvs_${t1[0]}$${k1}$(a, b, loc_) + end do + end function beta_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_REAL_KINDS_TYPES + elemental function beta_dist_pdf_${t1[0]}$${k1}$(x, a, b, loc) & + result(res) + ! Beta distribution probability density function + ! + ${t1}$, intent(in) :: x, a, b + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + ${t1}$ :: xs + ${t1}$ :: loc_ + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + if(a <= zero .or. b <= zero) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + return + end if + + loc_ = optval(loc, 0.0_${k1}$) + xs = x - loc_ + + if(xs <= zero .or. xs >= one) then + res = zero + return + end if + + ! Use log formulation for numerical stability + res = exp((a - one) * log(xs) + (b - one) * log(one - xs) - log_beta(a, b)) + end function beta_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_REAL_KINDS_TYPES + impure elemental function beta_dist_pdf_impure_${t1[0]}$${k1}$(x, a, b, err, loc) & + result(res) + ! Beta distribution probability density function (impure wrapper) + ! + ${t1}$, intent(in) :: x, a, b + integer, intent(out) :: err + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + + res = beta_dist_pdf_${t1[0]}$${k1}$(x, a, b, loc) + err = merge(1, 0, ieee_is_nan(res)) + end function beta_dist_pdf_impure_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_CMPLX_KINDS_TYPES + elemental function beta_dist_pdf_${t1[0]}$${k1}$(x, a, b, loc) & + result(res) + ! Complex parameter beta distributed. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: x, a, b + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + ${t1}$ :: loc_ + + loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$)) + + res = beta_dist_pdf_r${k1}$(x%re, a%re, b%re, loc_%re) + res = res * beta_dist_pdf_r${k1}$(x%im, a%im, b%im, loc_%im) + end function beta_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_CMPLX_KINDS_TYPES + impure elemental function beta_dist_pdf_impure_${t1[0]}$${k1}$(x, a, b, err, loc) & + result(res) + ! Complex parameter beta distributed. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: x, a, b + integer, intent(out) :: err + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + ${t1}$ :: loc_ + + loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$)) + + res = beta_dist_pdf_r${k1}$(x%re, a%re, b%re, loc_%re) + res = res * beta_dist_pdf_r${k1}$(x%im, a%im, b%im, loc_%im) + err = merge(1, 0, ieee_is_nan(res)) + end function beta_dist_pdf_impure_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_REAL_KINDS_TYPES + elemental function beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc) & + result(res) + ! Beta distribution cumulative distribution function + ! + ${t1}$, intent(in) :: x, a, b + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + ${t1}$ :: xs + ${t1}$ :: loc_ + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + loc_ = optval(loc, 0.0_${k1}$) + xs = x - loc_ + + if(a <= zero .or. b <= zero) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + return + end if + + if(xs <= zero) then + res = zero + return + end if + + if(xs >= one) then + res = one + return + end if + + res = real(incomplete_beta(xs, a, b), kind=${k1}$) + end function beta_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_REAL_KINDS_TYPES + impure elemental function beta_dist_cdf_impure_${t1[0]}$${k1}$(x, a, b, err, loc) & + result(res) + ! Beta distribution cumulative distribution function (impure wrapper) + ! + ${t1}$, intent(in) :: x, a, b + integer, intent(out) :: err + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + + res = beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc) + err = merge(1, 0, ieee_is_nan(res)) + end function beta_dist_cdf_impure_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_CMPLX_KINDS_TYPES + elemental function beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc) & + result(res) + ! Complex parameter beta distributed. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: x, a, b + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + ${t1}$ :: loc_ + + loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$)) + + res = beta_dist_cdf_r${k1}$(x%re, a%re, b%re, loc_%re) + res = res * beta_dist_cdf_r${k1}$(x%im, a%im, b%im, loc_%im) + end function beta_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_CMPLX_KINDS_TYPES + impure elemental function beta_dist_cdf_impure_${t1[0]}$${k1}$(x, a, b, err, loc) & + result(res) + ! Complex parameter beta distributed. The real part and imaginary part are + ! independent of each other. + ! + ${t1}$, intent(in) :: x, a, b + integer, intent(out) :: err + ${t1}$, intent(in), optional :: loc + real(${k1}$) :: res + ${t1}$ :: loc_ + + loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$)) + + res = beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc_) + err = merge(1, 0, ieee_is_nan(res)) + end function beta_dist_cdf_impure_${t1[0]}$${k1}$ + + #:endfor + + +end module stdlib_stats_distribution_beta diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index fe42aa3c8..96979a9a4 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -12,7 +12,8 @@ module test_specialfunctions_gamma log_lower_incomplete_gamma, & log_upper_incomplete_gamma, & regularized_gamma_p, & - regularized_gamma_q + regularized_gamma_q, & + beta, log_beta, incomplete_beta implicit none private @@ -44,7 +45,7 @@ contains #:endfor #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in GAMMA_REAL_KINDS_TYPES + #:for k2, t2 in GAMMA_REAL_KINDS_TYPES , new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$", & test_lincgamma_${t1[0]}$${k1}$${k2}$) & , new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$", & @@ -73,6 +74,12 @@ contains test_gamma_p_${t1[0]}$${k1}$) & , new_unittest("regularized_gamma_q_${t1[0]}$${k1}$", & test_gamma_q_${t1[0]}$${k1}$) & + , new_unittest("beta_${t1[0]}$${k1}$", & + test_beta_${t1[0]}$${k1}$) & + , new_unittest("log_beta_${t1[0]}$${k1}$", & + test_log_beta_${t1[0]}$${k1}$) & + , new_unittest("incomplete_beta_${t1[0]}$${k1}$", & + test_incomplete_beta_${t1[0]}$${k1}$) & #:endfor ] end subroutine collect_specialfunctions_gamma @@ -397,12 +404,83 @@ contains end do end subroutine test_gamma_q_${t1[0]}$${k1}$${k2}$ - #:endfor - #:endfor + #:endfor + #:endfor + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES - #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + + subroutine test_beta_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 3 + integer :: i + ${t1}$, parameter :: a(n) = [2.0_${k1}$, 1.0_${k1}$, 0.5_${k1}$] + ${t1}$, parameter :: b(n) = [3.0_${k1}$, 1.0_${k1}$, 0.5_${k1}$] + + ${t1}$, parameter :: ans(n) = [1.0_${k1}$/12.0_${k1}$, & + 1.0_${k1}$, & + acos(-1.0_${k1}$)] + + do i = 1, n + + call check(error, beta(a(i), b(i)), ans(i), & + "Beta function with a(kind=${k1}$) and b(kind=${k1}$) failed", & + thr = tol_${k1}$, rel = .true.) + if (allocated(error)) return + + end do + end subroutine test_beta_${t1[0]}$${k1}$ + + + + subroutine test_log_beta_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 3 + integer :: i + character(len=120) :: msg + ${t1}$ :: a(n) = [2.0_${k1}$, 1.0_${k1}$, 0.5_${k1}$] + ${t1}$ :: b(n) = [3.0_${k1}$, 1.0_${k1}$, 0.5_${k1}$] + + ${t1}$, parameter :: ans(n) = [log(1.0_${k1}$/12.0_${k1}$), & + 0.0_${k1}$, & + log(acos(-1.0_${k1}$))] + + do i = 1, n + write(msg, '(a,i0)') "Log-beta function with a(kind=${k1}$) and b(kind=${k1}$) failed at i=", i + + call check(error, log_beta(a(i), b(i)), ans(i), & + msg, & + thr = tol_${k1}$, rel = .true.) + if (allocated(error)) return + + end do + end subroutine test_log_beta_${t1[0]}$${k1}$ + + + + subroutine test_incomplete_beta_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: x(n) = [0.25_${k1}$, 0.5_${k1}$, 0.0_${k1}$, 1.0_${k1}$] + ${t1}$ :: a(n) = [2.0_${k1}$, 2.0_${k1}$, 2.0_${k1}$, 2.0_${k1}$] + ${t1}$ :: b(n) = [3.0_${k1}$, 3.0_${k1}$, 3.0_${k1}$, 3.0_${k1}$] + + ${t1}$, parameter :: ans(n) = [67.0_${k1}$/256.0_${k1}$, & + 11.0_${k1}$/16.0_${k1}$, & + 0.0_${k1}$, & + 1.0_${k1}$] + + do i = 1, n + + call check(error, incomplete_beta(x(i), a(i), b(i)), ans(i), & + "Incomplete beta I_x(a,b) with kind=${k1}$ failed", & + thr = tol_${k1}$, rel = .true.) + if (allocated(error)) return + + end do + end subroutine test_incomplete_beta_${t1[0]}$${k1}$ subroutine test_lincgamma_${t1[0]}$${k1}$(error) type(error_type), allocatable, intent(out) :: error diff --git a/test/stats/CMakeLists.txt b/test/stats/CMakeLists.txt index 94990ccf1..e38134678 100644 --- a/test/stats/CMakeLists.txt +++ b/test/stats/CMakeLists.txt @@ -8,6 +8,7 @@ set(fppFiles test_distribution_uniform.fypp test_distribution_normal.fypp test_distribution_exponential.fypp + test_distribution_beta.fypp test_distribution_gamma.fypp ) @@ -25,6 +26,7 @@ ADDTEST(random) ADDTEST(distribution_uniform) ADDTEST(distribution_normal) ADDTEST(distribution_exponential) +ADDTEST(distribution_beta) ADDTEST(distribution_gamma) if(DEFINED CMAKE_MAXIMUM_RANK) diff --git a/test/stats/test_distribution_beta.fypp b/test/stats/test_distribution_beta.fypp new file mode 100644 index 000000000..b33532a91 --- /dev/null +++ b/test/stats/test_distribution_beta.fypp @@ -0,0 +1,256 @@ +#:include "common.fypp" +#:set BETA_REAL_KINDS_TYPES = [item for item in REAL_KINDS_TYPES if item[0] in ('sp', 'dp', 'xdp')] +#:set BETA_CMPLX_KINDS_TYPES = [item for item in CMPLX_KINDS_TYPES if item[0] in ('sp', 'dp', 'xdp')] +#:set RC_KINDS_TYPES = BETA_REAL_KINDS_TYPES + BETA_CMPLX_KINDS_TYPES + +program test_distribution_beta + use stdlib_kinds, only : sp, dp, xdp + use stdlib_error, only : check + use stdlib_random, only : random_seed + use stdlib_stats_distribution_beta, only : rbeta => rvs_beta, & + beta_pdf => pdf_beta, beta_cdf => cdf_beta + + implicit none + #:for k1, t1 in BETA_REAL_KINDS_TYPES + ${t1}$, parameter :: ${k1}$tol = 1000 * epsilon(1.0_${k1}$) + #:endfor + logical :: warn = .true. + integer :: put, get + + put = 12345678 + call random_seed(put, get) + + call test_beta_random_generator + + #:for k1, t1 in RC_KINDS_TYPES + call test_beta_rvs_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + call test_beta_pdf_${t1[0]}$${k1}$ + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + call test_beta_cdf_${t1[0]}$${k1}$ + #:endfor + +contains + + subroutine test_beta_random_generator + integer, parameter :: num = 10000000, array_size = 1000 + integer :: i, j, freq(0:array_size-1) + real(dp) :: chisq, expct + + print *, "" + print *, "Test beta random generator with chi-squared" + + freq = 0 + do i = 1, num + j = min(array_size - 1, int(array_size * beta_cdf(rbeta(2.0_dp, 5.0_dp, 0.0_dp), 2.0_dp, 5.0_dp, 0.0_dp))) + freq(j) = freq(j) + 1 + end do + + chisq = 0.0_dp + expct = num / array_size + + do i = 0, array_size - 1 + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + + write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & + //" 1143.92" + write(*,*) "Chi-squared for 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 + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_beta_rvs_${t1[0]}$${k1}$ + integer, parameter :: k = 5, n = 10 + ${t1}$ :: res(n), ans(n), ba, bb, loc + integer :: i + integer :: seed, get + + print *, "Test beta_distribution_rvs_${t1[0]}$${k1}$" + seed = 639741825 + call random_seed(seed, get) + + #:if t1[0] == "r" + ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + #:else + ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) + #:endif + + do i = 1, k + res(i) = rbeta(ba, bb, loc) + end do + res(k + 1 : n) = rbeta(ba, bb, k, loc) + + seed = 639741825 + call random_seed(seed, get) + do i = 1, n + ans(i) = rbeta(ba, bb, loc) + end do + + do i = 1, n + call check(abs(res(i) - ans(i)) < ${k1}$tol, & + msg="beta_distribution_rvs_${t1[0]}$${k1}$ failed", & + warn=warn) + end do + end subroutine test_beta_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_beta_pdf_${t1[0]}$${k1}$ + ${t1}$ :: x1, x2(3,4), ba, bb, loc + integer, parameter :: n = 15 + integer :: i + real(${k1}$) :: res(n) + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(*) = & + [2.45759999999999978E+00_${k1}$, & + 2.45759999999999978E+00_${k1}$, & + 2.45759999999999978E+00_${k1}$, & + 1.22175937500000087E+00_${k1}$, & + 1.96829999999999972E+00_${k1}$, & + 2.34902812499999891E+00_${k1}$, & + 2.37304687500000044E+00_${k1}$, & + 2.16089999999999938E+00_${k1}$, & + 1.87431562500000037E+00_${k1}$, & + 1.55520000000000058E+00_${k1}$, & + 1.23533437500000032E+00_${k1}$, & + 9.37499999999999889E-01_${k1}$, & + 4.60800000000000043E-01_${k1}$, & + 1.70099999999999946E-01_${k1}$, & + 3.83999999999999897E-02_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(*) = & + [2.77620563793055197E+00_${k1}$, & + 2.77620563793055197E+00_${k1}$, & + 2.77620563793055197E+00_${k1}$, & + 4.35133599581893460E+00_${k1}$, & + 5.11042526155962751E+00_${k1}$, & + 3.91417214341107966E+00_${k1}$, & + 3.25033085838892477E+00_${k1}$, & + 2.44104116333175813E+00_${k1}$, & + 1.74312991236416215E+00_${k1}$, & + 1.18399939600181159E+00_${k1}$, & + 7.62828448855675023E-01_${k1}$, & + 4.63554726571548781E-01_${k1}$, & + 1.79353320295729979E-01_${k1}$, & + 5.09635475047793968E-02_${k1}$, & + 6.17909755246391895E-03_${k1}$] + #:endif + + print *, "Test beta_distribution_pdf_${t1[0]}$${k1}$" + #:if t1[0] == "r" + ba = 2.0_${k1}$ + bb = 5.0_${k1}$ + loc = 0._${k1}$ + x1 = 0.2_${k1}$ + x2 = reshape([0.05_${k1}$, 0.1_${k1}$, 0.15_${k1}$, 0.25_${k1}$, & + 0.3_${k1}$, 0.35_${k1}$, 0.4_${k1}$, 0.45_${k1}$, & + 0.5_${k1}$, 0.6_${k1}$, 0.7_${k1}$, 0.8_${k1}$], [3,4]) + #:else + ba = (2.0_${k1}$, 0.7_${k1}$) + bb = (5.0_${k1}$, 3.0_${k1}$) + loc = (0._${k1}$, 0._${k1}$) + x1 = (0.2_${k1}$, 0.3_${k1}$) + x2 = reshape([(0.05_${k1}$, 0.05_${k1}$), (0.1_${k1}$, 0.1_${k1}$), & + (0.15_${k1}$, 0.2_${k1}$), (0.25_${k1}$, 0.25_${k1}$), & + (0.3_${k1}$, 0.3_${k1}$), (0.35_${k1}$, 0.35_${k1}$), & + (0.4_${k1}$, 0.4_${k1}$), (0.45_${k1}$, 0.45_${k1}$), & + (0.5_${k1}$, 0.5_${k1}$), (0.6_${k1}$, 0.55_${k1}$), & + (0.7_${k1}$, 0.6_${k1}$), (0.8_${k1}$, 0.7_${k1}$)], [3,4]) + #:endif + + res(1:3) = beta_pdf(x1, ba, bb, loc) + res(4:15) = reshape(beta_pdf(x2, ba, bb, loc), [12]) + + do i = 1, n + call check(abs(res(i) - ans(i)) < max(${k1}$tol, 1.0e-11_${k1}$), & + msg="beta_distribution_pdf_${t1[0]}$${k1}$ failed", & + warn=warn) + end do + end subroutine test_beta_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_beta_cdf_${t1[0]}$${k1}$ + ${t1}$ :: x1, x2(3,4), ba, bb, loc + integer :: i + real(${k1}$) :: res(15) + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [3.44640000000000224E-01_${k1}$, & + 3.44640000000000224E-01_${k1}$, & + 3.44640000000000224E-01_${k1}$, & + 3.27738281249999944E-02_${k1}$, & + 1.14265000000000019E-01_${k1}$, & + 2.23515703125000020E-01_${k1}$, & + 4.66064453125000000E-01_${k1}$, & + 5.79825000000000257E-01_${k1}$, & + 6.80920078124999995E-01_${k1}$, & + 7.66720000000000068E-01_${k1}$, & + 8.36432578124999937E-01_${k1}$, & + 8.90625000000000000E-01_${k1}$, & + 9.59040000000000004E-01_${k1}$, & + 9.89064999999999972E-01_${k1}$, & + 9.98399999999999954E-01_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [2.64331290012673414E-01_${k1}$, & + 2.64331290012673414E-01_${k1}$, & + 2.64331290012673414E-01_${k1}$, & + 8.86382195944738702E-03_${k1}$, & + 4.81500626030448228E-02_${k1}$, & + 1.40607931860570912E-01_${k1}$, & + 3.28430860699885308E-01_${k1}$, & + 4.44713005546652551E-01_${k1}$, & + 5.57213246702040532E-01_${k1}$, & + 6.59756977764087704E-01_${k1}$, & + 7.48497876310954546E-01_${k1}$, & + 8.21680689855777691E-01_${k1}$, & + 9.05920357339824456E-01_${k1}$, & + 9.51254168091690833E-01_${k1}$, & + 9.82801096697218157E-01_${k1}$] + #:endif + + print *, "Test beta_distribution_cdf_${t1[0]}$${k1}$" + #:if t1[0] == "r" + ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + x1 = 0.2_${k1}$ + x2 = reshape([0.05_${k1}$, 0.1_${k1}$, 0.15_${k1}$, 0.25_${k1}$, & + 0.3_${k1}$, 0.35_${k1}$, 0.4_${k1}$, 0.45_${k1}$, & + 0.5_${k1}$, 0.6_${k1}$, 0.7_${k1}$, 0.8_${k1}$], [3,4]) + #:else + ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) + x1 = (0.2_${k1}$, 0.3_${k1}$) + x2 = reshape([(0.05_${k1}$, 0.05_${k1}$), (0.1_${k1}$, 0.1_${k1}$), & + (0.15_${k1}$, 0.2_${k1}$), (0.25_${k1}$, 0.25_${k1}$), & + (0.3_${k1}$, 0.3_${k1}$), (0.35_${k1}$, 0.35_${k1}$), & + (0.4_${k1}$, 0.4_${k1}$), (0.45_${k1}$, 0.45_${k1}$), & + (0.5_${k1}$, 0.5_${k1}$), (0.6_${k1}$, 0.55_${k1}$), & + (0.7_${k1}$, 0.6_${k1}$), (0.8_${k1}$, 0.7_${k1}$)], [3,4]) + #:endif + + res(1:3) = beta_cdf(x1, ba, bb, loc) + res(4:15) = reshape(beta_cdf(x2, ba, bb, loc), [12]) + + do i = 1, 15 + call check(abs(res(i) - ans(i)) < max(${k1}$tol, 1.0e-11_${k1}$), & + msg="beta_distribution_cdf_${t1[0]}$${k1}$ failed", & + warn=warn) + end do + end subroutine test_beta_cdf_${t1[0]}$${k1}$ + + #:endfor + +end program test_distribution_beta