Merge remote-tracking branch 'upstream/master' into moment_dev

jvdp1 · Feb 28, 2020 · dfbdd97 · dfbdd97
2 parents 195c62d + 5f35833
commit dfbdd97
Show file tree

Hide file tree

Showing 6 changed files with 467 additions and 94 deletions.
diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp
@@ -110,9 +110,10 @@ module stdlib_experimental_stats
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var_all",rank, t1, k1)
-        module function ${RName}$(x, mask) result(res)
+        module function ${RName}$(x, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           logical, intent(in), optional :: mask
+          logical, intent(in), optional :: corrected
           real(${k1}$) :: res
         end function ${RName}$
       #:endfor
@@ -121,9 +122,10 @@ module stdlib_experimental_stats
     #:for k1, t1 in INT_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var_all",rank, t1, k1, 'dp')
-        module function ${RName}$(x, mask) result(res)
+        module function ${RName}$(x, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           logical, intent(in), optional :: mask
+          logical, intent(in), optional :: corrected
           real(dp) :: res
         end function ${RName}$
       #:endfor
@@ -132,10 +134,11 @@ module stdlib_experimental_stats
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var",rank, t1, k1)
-        module function ${RName}$(x, dim, mask) result(res)
+        module function ${RName}$(x, dim, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           integer, intent(in) :: dim
           logical, intent(in), optional :: mask
+          logical, intent(in), optional :: corrected
           real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
         end function ${RName}$
       #:endfor
@@ -144,10 +147,11 @@ module stdlib_experimental_stats
     #:for k1, t1 in INT_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var",rank, t1, k1, 'dp')
-        module function ${RName}$(x, dim, mask) result(res)
+        module function ${RName}$(x, dim, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           integer, intent(in) :: dim
           logical, intent(in), optional :: mask
+          logical, intent(in), optional :: corrected
           real(dp) :: res${reduced_shape('x', rank, 'dim')}$
         end function ${RName}$
       #:endfor
@@ -156,9 +160,10 @@ module stdlib_experimental_stats
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var_mask_all",rank, t1, k1)
-        module function ${RName}$(x, mask) result(res)
+        module function ${RName}$(x, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           logical, intent(in) :: mask${ranksuffix(rank)}$
+          logical, intent(in), optional :: corrected
           real(${k1}$) :: res
         end function ${RName}$
       #:endfor
@@ -167,9 +172,10 @@ module stdlib_experimental_stats
     #:for k1, t1 in INT_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
-        module function ${RName}$(x, mask) result(res)
+        module function ${RName}$(x, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           logical, intent(in) :: mask${ranksuffix(rank)}$
+          logical, intent(in), optional :: corrected
           real(dp) :: res
         end function ${RName}$
       #:endfor
@@ -178,10 +184,11 @@ module stdlib_experimental_stats
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var_mask",rank, t1, k1)
-        module function ${RName}$(x, dim, mask) result(res)
+        module function ${RName}$(x, dim, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           integer, intent(in) :: dim
           logical, intent(in) :: mask${ranksuffix(rank)}$
+          logical, intent(in), optional :: corrected
           real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
         end function ${RName}$
       #:endfor
@@ -190,10 +197,11 @@ module stdlib_experimental_stats
     #:for k1, t1 in INT_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("var_mask",rank, t1, k1, 'dp')
-        module function ${RName}$(x, dim, mask) result(res)
+        module function ${RName}$(x, dim, mask, corrected) result(res)
           ${t1}$, intent(in) :: x${ranksuffix(rank)}$
           integer, intent(in) :: dim
           logical, intent(in) :: mask${ranksuffix(rank)}$
+          logical, intent(in), optional :: corrected
           real(dp) :: res${reduced_shape('x', rank, 'dim')}$
         end function ${RName}$
       #:endfor

diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md
@@ -55,17 +55,22 @@ end program demo_mean
 
 Returns the variance of all the elements of `array`, or of the elements of `array` along dimension `dim` if provided, and if the corresponding element in `mask` is `true`.
 
-The variance is defined as the best unbiased estimator and is computed as:
+Per default, the variance is defined as the best unbiased estimator and is computed as:
 
 ```
  var(x) = 1/(n-1) sum_i (array(i) - mean(array))^2
 ```
 
+where n is the number of elements.
+
+The use of the term `n-1` for scaling is called Bessel 's correction. The scaling can be changed with the logical argument `corrected`. If `corrected` is `.false.`, then the sum is scaled with `n`, otherwise with `n-1`.
+
+
 ### Syntax
 
-`result = var(array [, mask])`
+`result = var(array [, mask [, corrected]])`
 
-`result = var(array, dim [, mask])`
+`result = var(array, dim [, mask [, corrected]])`
 
 ### Arguments
 
@@ -75,16 +80,19 @@ The variance is defined as the best unbiased estimator and is computed as:
 
 `mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
 
+`corrected` (optional): Shall be a scalar of type `logical`. If `corrected` is `.true.` (default value), the sum is scaled with `n-1`. If `corrected` is `.false.`, then the sum is scaled with `n`.
+
 ### Return value
 
-If `array` is of type `real` or `complex`, the result is of the same type as `array`.
-If `array` is of type `integer`, the result is of type `double precision`.
+If `array` is of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
+If `array` is of type `integer`, the result is of type `real(dp)`.
 
-If `dim` is absent, a scalar with the variance of all elements in `array` is returned. Otherwise, an array of rank n-1, where n equals the rank of `array`, and a shape similar to that of `ar
-ray` with dimension `dim` dropped is returned.
+If `dim` is absent, a scalar with the variance of all elements in `array` is returned. Otherwise, an array of rank n-1, where n equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
 
 If `mask` is specified, the result is the variance of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`.
 
+If the variance is computed with only one single element, then the result is IEEE `NaN` if `corrected` is `.true.` and is `0.` if `corrected` is `.false.`.
+
 ### Example
 
 ```fortran
@@ -93,10 +101,14 @@ program demo_var
     implicit none
     real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
     print *, var(x)                            !returns 3.5
+    print *, var(x, corrected = .false.)       !returns 2.9167
     print *, var( reshape(x, [ 2, 3 ] ))       !returns 3.5
     print *, var( reshape(x, [ 2, 3 ] ), 1)    !returns [0.5, 0.5, 0.5]
     print *, var( reshape(x, [ 2, 3 ] ), 1,&
-                  reshape(x, [ 2, 3 ] ) > 3.)  !returns [0., NaN, 0.5]
+                  reshape(x, [ 2, 3 ] ) > 3.)  !returns [NaN, NaN, 0.5]
+    print *, var( reshape(x, [ 2, 3 ] ), 1,&
+                  reshape(x, [ 2, 3 ] ) > 3.,&
+                  corrected=.false.)           !returns [NaN, 0., 0.5]
 end program demo_var
 ```
 
diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp
@@ -13,9 +13,10 @@ contains
   #:for k1, t1 in RC_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var_all",rank, t1, k1)
-      module function ${RName}$(x, mask) result(res)
+      module function ${RName}$(x, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
         real(${k1}$) :: res
 
         real(${k1}$) :: n
@@ -26,13 +27,13 @@ contains
           return
         end if
 
-        n = real(size(x, kind = int64), ${k1}$)
+        n = size(x, kind = int64)
         mean = sum(x) / n
 
         #:if t1[0] == 'r'
-          res = sum((x - mean)**2) / (n - 1._${k1}$)
+          res = sum((x - mean)**2) / (n - merge(1, 0 , optval(corrected, .true.)))
         #:else
-          res = sum(abs(x - mean)**2) / (n - 1._${k1}$)
+          res = sum(abs(x - mean)**2) / (n - merge(1, 0, optval(corrected, .true.)))
         #:endif
 
       end function ${RName}$
@@ -43,9 +44,10 @@ contains
   #:for k1, t1 in INT_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var_all",rank, t1, k1, 'dp')
-      module function ${RName}$(x, mask) result(res)
+      module function ${RName}$(x, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
         real(dp) :: res
 
         real(dp) :: n, mean
@@ -55,10 +57,10 @@ contains
           return
         end if
 
-        n = real(size(x, kind = int64), dp)
+        n = size(x, kind = int64)
         mean = sum(real(x, dp)) / n
 
-        res = sum((real(x, dp) - mean)**2) / (n - 1._dp)
+        res = sum((real(x, dp) - mean)**2) / (n - merge(1, 0, optval(corrected, .true.)))
 
       end function ${RName}$
     #:endfor
@@ -68,10 +70,11 @@ contains
   #:for k1, t1 in RC_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var",rank, t1, k1)
-      module function ${RName}$(x, dim, mask) result(res)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         integer, intent(in) :: dim
         logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
         real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
 
         integer :: i
@@ -87,7 +90,7 @@ contains
         select case(dim)
           #:for fi in range(1, rank+1)
           case(${fi}$)
-            n = real(size(x, dim), ${k1}$)
+            n = size(x, dim)
             mean = sum(x, dim) / n
             do i = 1, size(x, dim)
               #:if t1[0] == 'r'
@@ -100,7 +103,7 @@ contains
           case default
             call error_stop("ERROR (mean): wrong dimension")
         end select
-        res = res / (n - 1._${k1}$)
+        res = res / (n - merge(1, 0, optval(corrected, .true.)))
 
       end function ${RName}$
     #:endfor
@@ -110,10 +113,11 @@ contains
   #:for k1, t1 in INT_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var",rank, t1, k1, 'dp')
-      module function ${RName}$(x, dim, mask) result(res)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         integer, intent(in) :: dim
         logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
         real(dp) :: res${reduced_shape('x', rank, 'dim')}$
 
         integer :: i
@@ -129,7 +133,7 @@ contains
         select case(dim)
           #:for fi in range(1, rank+1)
           case(${fi}$)
-            n = real(size(x, dim), dp)
+            n = size(x, dim)
             mean = sum(real(x, dp), dim) / n
             do i = 1, size(x, dim)
               res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2
@@ -138,7 +142,7 @@ contains
           case default
             call error_stop("ERROR (mean): wrong dimension")
         end select
-        res = res / (n - 1._dp)
+        res = res / (n - merge(1, 0, optval(corrected, .true.)))
 
       end function ${RName}$
     #:endfor
@@ -148,22 +152,24 @@ contains
   #:for k1, t1 in RC_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var_mask_all",rank, t1, k1)
-      module function ${RName}$(x, mask) result(res)
+      module function ${RName}$(x, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         logical, intent(in) :: mask${ranksuffix(rank)}$
+        logical, intent(in), optional :: corrected
         real(${k1}$) :: res
 
         real(${k1}$) :: n
         ${t1}$ :: mean
 
-        n = real(count(mask, kind = int64), ${k1}$)
+        n = count(mask, kind = int64)
         mean = sum(x, mask) / n
 
         #:if t1[0] == 'r'
-          res = sum((x - mean)**2, mask) / (n - 1._${k1}$)
+          res = sum((x - mean)**2, mask) / (n -&
         #:else
-          res = sum(abs(x - mean)**2, mask) / (n - 1._${k1}$)
+          res = sum(abs(x - mean)**2, mask) / (n -&
         #:endif
+                merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
 
       end function ${RName}$
     #:endfor
@@ -173,17 +179,19 @@ contains
   #:for k1, t1 in INT_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
-      module function ${RName}$(x, mask) result(res)
+      module function ${RName}$(x, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         logical, intent(in) :: mask${ranksuffix(rank)}$
+        logical, intent(in), optional :: corrected
         real(dp) :: res
 
         real(dp) :: n, mean
 
-        n = real(count(mask, kind = int64), dp)
+        n = count(mask, kind = int64)
         mean = sum(real(x, dp), mask) / n
 
-        res = sum((real(x, dp) - mean)**2, mask) / (n - 1._dp)
+        res = sum((real(x, dp) - mean)**2, mask) / (n -&
+              merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
 
       end function ${RName}$
     #:endfor
@@ -193,10 +201,11 @@ contains
   #:for k1, t1 in RC_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var_mask",rank, t1, k1)
-      module function ${RName}$(x, dim, mask) result(res)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         integer, intent(in) :: dim
         logical, intent(in) :: mask${ranksuffix(rank)}$
+        logical, intent(in), optional :: corrected
         real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
 
         integer :: i
@@ -207,7 +216,7 @@ contains
         select case(dim)
           #:for fi in range(1, rank+1)
           case(${fi}$)
-            n = real(count(mask, dim), ${k1}$)
+            n = count(mask, dim)
             mean = sum(x, dim, mask) / n
             do i = 1, size(x, dim)
               #:if t1[0] == 'r'
@@ -222,7 +231,7 @@ contains
           case default
             call error_stop("ERROR (mean): wrong dimension")
         end select
-        res = res / (n - 1._${k1}$)
+        res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
 
       end function ${RName}$
     #:endfor
@@ -232,10 +241,11 @@ contains
   #:for k1, t1 in INT_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("var_mask",rank, t1, k1, 'dp')
-      module function ${RName}$(x, dim, mask) result(res)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
         ${t1}$, intent(in) :: x${ranksuffix(rank)}$
         integer, intent(in) :: dim
         logical, intent(in) :: mask${ranksuffix(rank)}$
+        logical, intent(in), optional :: corrected
         real(dp) :: res${reduced_shape('x', rank, 'dim')}$
 
         integer :: i
@@ -246,7 +256,7 @@ contains
         select case(dim)
           #:for fi in range(1, rank+1)
           case(${fi}$)
-            n = real(count(mask, dim), dp)
+            n = count(mask, dim)
             mean = sum(real(x, dp), dim, mask) / n
             do i = 1, size(x, dim)
               res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2,&
@@ -256,7 +266,7 @@ contains
           case default
             call error_stop("ERROR (mean): wrong dimension")
         end select
-        res = res / (n - 1._dp)
+        res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
 
       end function ${RName}$
     #:endfor