Merge pull request #153 from jvdp1/central_moment

Addition of a function to compute the k-th order central moment of an array

Author: jvdp1

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_moment.fypp
     stdlib_experimental_stats_var.fypp
     stdlib_experimental_quadrature.fypp
     stdlib_experimental_quadrature_trapz.fypp

src/Makefile.manual

@@ -8,6 +8,7 @@ SRC = f18estop.f90 \
       stdlib_experimental_quadrature_trapz.f90 \
       stdlib_experimental_stats.f90 \
       stdlib_experimental_stats_mean.f90 \
+      stdlib_experimental_stats_moment.f90 \
       stdlib_experimental_stats_var.f90
 
 LIB = libstdlib.a
@@ -46,6 +47,10 @@ stdlib_experimental_stats_mean.o: \
 	stdlib_experimental_optval.o \
 	stdlib_experimental_kinds.o \
 	stdlib_experimental_stats.o
+stdlib_experimental_stats_moment.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 \
@@ -56,4 +61,5 @@ stdlib_experimental_io.f90: stdlib_experimental_io.fypp
 stdlib_experimental_quadrature.f90: stdlib_experimental_quadrature.fypp
 stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp
 stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp
+stdlib_experimental_stats_moment.f90: stdlib_experimental_stats_moment.fypp
 stdlib_experimental_stats_var.f90: stdlib_experimental_stats_var.fypp

src/stdlib_experimental_stats.fypp

@@ -7,7 +7,7 @@ module stdlib_experimental_stats
   implicit none
   private
   ! Public API
-  public :: mean, var
+  public :: mean, moment, var
 
   interface mean
     #:for k1, t1 in RC_KINDS_TYPES
@@ -209,5 +209,106 @@ module stdlib_experimental_stats
 
   end interface var
 
+  interface moment
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:for rank in RANKS
+        #:set RName = rname("moment_all",rank, t1, k1)
+        module function ${RName}$(x, order, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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("moment_all",rank, t1, k1, 'dp')
+        module function ${RName}$(x, order, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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("moment",rank, t1, k1)
+        module function ${RName}$(x, order, dim, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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("moment",rank, t1, k1, 'dp')
+        module function ${RName}$(x, order, dim, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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("moment_mask_all",rank, t1, k1)
+        module function ${RName}$(x, order, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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("moment_mask_all",rank, t1, k1, 'dp')
+        module function ${RName}$(x, order, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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("moment_mask",rank, t1, k1)
+        module function ${RName}$(x, order, dim, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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("moment_mask",rank, t1, k1, 'dp')
+        module function ${RName}$(x, order, dim, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: order
+          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 moment
 
 end module stdlib_experimental_stats

src/stdlib_experimental_stats.md

@@ -1,9 +1,29 @@
 # Descriptive statistics
 
-## Implemented
 
- * `mean`
- * `var`
+## Implemented
+<!-- vim-markdown-toc GFM -->
+
+* [`mean` - mean of array elements](#mean---mean-of-array-elements)
+	* [Description](#description)
+	* [Syntax](#syntax)
+	* [Arguments](#arguments)
+	* [Return value](#return-value)
+	* [Example](#example)
+* [`moment` - central moment of array elements](#moment---central-moment-of-array-elements)
+	* [Description](#description-1)
+	* [Syntax](#syntax-1)
+	* [Arguments](#arguments-1)
+	* [Return value](#return-value-1)
+	* [Example](#example-1)
+* [`var` - variance of array elements](#var---variance-of-array-elements)
+	* [Description](#description-2)
+	* [Syntax](#syntax-2)
+	* [Arguments](#arguments-2)
+	* [Return value](#return-value-2)
+	* [Example](#example-2)
+
+<!-- vim-markdown-toc -->
 
 ## `mean` - mean of array elements
 
@@ -28,7 +48,7 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a
 ### 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 `integer`, the result is of type `real(dp)`.
 
 If `dim` is absent, a scalar with the mean 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.
 
@@ -49,6 +69,60 @@ program demo_mean
 end program demo_mean
 ```
 
+## `moment` - central moment of array elements
+
+### Description
+
+Returns the _k_-th order central moment 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 _k_-th order central moment is defined as :
+
+```
+ moment(array) = 1/n sum_i (array(i) - mean(array))^k
+```
+
+where n is the number of elements.
+
+### Syntax
+
+`result = moment(array, order [, mask])`
+
+`result = moment(array, order, dim [, mask])`
+
+### Arguments
+
+`array`: Shall be an array of type `integer`, `real`, or `complex`.
+
+`order`: Shall be an scalar of type `integer`.
+
+`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 `real(dp)`.
+
+If `dim` is absent, a scalar with the _k_-th central moment 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 _k_-th central moment 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_moment
+    use stdlib_experimental_stats, only: moment
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    print *, moment(x, 2)                            !returns 2.9167
+    print *, moment( reshape(x, [ 2, 3 ] ), 2)       !returns 2.9167
+    print *, moment( reshape(x, [ 2, 3 ] ), 2, 1)    !returns [0.25, 0.25, 0.25]
+    print *, moment( reshape(x, [ 2, 3 ] ), 2, 1,&
+                     reshape(x, [ 2, 3 ] ) > 3.)     !returns [NaN, 0., 0.25]
+end program demo_moment
+```
+
 ## `var` - variance of array elements
 
 ### Description
@@ -58,7 +132,7 @@ Returns the variance of all the elements of `array`, or of the elements of `arra
 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
+ var(array) = 1/(n-1) sum_i (array(i) - mean(array))^2
 ```
 
 where n is the number of elements.
@@ -108,7 +182,7 @@ program demo_var
                   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]
+                  corrected=.false.)           !returns [NaN, 0., 0.25]
 end program demo_var
 ```