From 154090ff08cad5a83fb29714714514c16fe6b882 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:43:47 +0000 Subject: [PATCH 01/33] beta distribution documentation initial --- doc/specs/stdlib_stats_distribution_beta.md | 167 ++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 doc/specs/stdlib_stats_distribution_beta.md diff --git a/doc/specs/stdlib_stats_distribution_beta.md b/doc/specs/stdlib_stats_distribution_beta.md new file mode 100644 index 000000000..419853d2d --- /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. + +`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: + +$$ 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). From a8278e36307f94b3ccd5e23d683d56b5da4be589 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:46:14 +0000 Subject: [PATCH 02/33] beta subdirectory build --- example/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) 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) From b9c8c8f3f64dacd3e6d9c49f384baf2580554a2b Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:47:01 +0000 Subject: [PATCH 03/33] makes dir --- example/stats_distribution_beta/initial | 1 + 1 file changed, 1 insertion(+) create mode 100644 example/stats_distribution_beta/initial diff --git a/example/stats_distribution_beta/initial b/example/stats_distribution_beta/initial new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/example/stats_distribution_beta/initial @@ -0,0 +1 @@ + From 15144d16484e6fbadadc7f906ba2c216cd95a296 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:47:33 +0000 Subject: [PATCH 04/33] uploaded examples initial --- .../stats_distribution_beta/CMakeLists.txt | 4 ++ .../example_beta_cdf.f90 | 36 ++++++++++++++ .../example_beta_pdf.f90 | 36 ++++++++++++++ .../example_beta_rvs.f90 | 47 +++++++++++++++++++ 4 files changed, 123 insertions(+) create mode 100644 example/stats_distribution_beta/CMakeLists.txt create mode 100644 example/stats_distribution_beta/example_beta_cdf.f90 create mode 100644 example/stats_distribution_beta/example_beta_pdf.f90 create mode 100644 example/stats_distribution_beta/example_beta_rvs.f90 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..70e022e49 --- /dev/null +++ b/example/stats_distribution_beta/example_beta_cdf.f90 @@ -0,0 +1,36 @@ +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 :: x, a, b + real :: xarr(3, 4, 5) + integer :: seed_put, seed_get, i + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + a = 2.0; b = 5.0 + + ! cumulative probability at x=0.3 for beta(2,5) distribution + print *, beta_cdf(0.3, a, b) + ! 0.471808970 + + ! 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.471808970 + + ! generate random variates and compute their cdf + do i = 1, 60 + xarr(i, 1, 1) = rbeta(a, b) + end do + + print *, beta_cdf(xarr, a, b) + ! 0.374293357 0.136472717 0.153627276 0.166885555 0.535913110 + ! 4.47619371E-02 0.161991328 0.524897814 7.37934634E-02 8.41872990E-02 + ! 0.102836564 0.138294637 0.122145072 0.171486780 0.532901883 + ! 0.111447386 1.87034653E-02 0.119697235 0.550687969 9.98932496E-02 + ! 0.234867483 0.132088020 0.210879743 0.520989060 ... (60 elements total) + +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..8a15d8bf4 --- /dev/null +++ b/example/stats_distribution_beta/example_beta_pdf.f90 @@ -0,0 +1,36 @@ +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 :: x, a, b + real :: xarr(3, 4, 5) + integer :: seed_put, seed_get, i + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + a = 2.0; b = 5.0 + + ! probability density at x=0.3 for beta(2,5) distribution + print *, beta_pdf(0.3, a, b) + ! 1.32859993 + + ! probability density at x=1.3 with loc=1.0 for beta(2,5) distribution + print *, beta_pdf(1.3, a, b, 1.0) + ! 1.32859993 + + ! generate random variates and compute their pdf + do i = 1, 60 + xarr(i, 1, 1) = rbeta(a, b) + end do + + print *, beta_pdf(xarr, a, b) + ! 1.23914337 0.925268888 0.969913423 0.998693407 1.44889069 + ! 0.331506640 0.988773048 1.41048801 0.502468705 0.564825535 + ! 0.661424160 0.931534827 0.861649513 1.01826906 1.44227338 + ! 0.782453180 0.145308748 0.841068327 1.49033546 0.727149904 + ! 1.14145494 0.905850768 1.08809102 1.40226245 ... (60 elements total) + +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 From 4e7e2ab87e2079abe503ccc55fbbd6d044c54839 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:48:02 +0000 Subject: [PATCH 05/33] deleted initial file dir --- example/stats_distribution_beta/initial | 1 - 1 file changed, 1 deletion(-) delete mode 100644 example/stats_distribution_beta/initial diff --git a/example/stats_distribution_beta/initial b/example/stats_distribution_beta/initial deleted file mode 100644 index 8b1378917..000000000 --- a/example/stats_distribution_beta/initial +++ /dev/null @@ -1 +0,0 @@ - From e76aa94696bec4798ccc6ccf5baf58b0594f77c0 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:49:31 +0000 Subject: [PATCH 06/33] added beta to doc index --- doc/specs/index.md | 1 + 1 file changed, 1 insertion(+) 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 From 3efe58a2479e4c477c69a344e8df028d3eaebb2a Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:55:57 +0000 Subject: [PATCH 07/33] added beta functions to gamma special function more appropriate and efficient than an entirely new beta special function --- .../stdlib_specialfunctions_gamma.fypp | 170 ++++++++++++++++++ 1 file changed, 170 insertions(+) diff --git a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp index 9ee722b2d..e0f37b8ec 100644 --- a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp +++ b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp @@ -23,6 +23,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 @@ -182,6 +183,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 @@ -1256,4 +1285,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) + return + end if + + ! Use log-gamma for numerical stability + res = exp(log_gamma(a) + log_gamma(b) - log_gamma(a + b)) + 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) + return + end if + + res = log_gamma(a) + log_gamma(b) - log_gamma(a + b) + end function log_beta_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in GAMMA_REAL_KINDS_TYPES + impure 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) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + return + end if + + if(x < zero .or. x > one) then + res = ieee_value(1.0_${k1}$, ieee_quiet_nan) + return + end if + + if(x == zero) then + res = zero + return + end if + + if(x == one) then + res = one + 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 + integer :: m, m2 + ${t1}$, parameter :: one = 1.0_${k1}$ + ${t1}$, parameter :: 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 + + 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) exit + end do + + res = h + end function betacf_${t1[0]}$${k1}$ + + #:endfor + + end module stdlib_specialfunctions_gamma From 773566c881af0186daefc956e16426e348f5f2f5 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:57:29 +0000 Subject: [PATCH 08/33] added beta to build file --- src/stats/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) 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 From d6254367db7e41b7352206716e3693a2ecf94562 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:58:23 +0000 Subject: [PATCH 09/33] initial beta function --- src/stats/stdlib_stats_distribution_beta.fypp | 374 ++++++++++++++++++ 1 file changed, 374 insertions(+) create mode 100644 src/stats/stdlib_stats_distribution_beta.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..fd3eff4aa --- /dev/null +++ b/src/stats/stdlib_stats_distribution_beta.fypp @@ -0,0 +1,374 @@ +#: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 + + 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}$ + + 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 .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 = 0 + if(ieee_is_nan(res)) err = 1 + 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 = 0 + if(ieee_is_nan(res)) err = 1 + end function beta_dist_pdf_impure_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_REAL_KINDS_TYPES + impure 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 = 0 + if(ieee_is_nan(res)) err = 1 + end function beta_dist_cdf_impure_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in BETA_CMPLX_KINDS_TYPES + impure 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 = 0 + if(ieee_is_nan(res)) err = 1 + end function beta_dist_cdf_impure_${t1[0]}$${k1}$ + + #:endfor + + +end module stdlib_stats_distribution_beta From c7453059a21138f29c1f3ae527c60fc547467462 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:00:37 +0000 Subject: [PATCH 10/33] added beta to build file --- test/stats/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) 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) From 1619fe2baa9009d631c9a24a456bd2003189ed8d Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:01:06 +0000 Subject: [PATCH 11/33] added beta test --- test/stats/test_distribution_beta.fypp | 267 +++++++++++++++++++++++++ 1 file changed, 267 insertions(+) create mode 100644 test/stats/test_distribution_beta.fypp diff --git a/test/stats/test_distribution_beta.fypp b/test/stats/test_distribution_beta.fypp new file mode 100644 index 000000000..000322553 --- /dev/null +++ b/test/stats/test_distribution_beta.fypp @@ -0,0 +1,267 @@ +#: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), ba, bb, loc + integer :: i + integer :: seed, get + #:if t1[0] == "r" + #! for real type + ${t1}$ :: ans(n) = [0.23516498458385468E+00_${k1}$, & + 0.38906532525908948E+00_${k1}$, & + 0.34268730878829956E+00_${k1}$, & + 0.29513829946517944E+00_${k1}$, & + 0.19447886384116858E+00_${k1}$, & + 0.60383307933807373E+00_${k1}$, & + 0.58347225189208984E+00_${k1}$, & + 0.74638533592224121E+00_${k1}$, & + 0.55743497610092163E+00_${k1}$, & + 0.39149415493011475E+00_${k1}$] + #:else + #! for complex type + ${t1}$ :: ans(n) = & + [(0.44851309061050415E+00_${k1}$, 0.34970587491989136E+00_${k1}$), & + (0.25030922889709473E+00_${k1}$, 0.55935305356979370E+00_${k1}$), & + (0.69629210233688354E+00_${k1}$, 0.25098574161529541E+00_${k1}$), & + (0.54261916875839233E+00_${k1}$, 0.45337703824043274E+00_${k1}$), & + (0.25677326321601868E+00_${k1}$, 0.18765588104724884E+00_${k1}$), & + (0.59211272001266479E+00_${k1}$, 0.72924208641052246E+00_${k1}$), & + (0.75178939104080200E+00_${k1}$, 0.18360115587711334E+00_${k1}$), & + (0.41408637166023254E+00_${k1}$, 0.36145636439323425E+00_${k1}$), & + (0.47387075424194336E+00_${k1}$, 0.26313120126724243E+00_${k1}$), & + (0.71398687362670898E+00_${k1}$, 0.51318806409835815E+00_${k1}$)] + #:endif + + print *, "Test beta_distribution_rvs_${t1[0]}$${k1}$" + seed = 639741825 + call random_seed(seed, get) + + #:if t1[0] == "r" + #! for real type + ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + #:else + #! for complex type + 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) + + 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 :: i + integer :: seed, get + real(${k1}$) :: res(15) + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [1.23914337158203125E+00_${k1}$, & + 1.23914337158203125E+00_${k1}$, & + 1.23914337158203125E+00_${k1}$, & + 1.08365607261657715E+00_${k1}$, & + 1.04956150054931641E+00_${k1}$, & + 1.03561067581176758E+00_${k1}$, & + 1.22067928314208984E+00_${k1}$, & + 8.46962928771972656E-01_${k1}$, & + 2.52243280410766602E-01_${k1}$, & + 9.41652893066406250E-01_${k1}$, & + 6.57052993774414062E-01_${k1}$, & + 1.14095616340637207E+00_${k1}$, & + 9.09092426300048828E-01_${k1}$, & + 1.22890543937683105E+00_${k1}$, & + 4.40967559814453125E-01_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [1.15128993988037109E+00_${k1}$, & + 1.15128993988037109E+00_${k1}$, & + 1.15128993988037109E+00_${k1}$, & + 8.46223831176757812E-01_${k1}$, & + 3.30976486206054688E+00_${k1}$, & + 3.01892852783203125E+00_${k1}$, & + 1.23958110809326172E+00_${k1}$, & + 1.97296905517578125E+00_${k1}$, & + 1.17863655090332031E+00_${k1}$, & + 4.00447845458984375E-01_${k1}$, & + 2.92449951171875000E-02_${k1}$, & + 9.25743579864501953E-01_${k1}$, & + 2.15837764739990234E+00_${k1}$, & + 5.10101318359375000E-02_${k1}$, & + 1.03166580200195312E+00_${k1}$] + #:endif + + print *, "Test beta_distribution_pdf_${t1[0]}$${k1}$" + seed = 345987126 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + #:else + #! for complex type + ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) + #:endif + + x1 = rbeta(ba, bb, loc) + x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) + + res(1:3) = beta_pdf(x1, ba, bb, loc) + res(4:15) = reshape(beta_pdf(x2, ba, bb, loc), [12]) + + do i = 1, 15 + call check(abs(res(i) - ans(i)) < ${k1}$tol, & + 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 + integer :: seed, get + real(${k1}$) :: res(15) + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [0.37429335713386536E+00_${k1}$, & + 0.37429335713386536E+00_${k1}$, & + 0.37429335713386536E+00_${k1}$, & + 0.10568153858184814E+00_${k1}$, & + 0.61683785915374756E+00_${k1}$, & + 0.95208930969238281E-01_${k1}$, & + 0.34991249442100525E+00_${k1}$, & + 0.73912835121154785E+00_${k1}$, & + 0.95353639125823975E+00_${k1}$, & + 0.68856298923492432E+00_${k1}$, & + 0.36885082721710205E-01_${k1}$, & + 0.14866244792938232E+00_${k1}$, & + 0.74632394313812256E+00_${k1}$, & + 0.29785108566284180E+00_${k1}$, & + 0.14457702636718750E-01_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [0.14019477367401123E+00_${k1}$, & + 0.14019477367401123E+00_${k1}$, & + 0.14019477367401123E+00_${k1}$, & + 0.48782587051391602E+00_${k1}$, & + 0.19825363159179688E+00_${k1}$, & + 0.92399120330810547E-02_${k1}$, & + 0.55928134918212891E+00_${k1}$, & + 0.15336728096008301E+00_${k1}$, & + 0.31033802032470703E+00_${k1}$, & + 0.69235420227050781E+00_${k1}$, & + 0.92155790328979492E+00_${k1}$, & + 0.63153839111328125E+00_${k1}$, & + 0.86834430694580078E-01_${k1}$, & + 0.69353580474853516E+00_${k1}$, & + 0.59732055664062500E+00_${k1}$] + #:endif + + print *, "Test beta_distribution_cdf_${t1[0]}$${k1}$" + seed = 345987126 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + #:else + #! for complex type + ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) + #:endif + + x1 = rbeta(ba, bb, loc) + x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) + + 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)) < ${k1}$tol, & + 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 From 4320f5ed78bcedd2a035136779547024ede3fc05 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:04:32 +0000 Subject: [PATCH 12/33] added beta tests for special functions changes --- .../test_specialfunctions_gamma.fypp | 80 ++++++++++++++++++- 1 file changed, 78 insertions(+), 2 deletions(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 6249e51a3..6272c31df 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 @@ -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,7 +404,76 @@ contains end do end subroutine test_gamma_q_${t1[0]}$${k1}$${k2}$ - #:endfor + + + subroutine test_beta_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 3 + integer :: i + ${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) = [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.) + + 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 + ${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 + + call check(error, log_beta(a(i), b(i)), ans(i), & + "Log-beta function with a(kind=${k1}$) and b(kind=${k1}$) failed", & + thr = tol_${k1}$, rel = .true.) + + 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.) + + end do + end subroutine test_incomplete_beta_${t1[0]}$${k1}$ + + #:endfor +end module test_specialfunctions_gamma #:endfor From 00e93f3ac67391457155010a229e7b87de7f4644 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:28:38 +0000 Subject: [PATCH 13/33] use fpmin --- src/specialfunctions/stdlib_specialfunctions_gamma.fypp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp index e0f37b8ec..d39f926db 100644 --- a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp +++ b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp @@ -1386,10 +1386,11 @@ contains ${t1}$, intent(in) :: x, a, b, eps integer, intent(in) :: maxit ${t1}$ :: res - ${t1}$ :: aa, c, d, del, h, qab, qam, qap + ${t1}$ :: aa, c, d, del, h, qab, qam, qap, fpmin integer :: m, m2 ${t1}$, parameter :: one = 1.0_${k1}$ - ${t1}$, parameter :: fpmin = tiny(1.0_${k1}$) / eps + + fpmin = tiny(1.0_${k1}$) / eps qab = a + b qap = a + one From a0c24f5c127e88f5bb06aa5db50bae83e89dd655 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:45:21 +0000 Subject: [PATCH 14/33] lined up values to scipy output + refactor --- test/stats/test_distribution_beta.fypp | 268 ++++++++++--------------- 1 file changed, 106 insertions(+), 162 deletions(-) diff --git a/test/stats/test_distribution_beta.fypp b/test/stats/test_distribution_beta.fypp index 000322553..d6777c207 100644 --- a/test/stats/test_distribution_beta.fypp +++ b/test/stats/test_distribution_beta.fypp @@ -2,6 +2,7 @@ #: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 @@ -10,10 +11,12 @@ program test_distribution_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. + + logical :: warn = .true. integer :: put, get put = 12345678 @@ -25,11 +28,11 @@ program test_distribution_beta call test_beta_rvs_${t1[0]}$${k1}$ #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in BETA_REAL_KINDS_TYPES call test_beta_pdf_${t1[0]}$${k1}$ #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in BETA_REAL_KINDS_TYPES call test_beta_cdf_${t1[0]}$${k1}$ #:endfor @@ -53,215 +56,156 @@ contains expct = num / array_size do i = 0, array_size - 1 - chisq = chisq + (freq(i) - expct) ** 2 / expct + 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(*,*) "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) + 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 + integer, parameter :: n = 100 ${t1}$ :: res(n), ba, bb, loc integer :: i integer :: seed, get - #:if t1[0] == "r" - #! for real type - ${t1}$ :: ans(n) = [0.23516498458385468E+00_${k1}$, & - 0.38906532525908948E+00_${k1}$, & - 0.34268730878829956E+00_${k1}$, & - 0.29513829946517944E+00_${k1}$, & - 0.19447886384116858E+00_${k1}$, & - 0.60383307933807373E+00_${k1}$, & - 0.58347225189208984E+00_${k1}$, & - 0.74638533592224121E+00_${k1}$, & - 0.55743497610092163E+00_${k1}$, & - 0.39149415493011475E+00_${k1}$] - #:else - #! for complex type - ${t1}$ :: ans(n) = & - [(0.44851309061050415E+00_${k1}$, 0.34970587491989136E+00_${k1}$), & - (0.25030922889709473E+00_${k1}$, 0.55935305356979370E+00_${k1}$), & - (0.69629210233688354E+00_${k1}$, 0.25098574161529541E+00_${k1}$), & - (0.54261916875839233E+00_${k1}$, 0.45337703824043274E+00_${k1}$), & - (0.25677326321601868E+00_${k1}$, 0.18765588104724884E+00_${k1}$), & - (0.59211272001266479E+00_${k1}$, 0.72924208641052246E+00_${k1}$), & - (0.75178939104080200E+00_${k1}$, 0.18360115587711334E+00_${k1}$), & - (0.41408637166023254E+00_${k1}$, 0.36145636439323425E+00_${k1}$), & - (0.47387075424194336E+00_${k1}$, 0.26313120126724243E+00_${k1}$), & - (0.71398687362670898E+00_${k1}$, 0.51318806409835815E+00_${k1}$)] - #:endif print *, "Test beta_distribution_rvs_${t1[0]}$${k1}$" seed = 639741825 call random_seed(seed, get) #:if t1[0] == "r" - #! for real type - ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + ba = 2.0_${k1}$ + bb = 5.0_${k1}$ + loc = 0.0_${k1}$ #:else - #! for complex type - ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) + ba = (2.0_${k1}$, 0.7_${k1}$) + bb = (5.0_${k1}$, 3.0_${k1}$) + loc = (0.0_${k1}$, 0.0_${k1}$) #:endif - do i = 1, k + do i = 1, n res(i) = rbeta(ba, bb, loc) end do - res(k + 1 : n) = rbeta(ba, bb, k, loc) - + #:if t1[0] == "r" do i = 1, n - call check(abs(res(i) - ans(i)) < ${k1}$tol, & - msg="beta_distribution_rvs_${t1[0]}$${k1}$ failed", & - warn=warn) + call check((res(i) >= 0.0_${k1}$) .and. (res(i) <= 1.0_${k1}$), & + msg="beta_distribution_rvs_${t1[0]}$${k1}$ out of range", warn=warn) end do - end subroutine test_beta_rvs_${t1[0]}$${k1}$ + #:else + do i = 1, n + call check((real(res(i)) >= 0.0_${k1}$) .and. (real(res(i)) <= 1.0_${k1}$) .and. & + (aimag(res(i)) >= 0.0_${k1}$) .and. (aimag(res(i)) <= 1.0_${k1}$), & + msg="beta_distribution_rvs_${t1[0]}$${k1}$ out of range", warn=warn) + end do + #:endif + end subroutine test_beta_rvs_${t1[0]}$${k1}$ #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in BETA_REAL_KINDS_TYPES subroutine test_beta_pdf_${t1[0]}$${k1}$ - ${t1}$ :: x1, x2(3,4), ba, bb, loc integer :: i - integer :: seed, get real(${k1}$) :: res(15) - #:if t1[0] == "r" - #! for real type - real(${k1}$), parameter :: ans(15) = & - [1.23914337158203125E+00_${k1}$, & - 1.23914337158203125E+00_${k1}$, & - 1.23914337158203125E+00_${k1}$, & - 1.08365607261657715E+00_${k1}$, & - 1.04956150054931641E+00_${k1}$, & - 1.03561067581176758E+00_${k1}$, & - 1.22067928314208984E+00_${k1}$, & - 8.46962928771972656E-01_${k1}$, & - 2.52243280410766602E-01_${k1}$, & - 9.41652893066406250E-01_${k1}$, & - 6.57052993774414062E-01_${k1}$, & - 1.14095616340637207E+00_${k1}$, & - 9.09092426300048828E-01_${k1}$, & - 1.22890543937683105E+00_${k1}$, & - 4.40967559814453125E-01_${k1}$] - #:else - #! for complex type - real(${k1}$), parameter :: ans(15) = & - [1.15128993988037109E+00_${k1}$, & - 1.15128993988037109E+00_${k1}$, & - 1.15128993988037109E+00_${k1}$, & - 8.46223831176757812E-01_${k1}$, & - 3.30976486206054688E+00_${k1}$, & - 3.01892852783203125E+00_${k1}$, & - 1.23958110809326172E+00_${k1}$, & - 1.97296905517578125E+00_${k1}$, & - 1.17863655090332031E+00_${k1}$, & - 4.00447845458984375E-01_${k1}$, & - 2.92449951171875000E-02_${k1}$, & - 9.25743579864501953E-01_${k1}$, & - 2.15837764739990234E+00_${k1}$, & - 5.10101318359375000E-02_${k1}$, & - 1.03166580200195312E+00_${k1}$] - #:endif + real(${k1}$), parameter :: x_vals(15) = [ & + 0.05000000000000000E+00_${k1}$, & + 0.10000000000000001E+00_${k1}$, & + 0.15000000000000002E+00_${k1}$, & + 0.20000000000000001E+00_${k1}$, & + 0.25000000000000000E+00_${k1}$, & + 0.30000000000000004E+00_${k1}$, & + 0.35000000000000003E+00_${k1}$, & + 0.40000000000000002E+00_${k1}$, & + 0.45000000000000001E+00_${k1}$, & + 0.50000000000000000E+00_${k1}$, & + 0.55000000000000004E+00_${k1}$, & + 0.59999999999999998E+00_${k1}$, & + 0.70000000000000007E+00_${k1}$, & + 0.80000000000000004E+00_${k1}$, & + 0.90000000000000002E+00_${k1}$] + real(${k1}$), parameter :: ans(15) = [ & + 1.22175937500000087E+00_${k1}$, & + 1.96829999999999972E+00_${k1}$, & + 2.34902812499999891E+00_${k1}$, & + 2.45759999999999978E+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}$, & + 6.76603124999999861E-01_${k1}$, & + 4.60800000000000043E-01_${k1}$, & + 1.70099999999999946E-01_${k1}$, & + 3.83999999999999897E-02_${k1}$, & + 2.69999999999999928E-03_${k1}$] print *, "Test beta_distribution_pdf_${t1[0]}$${k1}$" - seed = 345987126 - call random_seed(seed, get) - #:if t1[0] == "r" - #! for real type - ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ - #:else - #! for complex type - ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) - #:endif - - x1 = rbeta(ba, bb, loc) - x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) - res(1:3) = beta_pdf(x1, ba, bb, loc) - res(4:15) = reshape(beta_pdf(x2, ba, bb, loc), [12]) + do i = 1, 15 + res(i) = beta_pdf(x_vals(i), 2.0_${k1}$, 5.0_${k1}$, 0.0_${k1}$) + end do do i = 1, 15 - call check(abs(res(i) - ans(i)) < ${k1}$tol, & - msg="beta_distribution_pdf_${t1[0]}$${k1}$ failed", & - warn=warn) + call check(abs(res(i) - ans(i)) < max(${k1}$tol, 1.0e-12_${k1}$), & + msg="beta_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) end do - end subroutine test_beta_pdf_${t1[0]}$${k1}$ + end subroutine test_beta_pdf_${t1[0]}$${k1}$ #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in BETA_REAL_KINDS_TYPES subroutine test_beta_cdf_${t1[0]}$${k1}$ - ${t1}$ :: x1, x2(3,4), ba, bb, loc integer :: i - integer :: seed, get real(${k1}$) :: res(15) - #:if t1[0] == "r" - #! for real type - real(${k1}$), parameter :: ans(15) = & - [0.37429335713386536E+00_${k1}$, & - 0.37429335713386536E+00_${k1}$, & - 0.37429335713386536E+00_${k1}$, & - 0.10568153858184814E+00_${k1}$, & - 0.61683785915374756E+00_${k1}$, & - 0.95208930969238281E-01_${k1}$, & - 0.34991249442100525E+00_${k1}$, & - 0.73912835121154785E+00_${k1}$, & - 0.95353639125823975E+00_${k1}$, & - 0.68856298923492432E+00_${k1}$, & - 0.36885082721710205E-01_${k1}$, & - 0.14866244792938232E+00_${k1}$, & - 0.74632394313812256E+00_${k1}$, & - 0.29785108566284180E+00_${k1}$, & - 0.14457702636718750E-01_${k1}$] - #:else - #! for complex type - real(${k1}$), parameter :: ans(15) = & - [0.14019477367401123E+00_${k1}$, & - 0.14019477367401123E+00_${k1}$, & - 0.14019477367401123E+00_${k1}$, & - 0.48782587051391602E+00_${k1}$, & - 0.19825363159179688E+00_${k1}$, & - 0.92399120330810547E-02_${k1}$, & - 0.55928134918212891E+00_${k1}$, & - 0.15336728096008301E+00_${k1}$, & - 0.31033802032470703E+00_${k1}$, & - 0.69235420227050781E+00_${k1}$, & - 0.92155790328979492E+00_${k1}$, & - 0.63153839111328125E+00_${k1}$, & - 0.86834430694580078E-01_${k1}$, & - 0.69353580474853516E+00_${k1}$, & - 0.59732055664062500E+00_${k1}$] - #:endif + real(${k1}$), parameter :: x_vals(15) = [ & + 0.05000000000000000E+00_${k1}$, & + 0.10000000000000001E+00_${k1}$, & + 0.15000000000000002E+00_${k1}$, & + 0.20000000000000001E+00_${k1}$, & + 0.25000000000000000E+00_${k1}$, & + 0.30000000000000004E+00_${k1}$, & + 0.35000000000000003E+00_${k1}$, & + 0.40000000000000002E+00_${k1}$, & + 0.45000000000000001E+00_${k1}$, & + 0.50000000000000000E+00_${k1}$, & + 0.55000000000000004E+00_${k1}$, & + 0.59999999999999998E+00_${k1}$, & + 0.70000000000000007E+00_${k1}$, & + 0.80000000000000004E+00_${k1}$, & + 0.90000000000000002E+00_${k1}$] + real(${k1}$), parameter :: ans(15) = [ & + 3.27738281249999944E-02_${k1}$, & + 1.14265000000000019E-01_${k1}$, & + 2.23515703125000020E-01_${k1}$, & + 3.44640000000000224E-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.30801953125000026E-01_${k1}$, & + 9.59040000000000004E-01_${k1}$, & + 9.89064999999999972E-01_${k1}$, & + 9.98399999999999954E-01_${k1}$, & + 9.99944999999999973E-01_${k1}$] print *, "Test beta_distribution_cdf_${t1[0]}$${k1}$" - seed = 345987126 - call random_seed(seed, get) - #:if t1[0] == "r" - #! for real type - ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ - #:else - #! for complex type - ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) - #:endif - - x1 = rbeta(ba, bb, loc) - x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) - res(1:3) = beta_cdf(x1, ba, bb, loc) - res(4:15) = reshape(beta_cdf(x2, ba, bb, loc), [12]) + do i = 1, 15 + res(i) = beta_cdf(x_vals(i), 2.0_${k1}$, 5.0_${k1}$, 0.0_${k1}$) + end do do i = 1, 15 - call check(abs(res(i) - ans(i)) < ${k1}$tol, & - msg="beta_distribution_cdf_${t1[0]}$${k1}$ failed", & - warn=warn) + call check(abs(res(i) - ans(i)) < max(${k1}$tol, 1.0e-12_${k1}$), & + msg="beta_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) end do - end subroutine test_beta_cdf_${t1[0]}$${k1}$ + end subroutine test_beta_cdf_${t1[0]}$${k1}$ #:endfor end program test_distribution_beta From b6a4460cb77321faba973587720d090d34104c66 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:48:38 +0000 Subject: [PATCH 15/33] last commit wrong draft | correct now --- test/stats/test_distribution_beta.fypp | 267 +++++++++++++++---------- 1 file changed, 162 insertions(+), 105 deletions(-) diff --git a/test/stats/test_distribution_beta.fypp b/test/stats/test_distribution_beta.fypp index d6777c207..682bff58f 100644 --- a/test/stats/test_distribution_beta.fypp +++ b/test/stats/test_distribution_beta.fypp @@ -11,12 +11,10 @@ program test_distribution_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. + logical :: warn = .true. integer :: put, get put = 12345678 @@ -28,11 +26,11 @@ program test_distribution_beta call test_beta_rvs_${t1[0]}$${k1}$ #:endfor - #:for k1, t1 in BETA_REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES call test_beta_pdf_${t1[0]}$${k1}$ #:endfor - #:for k1, t1 in BETA_REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES call test_beta_cdf_${t1[0]}$${k1}$ #:endfor @@ -56,156 +54,215 @@ contains expct = num / array_size do i = 0, array_size - 1 - chisq = chisq + (freq(i) - expct) ** 2 / expct + 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(*,*) "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) + 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 :: n = 100 + integer, parameter :: k = 5, n = 10 ${t1}$ :: res(n), ba, bb, loc integer :: i integer :: seed, get + #:if t1[0] == "r" + #! for real type + ${t1}$ :: ans(n) = [0.23516498458385468E+00_${k1}$, & + 0.38906532525908948E+00_${k1}$, & + 0.34268730878829956E+00_${k1}$, & + 0.29513829946517944E+00_${k1}$, & + 0.19447886384116858E+00_${k1}$, & + 0.60383307933807373E+00_${k1}$, & + 0.58347225189208984E+00_${k1}$, & + 0.74638533592224121E+00_${k1}$, & + 0.55743497610092163E+00_${k1}$, & + 0.39149415493011475E+00_${k1}$] + #:else + #! for complex type + ${t1}$ :: ans(n) = & + [(0.44851309061050415E+00_${k1}$, 0.34970587491989136E+00_${k1}$), & + (0.25030922889709473E+00_${k1}$, 0.55935305356979370E+00_${k1}$), & + (0.69629210233688354E+00_${k1}$, 0.25098574161529541E+00_${k1}$), & + (0.54261916875839233E+00_${k1}$, 0.45337703824043274E+00_${k1}$), & + (0.25677326321601868E+00_${k1}$, 0.18765588104724884E+00_${k1}$), & + (0.59211272001266479E+00_${k1}$, 0.72924208641052246E+00_${k1}$), & + (0.75178939104080200E+00_${k1}$, 0.18360115587711334E+00_${k1}$), & + (0.41408637166023254E+00_${k1}$, 0.36145636439323425E+00_${k1}$), & + (0.47387075424194336E+00_${k1}$, 0.26313120126724243E+00_${k1}$), & + (0.71398687362670898E+00_${k1}$, 0.51318806409835815E+00_${k1}$)] + #:endif 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.0_${k1}$ + #! for real type + 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.0_${k1}$, 0.0_${k1}$) + #! for complex type + ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) #:endif - do i = 1, n + do i = 1, k res(i) = rbeta(ba, bb, loc) end do - #:if t1[0] == "r" - do i = 1, n - call check((res(i) >= 0.0_${k1}$) .and. (res(i) <= 1.0_${k1}$), & - msg="beta_distribution_rvs_${t1[0]}$${k1}$ out of range", warn=warn) - end do - #:else + res(k + 1 : n) = rbeta(ba, bb, k, loc) + do i = 1, n - call check((real(res(i)) >= 0.0_${k1}$) .and. (real(res(i)) <= 1.0_${k1}$) .and. & - (aimag(res(i)) >= 0.0_${k1}$) .and. (aimag(res(i)) <= 1.0_${k1}$), & - msg="beta_distribution_rvs_${t1[0]}$${k1}$ out of range", warn=warn) + call check(abs(res(i) - ans(i)) < ${k1}$tol, & + msg="beta_distribution_rvs_${t1[0]}$${k1}$ failed", & + warn=warn) end do - #:endif - end subroutine test_beta_rvs_${t1[0]}$${k1}$ + #:endfor - #:for k1, t1 in BETA_REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES subroutine test_beta_pdf_${t1[0]}$${k1}$ + ${t1}$ :: x1, x2(3,4), ba, bb, loc integer :: i + integer :: seed, get real(${k1}$) :: res(15) - real(${k1}$), parameter :: x_vals(15) = [ & - 0.05000000000000000E+00_${k1}$, & - 0.10000000000000001E+00_${k1}$, & - 0.15000000000000002E+00_${k1}$, & - 0.20000000000000001E+00_${k1}$, & - 0.25000000000000000E+00_${k1}$, & - 0.30000000000000004E+00_${k1}$, & - 0.35000000000000003E+00_${k1}$, & - 0.40000000000000002E+00_${k1}$, & - 0.45000000000000001E+00_${k1}$, & - 0.50000000000000000E+00_${k1}$, & - 0.55000000000000004E+00_${k1}$, & - 0.59999999999999998E+00_${k1}$, & - 0.70000000000000007E+00_${k1}$, & - 0.80000000000000004E+00_${k1}$, & - 0.90000000000000002E+00_${k1}$] - real(${k1}$), parameter :: ans(15) = [ & - 1.22175937500000087E+00_${k1}$, & - 1.96829999999999972E+00_${k1}$, & - 2.34902812499999891E+00_${k1}$, & - 2.45759999999999978E+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}$, & - 6.76603124999999861E-01_${k1}$, & - 4.60800000000000043E-01_${k1}$, & - 1.70099999999999946E-01_${k1}$, & - 3.83999999999999897E-02_${k1}$, & - 2.69999999999999928E-03_${k1}$] + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [1.23914337158203125E+00_${k1}$, & + 1.23914337158203125E+00_${k1}$, & + 1.23914337158203125E+00_${k1}$, & + 1.08365607261657715E+00_${k1}$, & + 1.04956150054931641E+00_${k1}$, & + 1.03561067581176758E+00_${k1}$, & + 1.22067928314208984E+00_${k1}$, & + 8.46962928771972656E-01_${k1}$, & + 2.52243280410766602E-01_${k1}$, & + 9.41652893066406250E-01_${k1}$, & + 6.57052993774414062E-01_${k1}$, & + 1.14095616340637207E+00_${k1}$, & + 9.09092426300048828E-01_${k1}$, & + 1.22890543937683105E+00_${k1}$, & + 4.40967559814453125E-01_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [1.15128993988037109E+00_${k1}$, & + 1.15128993988037109E+00_${k1}$, & + 1.15128993988037109E+00_${k1}$, & + 8.46223831176757812E-01_${k1}$, & + 3.30976486206054688E+00_${k1}$, & + 3.01892852783203125E+00_${k1}$, & + 1.23958110809326172E+00_${k1}$, & + 1.97296905517578125E+00_${k1}$, & + 1.17863655090332031E+00_${k1}$, & + 4.00447845458984375E-01_${k1}$, & + 2.92449951171875000E-02_${k1}$, & + 9.25743579864501953E-01_${k1}$, & + 2.15837764739990234E+00_${k1}$, & + 5.10101318359375000E-02_${k1}$, & + 1.03166580200195312E+00_${k1}$] + #:endif print *, "Test beta_distribution_pdf_${t1[0]}$${k1}$" + seed = 345987126 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + #:else + #! for complex type + ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) + #:endif - do i = 1, 15 - res(i) = beta_pdf(x_vals(i), 2.0_${k1}$, 5.0_${k1}$, 0.0_${k1}$) - end do + x1 = rbeta(ba, bb, loc) + x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) + + res(1:3) = beta_pdf(x1, ba, bb, loc) + res(4:15) = reshape(beta_pdf(x2, ba, bb, loc), [12]) do i = 1, 15 - call check(abs(res(i) - ans(i)) < max(${k1}$tol, 1.0e-12_${k1}$), & - msg="beta_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) + call check(abs(res(i) - ans(i)) < ${k1}$tol, & + 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 BETA_REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES subroutine test_beta_cdf_${t1[0]}$${k1}$ + ${t1}$ :: x1, x2(3,4), ba, bb, loc integer :: i + integer :: seed, get real(${k1}$) :: res(15) - real(${k1}$), parameter :: x_vals(15) = [ & - 0.05000000000000000E+00_${k1}$, & - 0.10000000000000001E+00_${k1}$, & - 0.15000000000000002E+00_${k1}$, & - 0.20000000000000001E+00_${k1}$, & - 0.25000000000000000E+00_${k1}$, & - 0.30000000000000004E+00_${k1}$, & - 0.35000000000000003E+00_${k1}$, & - 0.40000000000000002E+00_${k1}$, & - 0.45000000000000001E+00_${k1}$, & - 0.50000000000000000E+00_${k1}$, & - 0.55000000000000004E+00_${k1}$, & - 0.59999999999999998E+00_${k1}$, & - 0.70000000000000007E+00_${k1}$, & - 0.80000000000000004E+00_${k1}$, & - 0.90000000000000002E+00_${k1}$] - real(${k1}$), parameter :: ans(15) = [ & - 3.27738281249999944E-02_${k1}$, & - 1.14265000000000019E-01_${k1}$, & - 2.23515703125000020E-01_${k1}$, & - 3.44640000000000224E-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.30801953125000026E-01_${k1}$, & - 9.59040000000000004E-01_${k1}$, & - 9.89064999999999972E-01_${k1}$, & - 9.98399999999999954E-01_${k1}$, & - 9.99944999999999973E-01_${k1}$] + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [0.37429335713386536E+00_${k1}$, & + 0.37429335713386536E+00_${k1}$, & + 0.37429335713386536E+00_${k1}$, & + 0.10568153858184814E+00_${k1}$, & + 0.61683785915374756E+00_${k1}$, & + 0.95208930969238281E-01_${k1}$, & + 0.34991249442100525E+00_${k1}$, & + 0.73912835121154785E+00_${k1}$, & + 0.95353639125823975E+00_${k1}$, & + 0.68856298923492432E+00_${k1}$, & + 0.36885082721710205E-01_${k1}$, & + 0.14866244792938232E+00_${k1}$, & + 0.74632394313812256E+00_${k1}$, & + 0.29785108566284180E+00_${k1}$, & + 0.14457702636718750E-01_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [0.14019477367401123E+00_${k1}$, & + 0.14019477367401123E+00_${k1}$, & + 0.14019477367401123E+00_${k1}$, & + 0.48782587051391602E+00_${k1}$, & + 0.19825363159179688E+00_${k1}$, & + 0.92399120330810547E-02_${k1}$, & + 0.55928134918212891E+00_${k1}$, & + 0.15336728096008301E+00_${k1}$, & + 0.31033802032470703E+00_${k1}$, & + 0.69235420227050781E+00_${k1}$, & + 0.92155790328979492E+00_${k1}$, & + 0.63153839111328125E+00_${k1}$, & + 0.86834430694580078E-01_${k1}$, & + 0.69353580474853516E+00_${k1}$, & + 0.59732055664062500E+00_${k1}$] + #:endif print *, "Test beta_distribution_cdf_${t1[0]}$${k1}$" + seed = 345987126 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + #:else + #! for complex type + ba = (2.0_${k1}$, 0.7_${k1}$); bb = (5.0_${k1}$, 3.0_${k1}$); loc = (0._${k1}$, 0._${k1}$) + #:endif - do i = 1, 15 - res(i) = beta_cdf(x_vals(i), 2.0_${k1}$, 5.0_${k1}$, 0.0_${k1}$) - end do + x1 = rbeta(ba, bb, loc) + x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) + + 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-12_${k1}$), & - msg="beta_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) + call check(abs(res(i) - ans(i)) < ${k1}$tol, & + 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 From 85cb4143bb1cd4d1b0b265dba097beb33dfa2d5c Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:52:10 +0000 Subject: [PATCH 16/33] expected results + cleanup --- test/stats/test_distribution_beta.fypp | 197 +++++++++++-------------- 1 file changed, 90 insertions(+), 107 deletions(-) diff --git a/test/stats/test_distribution_beta.fypp b/test/stats/test_distribution_beta.fypp index 682bff58f..a588a1b8e 100644 --- a/test/stats/test_distribution_beta.fypp +++ b/test/stats/test_distribution_beta.fypp @@ -14,7 +14,7 @@ program test_distribution_beta #:for k1, t1 in BETA_REAL_KINDS_TYPES ${t1}$, parameter :: ${k1}$tol = 1000 * epsilon(1.0_${k1}$) #:endfor - logical :: warn = .true. + logical :: warn = .true. integer :: put, get put = 12345678 @@ -68,54 +68,31 @@ contains #:for k1, t1 in RC_KINDS_TYPES subroutine test_beta_rvs_${t1[0]}$${k1}$ integer, parameter :: k = 5, n = 10 - ${t1}$ :: res(n), ba, bb, loc + ${t1}$ :: res(n), ans(n), ba, bb, loc integer :: i integer :: seed, get - #:if t1[0] == "r" - #! for real type - ${t1}$ :: ans(n) = [0.23516498458385468E+00_${k1}$, & - 0.38906532525908948E+00_${k1}$, & - 0.34268730878829956E+00_${k1}$, & - 0.29513829946517944E+00_${k1}$, & - 0.19447886384116858E+00_${k1}$, & - 0.60383307933807373E+00_${k1}$, & - 0.58347225189208984E+00_${k1}$, & - 0.74638533592224121E+00_${k1}$, & - 0.55743497610092163E+00_${k1}$, & - 0.39149415493011475E+00_${k1}$] - #:else - #! for complex type - ${t1}$ :: ans(n) = & - [(0.44851309061050415E+00_${k1}$, 0.34970587491989136E+00_${k1}$), & - (0.25030922889709473E+00_${k1}$, 0.55935305356979370E+00_${k1}$), & - (0.69629210233688354E+00_${k1}$, 0.25098574161529541E+00_${k1}$), & - (0.54261916875839233E+00_${k1}$, 0.45337703824043274E+00_${k1}$), & - (0.25677326321601868E+00_${k1}$, 0.18765588104724884E+00_${k1}$), & - (0.59211272001266479E+00_${k1}$, 0.72924208641052246E+00_${k1}$), & - (0.75178939104080200E+00_${k1}$, 0.18360115587711334E+00_${k1}$), & - (0.41408637166023254E+00_${k1}$, 0.36145636439323425E+00_${k1}$), & - (0.47387075424194336E+00_${k1}$, 0.26313120126724243E+00_${k1}$), & - (0.71398687362670898E+00_${k1}$, 0.51318806409835815E+00_${k1}$)] - #:endif print *, "Test beta_distribution_rvs_${t1[0]}$${k1}$" seed = 639741825 call random_seed(seed, get) #:if t1[0] == "r" - #! for real type ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ #:else - #! for complex type 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", & @@ -129,60 +106,63 @@ contains subroutine test_beta_pdf_${t1[0]}$${k1}$ ${t1}$ :: x1, x2(3,4), ba, bb, loc integer :: i - integer :: seed, get real(${k1}$) :: res(15) #:if t1[0] == "r" #! for real type real(${k1}$), parameter :: ans(15) = & - [1.23914337158203125E+00_${k1}$, & - 1.23914337158203125E+00_${k1}$, & - 1.23914337158203125E+00_${k1}$, & - 1.08365607261657715E+00_${k1}$, & - 1.04956150054931641E+00_${k1}$, & - 1.03561067581176758E+00_${k1}$, & - 1.22067928314208984E+00_${k1}$, & - 8.46962928771972656E-01_${k1}$, & - 2.52243280410766602E-01_${k1}$, & - 9.41652893066406250E-01_${k1}$, & - 6.57052993774414062E-01_${k1}$, & - 1.14095616340637207E+00_${k1}$, & - 9.09092426300048828E-01_${k1}$, & - 1.22890543937683105E+00_${k1}$, & - 4.40967559814453125E-01_${k1}$] + [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(15) = & - [1.15128993988037109E+00_${k1}$, & - 1.15128993988037109E+00_${k1}$, & - 1.15128993988037109E+00_${k1}$, & - 8.46223831176757812E-01_${k1}$, & - 3.30976486206054688E+00_${k1}$, & - 3.01892852783203125E+00_${k1}$, & - 1.23958110809326172E+00_${k1}$, & - 1.97296905517578125E+00_${k1}$, & - 1.17863655090332031E+00_${k1}$, & - 4.00447845458984375E-01_${k1}$, & - 2.92449951171875000E-02_${k1}$, & - 9.25743579864501953E-01_${k1}$, & - 2.15837764739990234E+00_${k1}$, & - 5.10101318359375000E-02_${k1}$, & - 1.03166580200195312E+00_${k1}$] + [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}$" - seed = 345987126 - call random_seed(seed, get) #:if t1[0] == "r" - #! for real type 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 - #! for complex type 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 - x1 = rbeta(ba, bb, loc) - x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) - res(1:3) = beta_pdf(x1, ba, bb, loc) res(4:15) = reshape(beta_pdf(x2, ba, bb, loc), [12]) @@ -199,60 +179,63 @@ contains subroutine test_beta_cdf_${t1[0]}$${k1}$ ${t1}$ :: x1, x2(3,4), ba, bb, loc integer :: i - integer :: seed, get real(${k1}$) :: res(15) #:if t1[0] == "r" #! for real type real(${k1}$), parameter :: ans(15) = & - [0.37429335713386536E+00_${k1}$, & - 0.37429335713386536E+00_${k1}$, & - 0.37429335713386536E+00_${k1}$, & - 0.10568153858184814E+00_${k1}$, & - 0.61683785915374756E+00_${k1}$, & - 0.95208930969238281E-01_${k1}$, & - 0.34991249442100525E+00_${k1}$, & - 0.73912835121154785E+00_${k1}$, & - 0.95353639125823975E+00_${k1}$, & - 0.68856298923492432E+00_${k1}$, & - 0.36885082721710205E-01_${k1}$, & - 0.14866244792938232E+00_${k1}$, & - 0.74632394313812256E+00_${k1}$, & - 0.29785108566284180E+00_${k1}$, & - 0.14457702636718750E-01_${k1}$] + [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) = & - [0.14019477367401123E+00_${k1}$, & - 0.14019477367401123E+00_${k1}$, & - 0.14019477367401123E+00_${k1}$, & - 0.48782587051391602E+00_${k1}$, & - 0.19825363159179688E+00_${k1}$, & - 0.92399120330810547E-02_${k1}$, & - 0.55928134918212891E+00_${k1}$, & - 0.15336728096008301E+00_${k1}$, & - 0.31033802032470703E+00_${k1}$, & - 0.69235420227050781E+00_${k1}$, & - 0.92155790328979492E+00_${k1}$, & - 0.63153839111328125E+00_${k1}$, & - 0.86834430694580078E-01_${k1}$, & - 0.69353580474853516E+00_${k1}$, & - 0.59732055664062500E+00_${k1}$] + [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}$" - seed = 345987126 - call random_seed(seed, get) #:if t1[0] == "r" - #! for real type 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 - #! for complex type 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 - x1 = rbeta(ba, bb, loc) - x2 = reshape(rbeta(ba, bb, 12, loc), [3,4]) - res(1:3) = beta_cdf(x1, ba, bb, loc) res(4:15) = reshape(beta_cdf(x2, ba, bb, loc), [12]) From 74c681ac9e812150336535138b869afefe2b5962 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:52:55 +0000 Subject: [PATCH 17/33] indent fix --- test/stats/test_distribution_beta.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/stats/test_distribution_beta.fypp b/test/stats/test_distribution_beta.fypp index a588a1b8e..ed47ffeba 100644 --- a/test/stats/test_distribution_beta.fypp +++ b/test/stats/test_distribution_beta.fypp @@ -167,7 +167,7 @@ contains res(4:15) = reshape(beta_pdf(x2, ba, bb, loc), [12]) do i = 1, 15 - call check(abs(res(i) - ans(i)) < ${k1}$tol, & + 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 @@ -240,7 +240,7 @@ contains res(4:15) = reshape(beta_cdf(x2, ba, bb, loc), [12]) do i = 1, 15 - call check(abs(res(i) - ans(i)) < ${k1}$tol, & + 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 From ba23bc8483f1a202262f67a5125b05c188c33b32 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 17:43:45 +0000 Subject: [PATCH 18/33] fixes bounds error --- example/stats_distribution_beta/example_beta_pdf.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/example/stats_distribution_beta/example_beta_pdf.f90 b/example/stats_distribution_beta/example_beta_pdf.f90 index 8a15d8bf4..452cc8519 100644 --- a/example/stats_distribution_beta/example_beta_pdf.f90 +++ b/example/stats_distribution_beta/example_beta_pdf.f90 @@ -6,7 +6,7 @@ program example_beta_pdf implicit none real :: x, a, b real :: xarr(3, 4, 5) - integer :: seed_put, seed_get, i + integer :: seed_put, seed_get seed_put = 1234567 call random_seed(seed_put, seed_get) @@ -22,9 +22,7 @@ program example_beta_pdf ! 1.32859993 ! generate random variates and compute their pdf - do i = 1, 60 - xarr(i, 1, 1) = rbeta(a, b) - end do + xarr = reshape(rbeta(a, b, 60), [3, 4, 5]) print *, beta_pdf(xarr, a, b) ! 1.23914337 0.925268888 0.969913423 0.998693407 1.44889069 From 777b0a3cd04554be857dcb9b89744189608a407e Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 9 Mar 2026 17:44:14 +0000 Subject: [PATCH 19/33] fixes bounds error --- example/stats_distribution_beta/example_beta_cdf.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/example/stats_distribution_beta/example_beta_cdf.f90 b/example/stats_distribution_beta/example_beta_cdf.f90 index 70e022e49..cb32e9334 100644 --- a/example/stats_distribution_beta/example_beta_cdf.f90 +++ b/example/stats_distribution_beta/example_beta_cdf.f90 @@ -6,7 +6,7 @@ program example_beta_cdf implicit none real :: x, a, b real :: xarr(3, 4, 5) - integer :: seed_put, seed_get, i + integer :: seed_put, seed_get seed_put = 1234567 call random_seed(seed_put, seed_get) @@ -22,9 +22,7 @@ program example_beta_cdf ! 0.471808970 ! generate random variates and compute their cdf - do i = 1, 60 - xarr(i, 1, 1) = rbeta(a, b) - end do + xarr = reshape(rbeta(a, b, 60), [3, 4, 5]) print *, beta_cdf(xarr, a, b) ! 0.374293357 0.136472717 0.153627276 0.166885555 0.535913110 From 6d769fde4f83370580b443ebc7018d81f390b7b7 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Thu, 12 Mar 2026 11:31:28 +0000 Subject: [PATCH 20/33] loc clarification --- doc/specs/stdlib_stats_distribution_beta.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_beta.md b/doc/specs/stdlib_stats_distribution_beta.md index 419853d2d..d0f4ae7a1 100644 --- a/doc/specs/stdlib_stats_distribution_beta.md +++ b/doc/specs/stdlib_stats_distribution_beta.md @@ -44,7 +44,7 @@ If `a` is `real`, its value must be positive. If `a` is `complex`, both the real 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. +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. @@ -118,7 +118,7 @@ Experimental ### Description -Cumulative distribution function (cdf) of the single real variable beta distribution: +Cumulative distribution function (cdf) of the single real variable beta distribution is: $$ F(x)= I_x(a, b),\\quad 00,\\ b>0 $$ From d325bd1da1f8d98ee175cf150b10d07641b0f5eb Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Thu, 12 Mar 2026 11:48:46 +0000 Subject: [PATCH 21/33] fixed fypp template from paste error --- .../test_specialfunctions_gamma.fypp | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 6272c31df..ec09e55d7 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -404,6 +404,11 @@ contains end do end subroutine test_gamma_q_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES[:-1] + subroutine test_beta_${t1[0]}$${k1}$(error) @@ -472,14 +477,6 @@ contains end do end subroutine test_incomplete_beta_${t1[0]}$${k1}$ - #:endfor -end module test_specialfunctions_gamma - #:endfor - - - - #:for k1, t1 in REAL_KINDS_TYPES[:-1] - subroutine test_lincgamma_${t1[0]}$${k1}$(error) type(error_type), allocatable, intent(out) :: error integer, parameter :: n = 4 From 72a00a2b582e87187635f295fed1373ff83f585d Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Sun, 15 Mar 2026 14:58:36 +0000 Subject: [PATCH 22/33] quiet fix to keep acuraccy while finding convergence Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../stdlib_specialfunctions_gamma.fypp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp index d04ce4f31..ae942b37d 100644 --- a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp +++ b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp @@ -1390,6 +1390,7 @@ contains ${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 @@ -1402,6 +1403,7 @@ contains if(abs(d) < fpmin) d = fpmin d = one / d h = d + converged = .false. do m = 1, maxit m2 = 2 * m @@ -1420,10 +1422,17 @@ contains d = one / d del = d * c h = h * del - if(abs(del - one) < eps) exit + if(abs(del - one) < eps) then + converged = .true. + exit + end if end do - res = h + if (converged) then + res = h + else + res = ieee_value(one, ieee_quiet_nan) + end if end function betacf_${t1[0]}$${k1}$ #:endfor From 3c736e894d4d47a76889fc9c955915c28ec08209 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Sun, 15 Mar 2026 15:10:29 +0000 Subject: [PATCH 23/33] drops the impure on incomplete beta --- src/specialfunctions/stdlib_specialfunctions_gamma.fypp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp index ae942b37d..fb4c7ddd3 100644 --- a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp +++ b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp @@ -1287,6 +1287,8 @@ contains #:endfor #:endfor + + ! ! Beta function implementations ! @@ -1329,7 +1331,7 @@ contains #:for k1, t1 in GAMMA_REAL_KINDS_TYPES - impure elemental function incomplete_beta_${t1[0]}$${k1}$(x, a, b) result(res) + 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 ! From c6604c1591b14ae091d2a06a2521e484c72ecc3a Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Sun, 15 Mar 2026 15:14:15 +0000 Subject: [PATCH 24/33] uses log beta same approach as gamma numerical stability change suggested by copilot --- src/stats/stdlib_stats_distribution_beta.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stats/stdlib_stats_distribution_beta.fypp b/src/stats/stdlib_stats_distribution_beta.fypp index fd3eff4aa..226589c3c 100644 --- a/src/stats/stdlib_stats_distribution_beta.fypp +++ b/src/stats/stdlib_stats_distribution_beta.fypp @@ -9,7 +9,7 @@ Module stdlib_stats_distribution_beta 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 + use stdlib_specialfunctions_gamma, only : beta, incomplete_beta, log_beta implicit none private @@ -208,7 +208,7 @@ contains end if ! Use log formulation for numerical stability - res = exp((a - one) * log(xs) + (b - one) * log(one - xs) - log(beta(a, b))) + res = exp((a - one) * log(xs) + (b - one) * log(one - xs) - log_beta(a, b)) end function beta_dist_pdf_${t1[0]}$${k1}$ #:endfor From fc3de44d534dada753b6b9903e3e6d8ca7ed970d Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Sun, 15 Mar 2026 15:17:50 +0000 Subject: [PATCH 25/33] error handling + elemental types should address the suggestions --- src/stats/stdlib_stats_distribution_beta.fypp | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/stats/stdlib_stats_distribution_beta.fypp b/src/stats/stdlib_stats_distribution_beta.fypp index 226589c3c..4001a5e1a 100644 --- a/src/stats/stdlib_stats_distribution_beta.fypp +++ b/src/stats/stdlib_stats_distribution_beta.fypp @@ -225,8 +225,7 @@ contains real(${k1}$) :: res res = beta_dist_pdf_${t1[0]}$${k1}$(x, a, b, loc) - err = 0 - if(ieee_is_nan(res)) err = 1 + err = merge(1, 0, ieee_is_nan(res)) end function beta_dist_pdf_impure_${t1[0]}$${k1}$ #:endfor @@ -268,15 +267,14 @@ contains 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 = 0 - if(ieee_is_nan(res)) err = 1 + 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 - impure elemental function beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc) & + elemental function beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc) & result(res) ! Beta distribution cumulative distribution function ! @@ -322,15 +320,14 @@ contains real(${k1}$) :: res res = beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc) - err = 0 - if(ieee_is_nan(res)) err = 1 + 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 - impure elemental function beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc) & + 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. @@ -364,8 +361,7 @@ contains loc_ = optval(loc, cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$)) res = beta_dist_cdf_${t1[0]}$${k1}$(x, a, b, loc_) - err = 0 - if(ieee_is_nan(res)) err = 1 + err = merge(1, 0, ieee_is_nan(res)) end function beta_dist_cdf_impure_${t1[0]}$${k1}$ #:endfor From 0191ae8fb1042449e8c8fc55dfe5c6087eeed62c Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 13:41:08 +0000 Subject: [PATCH 26/33] beta/log_beta changes for simplicity --- .../stdlib_specialfunctions_gamma.fypp | 28 ++++++------------- 1 file changed, 8 insertions(+), 20 deletions(-) diff --git a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp index fb4c7ddd3..d60f6ae07 100644 --- a/src/specialfunctions/stdlib_specialfunctions_gamma.fypp +++ b/src/specialfunctions/stdlib_specialfunctions_gamma.fypp @@ -1302,11 +1302,10 @@ contains if(a <= 0.0_${k1}$ .or. b <= 0.0_${k1}$) then res = ieee_value(1.0_${k1}$, ieee_quiet_nan) - return + else + ! Use log-gamma for numerical stability + res = exp(log_gamma(a) + log_gamma(b) - log_gamma(a + b)) end if - - ! Use log-gamma for numerical stability - res = exp(log_gamma(a) + log_gamma(b) - log_gamma(a + b)) end function beta_${t1[0]}$${k1}$ #:endfor @@ -1321,10 +1320,9 @@ contains if(a <= 0.0_${k1}$ .or. b <= 0.0_${k1}$) then res = ieee_value(1.0_${k1}$, ieee_quiet_nan) - return + else + res = log_gamma(a) + log_gamma(b) - log_gamma(a + b) end if - - res = log_gamma(a) + log_gamma(b) - log_gamma(a + b) end function log_beta_${t1[0]}$${k1}$ #:endfor @@ -1342,23 +1340,13 @@ contains ${t1}$, parameter :: eps = epsilon(1.0_${k1}$) integer, parameter :: maxit = 200 - if(a <= zero .or. b <= zero) then - res = ieee_value(1.0_${k1}$, ieee_quiet_nan) - return - end if - - if(x < zero .or. x > one) then + 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) then - res = zero - return - end if - - if(x == one) then - res = one + if(x == zero .or. x == one) then + res = x return end if From 9940c0ac0d1c7d6d7ae877974b6a047baf859397 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 13:47:06 +0000 Subject: [PATCH 27/33] tweaked beta cdf example made a and b paramters and reduced the elements --- example/stats_distribution_beta/example_beta_cdf.f90 | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/example/stats_distribution_beta/example_beta_cdf.f90 b/example/stats_distribution_beta/example_beta_cdf.f90 index cb32e9334..20a7a263f 100644 --- a/example/stats_distribution_beta/example_beta_cdf.f90 +++ b/example/stats_distribution_beta/example_beta_cdf.f90 @@ -4,15 +4,13 @@ program example_beta_cdf beta_cdf => cdf_beta implicit none - real :: x, a, b - real :: xarr(3, 4, 5) + 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) - a = 2.0; b = 5.0 - ! cumulative probability at x=0.3 for beta(2,5) distribution print *, beta_cdf(0.3, a, b) ! 0.471808970 @@ -22,13 +20,10 @@ program example_beta_cdf ! 0.471808970 ! generate random variates and compute their cdf - xarr = reshape(rbeta(a, b, 60), [3, 4, 5]) + xarr = reshape(rbeta(a, b, 10), [2, 5]) print *, beta_cdf(xarr, a, b) ! 0.374293357 0.136472717 0.153627276 0.166885555 0.535913110 ! 4.47619371E-02 0.161991328 0.524897814 7.37934634E-02 8.41872990E-02 - ! 0.102836564 0.138294637 0.122145072 0.171486780 0.532901883 - ! 0.111447386 1.87034653E-02 0.119697235 0.550687969 9.98932496E-02 - ! 0.234867483 0.132088020 0.210879743 0.520989060 ... (60 elements total) end program example_beta_cdf From 15534da6a56a695a4bbe7e40a40d9fef1a0c8042 Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 13:48:01 +0000 Subject: [PATCH 28/33] refreshed comments to match new output --- example/stats_distribution_beta/example_beta_cdf.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/example/stats_distribution_beta/example_beta_cdf.f90 b/example/stats_distribution_beta/example_beta_cdf.f90 index 20a7a263f..19de5c4e2 100644 --- a/example/stats_distribution_beta/example_beta_cdf.f90 +++ b/example/stats_distribution_beta/example_beta_cdf.f90 @@ -13,17 +13,17 @@ program example_beta_cdf ! cumulative probability at x=0.3 for beta(2,5) distribution print *, beta_cdf(0.3, a, b) - ! 0.471808970 + ! 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.471808970 + ! 0.579824746 ! generate random variates and compute their cdf xarr = reshape(rbeta(a, b, 10), [2, 5]) print *, beta_cdf(xarr, a, b) - ! 0.374293357 0.136472717 0.153627276 0.166885555 0.535913110 - ! 4.47619371E-02 0.161991328 0.524897814 7.37934634E-02 8.41872990E-02 + ! 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 From 50f7ba33eb290ccf636def2c24fbea9836a0e5bb Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 14:00:47 +0000 Subject: [PATCH 29/33] applied same as previous change to cdf --- .../example_beta_pdf.f90 | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/example/stats_distribution_beta/example_beta_pdf.f90 b/example/stats_distribution_beta/example_beta_pdf.f90 index 452cc8519..e7c34f974 100644 --- a/example/stats_distribution_beta/example_beta_pdf.f90 +++ b/example/stats_distribution_beta/example_beta_pdf.f90 @@ -4,31 +4,26 @@ program example_beta_pdf beta_pdf => pdf_beta implicit none - real :: x, a, b - real :: xarr(3, 4, 5) + 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) - a = 2.0; b = 5.0 - ! probability density at x=0.3 for beta(2,5) distribution print *, beta_pdf(0.3, a, b) - ! 1.32859993 + ! 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) - ! 1.32859993 + ! 2.16090035 ! generate random variates and compute their pdf - xarr = reshape(rbeta(a, b, 60), [3, 4, 5]) + xarr = reshape(rbeta(a, b, 10), [2, 5]) print *, beta_pdf(xarr, a, b) - ! 1.23914337 0.925268888 0.969913423 0.998693407 1.44889069 - ! 0.331506640 0.988773048 1.41048801 0.502468705 0.564825535 - ! 0.661424160 0.931534827 0.861649513 1.01826906 1.44227338 - ! 0.782453180 0.145308748 0.841068327 1.49033546 0.727149904 - ! 1.14145494 0.905850768 1.08809102 1.40226245 ... (60 elements total) + ! 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 From 42de7af21923889209331499ae808248985e115f Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 14:05:14 +0000 Subject: [PATCH 30/33] capitalisation fix Co-authored-by: Jeremie Vandenplas --- src/stats/stdlib_stats_distribution_beta.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stats/stdlib_stats_distribution_beta.fypp b/src/stats/stdlib_stats_distribution_beta.fypp index 4001a5e1a..5f99bd7c1 100644 --- a/src/stats/stdlib_stats_distribution_beta.fypp +++ b/src/stats/stdlib_stats_distribution_beta.fypp @@ -2,7 +2,7 @@ #: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 +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 From 2dbcdb523594b519cbf77e9ba60760535d6fb28b Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 14:08:18 +0000 Subject: [PATCH 31/33] moved error return before opt value parse --- src/stats/stdlib_stats_distribution_beta.fypp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stats/stdlib_stats_distribution_beta.fypp b/src/stats/stdlib_stats_distribution_beta.fypp index 5f99bd7c1..0ad5bac43 100644 --- a/src/stats/stdlib_stats_distribution_beta.fypp +++ b/src/stats/stdlib_stats_distribution_beta.fypp @@ -2,7 +2,7 @@ #: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 +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 @@ -194,14 +194,14 @@ contains ${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 + loc_ = optval(loc, 0.0_${k1}$) + xs = x - loc_ + if(xs <= zero .or. xs >= one) then res = zero return From 08ec3b85d1107f9dd384c5ba5d680f34d4662c1d Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 14:13:18 +0000 Subject: [PATCH 32/33] changed error returns + parameter definitions --- .../test_specialfunctions_gamma.fypp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index e60038ecc..96979a9a4 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -415,8 +415,8 @@ contains type(error_type), allocatable, intent(out) :: error integer, parameter :: n = 3 integer :: i - ${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 :: 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}$, & @@ -427,6 +427,7 @@ contains 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}$ @@ -437,6 +438,7 @@ contains 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}$] @@ -445,10 +447,12 @@ contains 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), & - "Log-beta function with a(kind=${k1}$) and b(kind=${k1}$) failed", & + msg, & thr = tol_${k1}$, rel = .true.) + if (allocated(error)) return end do end subroutine test_log_beta_${t1[0]}$${k1}$ @@ -473,6 +477,7 @@ contains 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}$ @@ -646,4 +651,4 @@ program tester write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" error stop end if -end program tester \ No newline at end of file +end program tester From 8bf15dc3a121d65b077f7df679ad0b0e4d60a11b Mon Sep 17 00:00:00 2001 From: Callum Winder <69708064+Calluumm@users.noreply.github.com> Date: Mon, 16 Mar 2026 14:23:02 +0000 Subject: [PATCH 33/33] parameters over values + misc change --- test/stats/test_distribution_beta.fypp | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/test/stats/test_distribution_beta.fypp b/test/stats/test_distribution_beta.fypp index ed47ffeba..b33532a91 100644 --- a/test/stats/test_distribution_beta.fypp +++ b/test/stats/test_distribution_beta.fypp @@ -105,11 +105,12 @@ contains #: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(15) + real(${k1}$) :: res(n) #:if t1[0] == "r" #! for real type - real(${k1}$), parameter :: ans(15) = & + real(${k1}$), parameter :: ans(*) = & [2.45759999999999978E+00_${k1}$, & 2.45759999999999978E+00_${k1}$, & 2.45759999999999978E+00_${k1}$, & @@ -127,7 +128,7 @@ contains 3.83999999999999897E-02_${k1}$] #:else #! for complex type - real(${k1}$), parameter :: ans(15) = & + real(${k1}$), parameter :: ans(*) = & [2.77620563793055197E+00_${k1}$, & 2.77620563793055197E+00_${k1}$, & 2.77620563793055197E+00_${k1}$, & @@ -147,13 +148,17 @@ contains print *, "Test beta_distribution_pdf_${t1[0]}$${k1}$" #:if t1[0] == "r" - ba = 2.0_${k1}$; bb = 5.0_${k1}$; loc = 0._${k1}$ + 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}$) + 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}$), & @@ -166,7 +171,7 @@ contains res(1:3) = beta_pdf(x1, ba, bb, loc) res(4:15) = reshape(beta_pdf(x2, ba, bb, loc), [12]) - do i = 1, 15 + 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)