Merge pull request #144 from jvdp1/variance_dev

Addition of variance function in stdlib_experimental_stats

src/CMakeLists.txt

@@ -6,6 +6,7 @@ set(fppFiles
     stdlib_experimental_optval.fypp
     stdlib_experimental_stats.fypp
     stdlib_experimental_stats_mean.fypp
+    stdlib_experimental_stats_var.fypp
 )
 
 

src/Makefile.manual

@@ -5,7 +5,8 @@ SRC = stdlib_experimental_ascii.f90 \
       stdlib_experimental_kinds.f90 \
       f18estop.f90 \
       stdlib_experimental_stats.f90 \
-      stdlib_experimental_stats_mean.f90
+      stdlib_experimental_stats_mean.f90 \
+      stdlib_experimental_stats_var.f90
 
 LIB = libstdlib.a
 
@@ -42,6 +43,11 @@ stdlib_experimental_stats_mean.o: \
 	stdlib_experimental_optval.o \
 	stdlib_experimental_kinds.o \
 	stdlib_experimental_stats.o
+stdlib_experimental_stats_var.o: \
+	stdlib_experimental_optval.o \
+	stdlib_experimental_kinds.o \
+	stdlib_experimental_stats.o
 stdlib_experimental_io.f90: stdlib_experimental_io.fypp
 stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp
 stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp
+stdlib_experimental_stats_var.f90: stdlib_experimental_stats_var.fypp

src/common.fypp

@@ -97,6 +97,7 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
   #:endif
 #:enddef
 
+
 #! Generates a routine name from a generic name, rank, type and kind
 #!
 #! Args:
@@ -113,4 +114,28 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
   $:"{0}_{1}_{2}{3}_{2}{3}".format(gname, rank, type[0], kind) if suffix == '' else "{0}_{1}_{2}{3}_{4}".format(gname, rank, type[0], kind, suffix)
 #:enddef
 
+
+#! Generates an array rank suffix for subarrays reducing the dimension
+#!
+#! Args:
+#!   rank (int): Rank of the original variable
+#!   selectors (array): Dimension and name of the variable(s)
+#!
+#! Returns:
+#!   Array rank suffix string enclosed in braces
+#!
+#! E.g.,
+#!   select_subarray(5 , [(4, 'i'), (5, 'j')])}$
+#!   -> (:, :, :, i, j)
+#!
+#:def select_subarray(rank, selectors)
+  #:assert rank > 0
+  #:set seldict = dict(selectors)
+  #:call join_lines(joinstr=", ", prefix="(", suffix=")")
+    #:for i in range(1, rank + 1)
+      $:seldict.get(i, ":")
+    #:endfor
+  #:endcall
+#:enddef
+
 #:endmute

src/stdlib_experimental_stats.fypp

@@ -7,7 +7,7 @@ module stdlib_experimental_stats
   implicit none
   private
   ! Public API
-  public :: mean
+  public :: mean, var
 
   interface mean
     #:for k1, t1 in RC_KINDS_TYPES
@@ -104,4 +104,102 @@ module stdlib_experimental_stats
 
   end interface mean
 
+
+  interface var
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          logical, intent(in), optional :: mask
+          real(${k1}$) :: res
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          logical, intent(in), optional :: mask
+          real(dp) :: res
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: dim
+          logical, intent(in), optional :: mask
+          real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: dim
+          logical, intent(in), optional :: mask
+          real(dp) :: res${reduced_shape('x', rank, 'dim')}$
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          logical, intent(in) :: mask${ranksuffix(rank)}$
+          real(${k1}$) :: res
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          logical, intent(in) :: mask${ranksuffix(rank)}$
+          real(dp) :: res
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: dim
+          logical, intent(in) :: mask${ranksuffix(rank)}$
+          real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+    #: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)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: dim
+          logical, intent(in) :: mask${ranksuffix(rank)}$
+          real(dp) :: res${reduced_shape('x', rank, 'dim')}$
+        end function ${RName}$
+      #:endfor
+    #:endfor
+
+  end interface var
+
+
 end module stdlib_experimental_stats

src/stdlib_experimental_stats.md

@@ -3,6 +3,7 @@
 ## Implemented
 
  * `mean`
+ * `var`
 
 ## `mean` - mean of array elements
 
@@ -47,3 +48,55 @@ program demo_mean
                    reshape(x, [ 2, 3 ] ) > 3.)  !returns [ NaN, 4.0, 5.5 ]
 end program demo_mean
 ```
+
+## `var` - variance of array elements
+
+### Description
+
+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:
+
+```
+ var(x) = 1/(n-1) sum_i (array(i) - mean(array))^2
+```
+
+### Syntax
+
+`result = var(array [, mask])`
+
+`result = var(array, dim [, mask])`
+
+### Arguments
+
+`array`: Shall be an array of type `integer`, `real`, or `complex`.
+
+`dim`: Shall be a scalar of type `integer` with a value in the range from 1 to n, where n is the rank of `array`.
+
+`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
+
+### 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 `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 `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`.
+
+### Example
+
+```fortran
+program demo_var
+    use stdlib_experimental_stats, only: var
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    print *, var(x)                            !returns 3.5
+    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]
+end program demo_var
+```
+