Covariance of array elements

src/CMakeLists.txt

@@ -5,6 +5,7 @@ set(fppFiles
     stdlib_experimental_io.fypp
     stdlib_experimental_optval.fypp
     stdlib_experimental_stats.fypp
+    stdlib_experimental_stats_cov.fypp
     stdlib_experimental_stats_mean.fypp
     stdlib_experimental_stats_moment.fypp
     stdlib_experimental_stats_var.fypp

src/stdlib_experimental_stats.fypp

@@ -8,7 +8,102 @@ module stdlib_experimental_stats
   implicit none
   private
   ! Public API
-  public :: mean, moment, var
+  public :: cov, mean, moment, var
+
+  interface cov
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("cov",1, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        real(${k1}$) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("cov",1, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        real(dp) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("cov_mask",1, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:)
+        logical, intent(in), optional :: corrected
+        real(${k1}$) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("cov_mask",1, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:)
+        logical, intent(in), optional :: corrected
+        real(dp) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("cov",2, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                            , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("cov",2, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                        , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("cov_mask",2, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:,:)
+        logical, intent(in), optional :: corrected
+        ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                            , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("cov_mask",2, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:,:)
+        logical, intent(in), optional :: corrected
+        real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                            , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+  end interface cov
+
 
   interface mean
     #:for k1, t1 in RC_KINDS_TYPES
@@ -210,6 +305,7 @@ module stdlib_experimental_stats
 
   end interface var
 
+
   interface moment
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS

src/stdlib_experimental_stats.md

@@ -1,9 +1,67 @@
 # Descriptive statistics
 
+* [`cov` - covariance of array elements](#cov---covariance-of-array-elements)
 * [`mean` - mean of array elements](#mean---mean-of-array-elements)
 * [`moment` - central moments of array elements](#moment---central-moments-of-array-elements)
 * [`var` - variance of array elements](#var---variance-of-array-elements)
 
+
+## `cov` - covariance of array elements
+
+### Description
+
+Returns the covariance of the elements of `array` along dimension `dim` if the corresponding element in `mask` is `true`.
+
+Per default, the covariance is defined as:
+
+```
+ cov(array) = 1/(n-1) sum_i (array(i) - mean(array) * (array(i) - mean(array)))
+```
+
+where `n` is the number of elements.
+
+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 = cov(array, dim [, mask [, corrected]])`
+
+### Arguments
+
+`array`: Shall be a rank-1 or a rank-2 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 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 rank 1 and of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
+If `array` is of rank 2 and 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 `array` is of rank 1, a scalar with the covariance (that is the variance) of all elements in `array` is returned.
+If `array` is of rank 2, a rank-2 array is returned.
+
+If `mask` is specified, the result is the covariance 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_cov
+    use stdlib_experimental_stats, only: cov
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    real :: y(1:2, 1:3) = reshape([ 1., 2., 3., 4., 5., 6. ], [ 2, 3])
+    print *, cov(x, 1)                          !returns 3.5
+    print *, cov(x, 1, corrected = .false.)     !returns 2.9167
+    print *, cov(y, 1)                          !returns a square matrix of size 3 with all elements equal to 0.5
+end program demo_cov
+```
+
 ## `mean` - mean of array elements
 
 ### Description
@@ -20,16 +78,16 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a
 
 `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`.
+`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`.
+`mask` (optional): Shall be of type `logical` and either 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 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.
+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.
 
 If `mask` is specified, the result is the mean of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`.
 
@@ -63,7 +121,7 @@ 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.
+where `n` is the number of elements.
 
 The _k_-th order moment about `center` is defined as :
 
@@ -83,18 +141,18 @@ The _k_-th order moment about `center` is defined as :
 
 `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`.
+`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`.
 
 `center` (optional): Shall be a scalar of the same type of `result` if `dim` is not provided. If `dim` is provided, `center` shall be a scalar or an array (with a shape similar to that of `array` with dimension `dim` dropped) of the same type of `result`.
 
-`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
+`mask` (optional): Shall be of type `logical` and either 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 `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`.
 
@@ -127,7 +185,7 @@ Per default, the variance is defined as the best unbiased estimator and is compu
  var(array) = 1/(n-1) sum_i (array(i) - mean(array))^2
 ```
 
-where n is the number of elements.
+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`.
 
@@ -142,9 +200,9 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin
 
 `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`.
+`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`.
+`mask` (optional): Shall be of type `logical` and either 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`.
 
@@ -153,7 +211,7 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin
 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 `array` 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`.