add: beta probability distribution

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

doc/specs/stdlib_stats_distribution_beta.md

@@ -0,0 +1,167 @@
+---
+title: stats_distribution_beta
+---
+
+# Statistical Distributions -- Beta Distribution Module
+
+[TOC]
+
+## `rvs_beta` - beta distribution random variates
+
+### Status
+
+Experimental
+
+### Description
+
+The beta distribution is a continuous probability distribution defined on the interval [0, 1], widely used for modeling random variables that represent proportions, probabilities, and other bounded quantities. It is defined by two shape parameters (\\(a\\) and \\(b\\)) that control the distribution's form.
+
+With two arguments (a, b), the function returns a random sample from the beta distribution \\(\\text{Beta}(a, b)\\).
+
+The optional `loc` parameter specifies the location (shift) of the distribution.
+
+With three or more arguments including `array_size`, the function returns a rank-1 array of beta distributed random variates.
+
+For complex shape parameters, the real and imaginary parts are sampled independently of each other.
+
+@note
+For shape parameters less than 1, the function uses a uniform method. For parameters greater than or equal to 1, it uses the gamma ratio method[^1], where \\(X \\sim \\text{Beta}(a,b)\\) is generated as \\(X = \\frac{Y_1}{Y_1 + Y_2}\\) where \\(Y_1 \\sim \\Gamma(a,1)\\) and \\(Y_2 \\sim \\Gamma(b,1)\\).
+
+### Syntax
+
+`result = [[stdlib_stats_distribution_beta(module):rvs_beta(interface)]](a, b [[, loc]] [[, array_size]])`
+
+### Class
+
+Impure elemental function
+
+### Arguments
+
+`a`: has `intent(in)` and is a scalar of type `real` or `complex`. 
+If `a` is `real`, its value must be positive. If `a` is `complex`, both the real and imaginary components must be positive. This is the first shape parameter of the distribution.
+
+`b`: has `intent(in)` and is a scalar of type `real` or `complex`. 
+If `b` is `real`, its value must be positive. If `b` is `complex`, both the real and imaginary components must be positive. This is the second shape parameter of the distribution.
+
+`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. 
+Specifies the location (shift) of the distribution with default value 0.0. The distribution support is loc < x < loc + 1.
+
+`array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind.
+
+### Return value
+
+The result is a scalar or rank-1 array with a size of `array_size`, and the same type as `a`. If `a` or `b` is non-positive, the result is `NaN`.
+
+### Example
+
+```fortran
+{!example/stats_distribution_beta/example_beta_rvs.f90!}
+```
+
+## `pdf_beta` - beta distribution probability density function
+
+### Status
+
+Experimental
+
+### Description
+
+The probability density function (pdf) of the single real variable beta distribution is:
+
+$$ f(x)= \\frac{x^{a-1}(1-x)^{b-1}}{B(a,b)} ,\\quad 0<x<1,\\ a>0,\\ b>0 $$
+
+where \\(a\\) and \\(b\\) are the shape parameters, and \\(B(a,b)\\) is the beta function.
+
+An optional `loc` parameter specifies the location (shift) of the distribution.
+
+For a complex variable \\(z=(x + y i)\\) with independent real \\(x\\) and imaginary \\(y\\) parts, the joint probability density function is the product of the corresponding real and imaginary marginal pdfs:[^2]
+
+$$f(x+\\mathit{i}y)=f(x)f(y)$$
+
+### Syntax
+
+`result = [[stdlib_stats_distribution_beta(module):pdf_beta(interface)]](x, a, b [[, loc]])`
+
+### Class
+
+Impure elemental function
+
+### Arguments
+
+`x`: has `intent(in)` and is a scalar of type `real` or `complex`. The point at which to evaluate the pdf.
+
+`a`: has `intent(in)` and is a scalar of type `real` or `complex`. The first shape parameter. 
+If `a` is `real`, its value must be positive. If `a` is `complex`, both the real and imaginary components must be positive.
+
+`b`: has `intent(in)` and is a scalar of type `real` or `complex`. The second shape parameter. 
+If `b` is `real`, its value must be positive. If `b` is `complex`, both the real and imaginary components must be positive.
+
+`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. The location (shift) parameter with default value 0.0.
+
+All arguments must have the same type.
+
+### Return value
+
+The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `a` or `b` is non-positive, the result is `NaN`.
+
+### Example
+
+```fortran
+{!example/stats_distribution_beta/example_beta_pdf.f90!}
+```
+
+## `cdf_beta` - beta distribution cumulative distribution function
+
+### Status
+
+Experimental
+
+### Description
+
+Cumulative distribution function (cdf) of the single real variable beta distribution is:
+
+$$ F(x)= I_x(a, b),\\quad 0<x<1,\\ a>0,\\ 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).

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)

example/stats_distribution_beta/CMakeLists.txt

@@ -0,0 +1,4 @@
+ADD_EXAMPLE(beta_rvs)
+ADD_EXAMPLE(beta_pdf)
+ADD_EXAMPLE(beta_cdf)
+

example/stats_distribution_beta/example_beta_cdf.f90

@@ -0,0 +1,34 @@
+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
+
+  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
+  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
+  ! 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

example/stats_distribution_beta/example_beta_pdf.f90

@@ -0,0 +1,34 @@
+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
+
+  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
+  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
+  ! 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

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