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

@@ -113,4 +113,33 @@ ${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 Fortran expressions.
+#!
+#! Args:
+#!   varname (str): Name of the variable to be used as origin
+#!   varname1 (str): Name of the variable to be used instead of varname
+#!   origrank (int): Rank of the original variable
+#!   dim (int): Index of the used expression varname1
+#!
+#! Returns:
+#!   Shape expression enclosed in braces, except for the index dim.
+#!
+#! E.g., (:, :, :, i, :, :)
+#!
+
+#:def rankindice(varname, varname1, origrank, dim)
+  #:assert origrank > 0
+  #:if origrank > 0
+    #:call join_lines(joinstr=", ", prefix="(", suffix=")")
+      #:for i in range(1, origrank+1)
+        #:if i == dim
+          ${varname1}$
+        #:else
+          ${varname}$
+        #:endif
+      #:endfor
+    #:endcall
+  #:endif
+#: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,106 @@ 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
+          ${t1}$ :: 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
+          ${t1}$ :: 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)}$
+          ${t1}$ :: 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)}$
+          ${t1}$ :: 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
+```
+