Probability Distribution and Statistical Functions -- Uniform Distribution Module (fortran-lang#593)

* add pdf/cdf with corresponding precisions
* correction for kinds
* add description of variable loc for pdf/cdf

Author: Jim-215-Fisher

doc/specs/stdlib_stats_distribution_uniform.md

@@ -36,7 +36,7 @@ Return a randomized rank one array of the input type.
 
 ```fortran
 program demo_shuffle
-    use stdlib_stats_distribution_PRNG, only : random_seed
+    use stdlib_random, only : random_seed
     use stdlib_stats_distribution_uniform, only : shuffle
     implicit none
     integer :: seed_put, seed_get, i
@@ -79,14 +79,16 @@ Without argument the function returns a scalar standard uniformly distributed va
 
 With single argument `scale` of `integer` type the function returns a scalar uniformly distributed variate of `integer` type on [0,scale]. This is the standard Rectangular distribution.
 
-With single argument `scale` of `real` or `complex` type the function returns a scalar uniformly distributed variate of `real` or `complex` type on [0, scale].
+With single argument `scale` of `real` or `complex` type the function returns a scalar uniformly distributed variate of `real` type on [0, scale] or `complex` type on [(0, 0i), (scale, i(scale))].
 
-With double arguments `loc` and `scale` the function returns a scalar uniformly distributed random variates of `integer`, `real` or `complex` type on [loc, loc + scale] dependent of input type.
+With double arguments `loc` and `scale` the function returns a scalar uniformly distributed random variates of `integer` or `real` type on [loc, loc + scale], or `complex` type on [(loc, i(loc)), ((loc + scale), i(loc + scale))], dependent of input type.
 
 With triple arguments `loc`, `scale` and `array_size` the function returns a rank one array of uniformly distributed variates of `integer`, `real` or `complex` type with an array size of `array_size`.
 
 For `complex` type, the real part and imaginary part are independent of each other.
 
+Note: the algorithm used for generating uniform random variates is fundamentally limited to double precision.
+
 ### Syntax
 
 `result = [[stdlib_stats_distribution_uniform(module):rvs_uniform(interface)]]([[loc,] scale] [[[,array_size]]])`
@@ -101,19 +103,19 @@ Elemental function (without the third argument).
 
 `scale`: optional argument has `intent(in)` and is a scalar of type `integer`, `real` or `complex`.
 
-`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`.
+`array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind.
 
 `loc` and `scale` must have the same type and kind when both are present.
 
 ### Return value
 
-The result is a scalar or a rank one array, with size of `array_size`, of type `integer`, `real` or `complex` depending on the input type.
+The result is a scalar or a rank one array with size of `array_size`, of type `integer`, `real` or `complex` depending on the input type.
 
 ### Example
 
 ```fortran
 program demo_uniform_rvs
-    use stdlib_stats_distribution_PRNG, only:random_seed
+    use stdlib_random, only:random_seed
     use stdlib_stats_distribution_uniform, only:uni=> rvs_uniform
 
     implicit none
@@ -193,15 +195,19 @@ program demo_uniform_rvs
 end program demo_uniform_rvs
 ```
 
-## `pdf_uniform` - Uniform probability density function
+## `pdf_uniform` - Uniform distribution probability density function
 
 ### Status
 
 Experimental
 
 ### Description
 
-The probability density function of the uniform distribution.
+The probability density function of the uniform distribution:
+
+f(x) = 0       x < loc or x > loc + scale  for all types uniform distributions
+
+For random variable x in [loc, loc + scale]:
 
 f(x) = 1 / (scale + 1);            for discrete uniform distribution.
 
@@ -235,7 +241,7 @@ The result is a scalar or an array, with a shape conformable to arguments, of ty
 
 ```fortran
 program demo_uniform_pdf
-    use stdlib_stats_distribution_PRNG, only : random_seed
+    use stdlib_random, only : random_seed
     use stdlib_stats_distribution_uniform,  only : uni_pdf => pdf_uniform,     &
                                                    uni => rvs_uniform
 
@@ -286,15 +292,21 @@ end program demo_uniform_pdf
 
 ```
 
-## `cdf_uniform` - Uniform cumulative distribution function
+## `cdf_uniform` - Uniform distribution cumulative distribution function
 
 ### Status
 
 Experimental
 
 ### Description
 
-Cumulative distribution function of the uniform distribution
+Cumulative distribution function of the uniform distribution:
+
+F(x) = 0             x < loc             for all types uniform distributions
+
+F(x) = 1             x > loc + scale     for all types uniform distributions
+
+For random variable x in [loc, loc + scale]:
 
 F(x) = (x - loc + 1) / (scale + 1);      for discrete uniform distribution.
 
@@ -328,7 +340,7 @@ The result is a scalar or an array, with a shape conformable to arguments, of ty
 
 ```fortran
 program demo_uniform_cdf
-    use stdlib_stats_distribution_PRNG, only : random_seed
+    use stdlib_random, only : random_seed
     use stdlib_stats_distribution_uniform, only : uni_cdf => cdf_uniform,      &
                                                   uni => rvs_uniform
 

src/stdlib_stats_distribution_uniform.fypp

@@ -2,7 +2,7 @@
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
 #:set ALL_KINDS_TYPES = INT_KINDS_TYPES + RC_KINDS_TYPES
 module stdlib_stats_distribution_uniform
-    use stdlib_kinds
+    use stdlib_kinds, only : sp, dp, xdp, qp, int8, int16, int32, int64
     use stdlib_error, only : error_stop
     use stdlib_random, only : dist_rand
 
@@ -383,14 +383,15 @@ contains
     elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
 
         ${t1}$, intent(in) :: x, loc, scale
-        real :: res
+        ${t1}$ :: res
+        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
 
-        if(scale == 0.0_${k1}$) then
-            res = 0.0
+        if(scale == zero) then
+            res = zero
         else if(x < loc .or. x > (loc + scale)) then
-            res = 0.0
+            res = zero
         else
-            res = 1.0 / scale
+            res = one / scale
         end if
     end function pdf_unif_${t1[0]}$${k1}$
 
@@ -402,17 +403,17 @@ contains
     elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
 
         ${t1}$, intent(in) :: x, loc, scale
-        real :: res
-        real(${k1}$) :: tr, ti
+        real(${k1}$) :: res, tr, ti
+        real(${k1}$), parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
 
         tr = loc % re + scale % re; ti = loc % im + scale % im
-        if(scale == (0.0_${k1}$,0.0_${k1}$)) then
-            res = 0.0
+        if(scale == (zero, zero)) then
+            res = zero
         else if((x % re >= loc % re .and. x % re <= tr) .and.                  &
             (x % im >= loc % im .and. x % im <= ti)) then
-            res = 1.0 / (scale % re * scale % im)
+            res = one / (scale % re * scale % im)
         else
-            res = 0.0
+            res = zero
         end if
     end function pdf_unif_${t1[0]}$${k1}$
 
@@ -445,16 +446,17 @@ contains
     elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
 
         ${t1}$, intent(in) :: x, loc, scale
-        real :: res
+        ${t1}$ :: res
+        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
 
-        if(scale == 0.0_${k1}$) then
-            res = 0.0
+        if(scale == zero) then
+            res = zero
         else if(x < loc) then
-            res = 0.0
+            res = zero
         else if(x >= loc .and. x <= (loc + scale)) then
             res = (x - loc) / scale
         else
-            res = 1.0
+            res = one
         end if
     end function cdf_unif_${t1[0]}$${k1}$
 
@@ -466,29 +468,30 @@ contains
     elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
 
         ${t1}$, intent(in) :: x, loc, scale
-        real :: res
+        real(${k1}$) :: res
+        real(${k1}$), parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
         logical :: r1, r2, i1, i2
 
-        if(scale == (0.0_${k1}$,0.0_${k1}$)) then
-            res = 0.0
+        if(scale == (zero, zero)) then
+            res = zero
             return
         end if
         r1 = x % re < loc % re
         r2 = x % re > (loc % re + scale % re)
         i1 = x % im < loc % im
         i2 = x % im > (loc % im + scale % im)
         if(r1 .or. i1) then
-            res = 0.0
+            res = zero
         else if((.not. r1) .and. (.not. r2) .and. i2) then
             res = (x % re - loc % re) / scale % re
         else if((.not. i1) .and. (.not. i2) .and. r2) then
             res = (x % im - loc % im) / scale % im
         else if((.not. r1) .and. (.not. r2) .and. (.not. i1) .and. (.not. i2)) &
             then
-            res = (x % re - loc % re) * (x % im - loc % im) /                  &
-                   (scale % re * scale % im)
+            res = ((x % re - loc % re) / scale % re) * ((x % im - loc % im) /  &
+                  scale % im)
         else if(r2 .and. i2)then
-             res = 1.0
+             res = one
         end if
     end function cdf_unif_${t1[0]}$${k1}$
 

src/tests/stats/test_distribution_uniform.fypp

@@ -207,16 +207,18 @@ contains
     subroutine test_uni_pdf_${t1[0]}$${k1}$
         ${t1}$ :: x1, x2(3,4), loc, scale
         integer :: seed, get, i
-        real :: res(3,5)
         #:if t1[0] == "i"
             #! for integer type
+        real :: res(3, 5)
         real, parameter :: ans(15) = [(1.96078438E-02, i=1,15)]
         #:elif t1[0] == "r"
             #! for real type
-        real, parameter :: ans(15) = [(0.5, i=1,15)]
+        ${t1}$ :: res(3, 5)
+        ${t1}$, parameter :: ans(15) = [(0.5_${k1}$, i=1,15)]
         #:else
             #! for complex type
-        real, parameter :: ans(15) = [(1.0, i=1,15)]
+        real(${k1}$) :: res(3, 5)
+        real(${k1}$), parameter :: ans(15) = [(1.0_${k1}$, i=1,15)]
         #:endif
 
         print *, "Test pdf_uniform_${t1[0]}$${k1}$"
@@ -236,8 +238,15 @@ contains
         x2 = reshape(uni_rvs(loc, scale, 12), [3,4])
         res(:,1) = uni_pdf(x1, loc, scale)
         res(:, 2:5) = uni_pdf(x2, loc, scale)
+        #:if t1[0] == "i"
+            #! for integer type
         call check(all(abs(res - reshape(ans,[3,5])) < sptol),                 &
             msg = "pdf_uniform_${t1[0]}$${k1}$ failed", warn=warn)
+        #:else
+            #! for real and complex types
+        call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol),             &
+            msg = "pdf_uniform_${t1[0]}$${k1}$ failed", warn=warn)
+        #:endif
     end subroutine test_uni_pdf_${t1[0]}$${k1}$
 
     #:endfor
@@ -247,47 +256,73 @@ contains
     #:for k1, t1 in ALL_KINDS_TYPES
     subroutine test_uni_cdf_${t1[0]}$${k1}$
         ${t1}$ :: x1, x2(3,4), loc, scale
-        real :: res(3,5)
         integer :: seed, get
         #:if k1 == "int8"
+        real :: res(3,5)
         real, parameter :: ans(15) = [0.435643554, 0.435643554, 0.435643554,   &
                                       0.702970326, 0.653465331, 0.485148519,   &
                                       0.386138618, 0.386138618, 0.336633652,   &
                                       0.277227730, 0.237623766, 0.524752498,   &
                                       0.732673287, 0.534653485, 0.415841579]
         #:elif k1 == "int16"
+        real :: res(3,5)
         real, parameter :: ans(15) = [0.178217828, 0.178217828, 0.178217828,   &
                                       0.465346545, 0.673267305, 0.247524753,   &
                                       0.158415839, 0.792079210, 0.742574275,   &
                                       0.574257433, 0.881188095, 0.663366318,   &
                                       0.524752498, 0.623762369, 0.178217828]
         #:elif k1 == "int32"
+        real :: res(3,5)
         real, parameter :: ans(15) = [0.732673287, 0.732673287, 0.732673287,   &
                                       0.722772300, 0.792079210, 5.94059415E-02,&
                                       0.841584146, 0.405940592, 0.960396051,   &
                                       0.534653485, 0.782178223, 0.861386120,   &
                                       0.564356446, 0.613861382, 0.306930691]
         #:elif k1 == "int64"
+        real :: res(3,5)
         real, parameter :: ans(15) = [0.455445558, 0.455445558, 0.455445558,   &
                                       0.277227730, 0.455445558, 0.930693090,   &
                                       0.851485133, 0.623762369, 5.94059415E-02,&
                                       0.693069279, 0.544554472, 0.207920790,   &
                                       0.306930691, 0.356435657, 0.128712878]
         #:elif t1[0] == "r"
             #! for real type
-        real, parameter :: ans(15) = [0.170192942, 0.170192942, 0.170192942,   &
-                                      0.276106149, 0.754930079, 0.406620681,   &
-                                      0.187742814, 0.651605546, 0.934733927,   &
-                                      0.151271492, 0.987674534, 0.130533904,   &
-                                      0.106271908, 9.27578658E-02, 0.203399420]
+        ${t1}$ :: res(3,5)
+        ${t1}$, parameter :: ans(15) =                                         &
+                             [0.170192944297557408050991512027394492_${k1}$,   &
+                              0.170192944297557408050991512027394492_${k1}$,   &
+                              0.170192944297557408050991512027394492_${k1}$,   &
+                              0.276106146274646191418611351764411665_${k1}$,   &
+                              0.754930097473875072466853453079238534_${k1}$,   &
+                              0.406620682573118008562573777453508228_${k1}$,   &
+                              0.187742819191801080247472555129206739_${k1}$,   &
+                              0.651605526090507591874256831943057477_${k1}$,   &
+                              0.934733949732104885121941606485052034_${k1}$,   &
+                              0.151271491851613287815681019310432021_${k1}$,   &
+                              0.987674522797719611766353864368284121_${k1}$,   &
+                              0.130533899463404684526679488953959662_${k1}$,   &
+                              0.106271905921876880229959283497009892_${k1}$,   &
+                         9.27578652240113182836367400341259781E-0002_${k1}$,   &
+                              0.203399426853420439709196898547816090_${k1}$]
         #:else
             #! for complex type
-        real, parameter :: ans(15) = [4.69913185E-02, 4.69913185E-02,          &
-                                      4.69913185E-02, 0.306970179, 0.122334257,&
-                                      0.141398594, 0.128925011, 9.85755492E-03,&
-                                      8.16527531E-02, 0.163921610,             &
-                                      7.81712309E-02, 0.446415812,             &
-                                      5.31753292E-04, 0.101455867, 0.155276477]
+        real(${k1}$) :: res(3,5)
+        real(${k1}$), parameter :: ans(15) =                                   &
+                        [4.69913179731340971083526490627998168E-0002_${k1}$,   &
+                         4.69913179731340971083526490627998168E-0002_${k1}$,   &
+                         4.69913179731340971083526490627998168E-0002_${k1}$,   &
+                              0.306970191529817593217448363707416739_${k1}$,   &
+                              0.122334258469188588238756489506609443_${k1}$,   &
+                              0.141398599060326408705075175176932616_${k1}$,   &
+                              0.128925006861443729884744412460848140_${k1}$,   &
+                         9.85755512660026594506599410104817938E-0003_${k1}$,   &
+                         8.16527497645585445208592050401597260E-0002_${k1}$,   &
+                              0.163921605454974749736935624944263178_${k1}$,   &
+                         7.81712317416218284294000447064256003E-0002_${k1}$,   &
+                              0.446415807686727375005224206895756087_${k1}$,   &
+                         5.31753272901435018841591264266743165E-0004_${k1}$,   &
+                              0.101455865191097416942685556683943046_${k1}$,   &
+                              0.155276470981788516449112374966730510_${k1}$]
         #:endif
 
         print *, "Test cdf_uniform_${t1[0]}$${k1}$"
@@ -307,8 +342,15 @@ contains
         x2 = reshape(uni_rvs(loc, scale, 12), [3,4])
         res(:,1) = uni_cdf(x1, loc, scale)
         res(:, 2:5) = uni_cdf(x2, loc, scale)
+        #:if t1[0] == "i"
+            #! for integer type
         call check(all(abs(res - reshape(ans,[3,5])) < sptol),                 &
             msg = "cdf_uniform_${t1[0]}$${k1}$ failed", warn=warn)
+        #:else
+            #! for real and complex types
+        call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol),             &
+            msg = "cdf_uniform_${t1[0]}$${k1}$ failed", warn=warn)
+        #:endif
     end subroutine test_uni_cdf_${t1[0]}$${k1}$
 
     #:endfor