Skip to content

Commit 639cf2f

Browse files
jvdp1jvdp1
authored
Merge pull request #172 from jvdp1/covariance_dev
Covariance of array elements
2 parents 6a3515d + 8f626cd commit 639cf2f

6 files changed

+1024
-12
lines changed

‎src/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ set(fppFiles
55
stdlib_experimental_io.fypp
66
stdlib_experimental_optval.fypp
77
stdlib_experimental_stats.fypp
8+
stdlib_experimental_stats_cov.fypp
89
stdlib_experimental_stats_mean.fypp
910
stdlib_experimental_stats_moment.fypp
1011
stdlib_experimental_stats_var.fypp

‎src/stdlib_experimental_stats.fypp

+97-1
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,102 @@ module stdlib_experimental_stats
88
implicit none
99
private
1010
! Public API
11-
public :: mean, moment, var
11+
public :: cov, mean, moment, var
12+
13+
interface cov
14+
#:for k1, t1 in RC_KINDS_TYPES
15+
#:set RName = rname("cov",1, t1, k1)
16+
module function ${RName}$(x, dim, mask, corrected) result(res)
17+
${t1}$, intent(in) :: x(:)
18+
integer, intent(in) :: dim
19+
logical, intent(in), optional :: mask
20+
logical, intent(in), optional :: corrected
21+
real(${k1}$) :: res
22+
end function ${RName}$
23+
#:endfor
24+
25+
#:for k1, t1 in INT_KINDS_TYPES
26+
#:set RName = rname("cov",1, t1, k1, 'dp')
27+
module function ${RName}$(x, dim, mask, corrected) result(res)
28+
${t1}$, intent(in) :: x(:)
29+
integer, intent(in) :: dim
30+
logical, intent(in), optional :: mask
31+
logical, intent(in), optional :: corrected
32+
real(dp) :: res
33+
end function ${RName}$
34+
#:endfor
35+
36+
#:for k1, t1 in RC_KINDS_TYPES
37+
#:set RName = rname("cov_mask",1, t1, k1)
38+
module function ${RName}$(x, dim, mask, corrected) result(res)
39+
${t1}$, intent(in) :: x(:)
40+
integer, intent(in) :: dim
41+
logical, intent(in) :: mask(:)
42+
logical, intent(in), optional :: corrected
43+
real(${k1}$) :: res
44+
end function ${RName}$
45+
#:endfor
46+
47+
#:for k1, t1 in INT_KINDS_TYPES
48+
#:set RName = rname("cov_mask",1, t1, k1, 'dp')
49+
module function ${RName}$(x, dim, mask, corrected) result(res)
50+
${t1}$, intent(in) :: x(:)
51+
integer, intent(in) :: dim
52+
logical, intent(in) :: mask(:)
53+
logical, intent(in), optional :: corrected
54+
real(dp) :: res
55+
end function ${RName}$
56+
#:endfor
57+
58+
#:for k1, t1 in RC_KINDS_TYPES
59+
#:set RName = rname("cov",2, t1, k1)
60+
module function ${RName}$(x, dim, mask, corrected) result(res)
61+
${t1}$, intent(in) :: x(:, :)
62+
integer, intent(in) :: dim
63+
logical, intent(in), optional :: mask
64+
logical, intent(in), optional :: corrected
65+
${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
66+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
67+
end function ${RName}$
68+
#:endfor
69+
70+
#:for k1, t1 in INT_KINDS_TYPES
71+
#:set RName = rname("cov",2, t1, k1, 'dp')
72+
module function ${RName}$(x, dim, mask, corrected) result(res)
73+
${t1}$, intent(in) :: x(:, :)
74+
integer, intent(in) :: dim
75+
logical, intent(in), optional :: mask
76+
logical, intent(in), optional :: corrected
77+
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
78+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
79+
end function ${RName}$
80+
#:endfor
81+
82+
#:for k1, t1 in RC_KINDS_TYPES
83+
#:set RName = rname("cov_mask",2, t1, k1)
84+
module function ${RName}$(x, dim, mask, corrected) result(res)
85+
${t1}$, intent(in) :: x(:, :)
86+
integer, intent(in) :: dim
87+
logical, intent(in) :: mask(:,:)
88+
logical, intent(in), optional :: corrected
89+
${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
90+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
91+
end function ${RName}$
92+
#:endfor
93+
94+
#:for k1, t1 in INT_KINDS_TYPES
95+
#:set RName = rname("cov_mask",2, t1, k1, 'dp')
96+
module function ${RName}$(x, dim, mask, corrected) result(res)
97+
${t1}$, intent(in) :: x(:, :)
98+
integer, intent(in) :: dim
99+
logical, intent(in) :: mask(:,:)
100+
logical, intent(in), optional :: corrected
101+
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
102+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
103+
end function ${RName}$
104+
#:endfor
105+
end interface cov
106+
12107

13108
interface mean
14109
#:for k1, t1 in RC_KINDS_TYPES
@@ -210,6 +305,7 @@ module stdlib_experimental_stats
210305

211306
end interface var
212307

308+
213309
interface moment
214310
#:for k1, t1 in RC_KINDS_TYPES
215311
#:for rank in RANKS

‎src/stdlib_experimental_stats.md

+69-11
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,67 @@
11
# Descriptive statistics
22

3+
* [`cov` - covariance of array elements](#cov---covariance-of-array-elements)
34
* [`mean` - mean of array elements](#mean---mean-of-array-elements)
45
* [`moment` - central moments of array elements](#moment---central-moments-of-array-elements)
56
* [`var` - variance of array elements](#var---variance-of-array-elements)
67

8+
9+
## `cov` - covariance of array elements
10+
11+
### Description
12+
13+
Returns the covariance of the elements of `array` along dimension `dim` if the corresponding element in `mask` is `true`.
14+
15+
Per default, the covariance is defined as:
16+
17+
```
18+
cov(array) = 1/(n-1) sum_i (array(i) - mean(array) * (array(i) - mean(array)))
19+
```
20+
21+
where `n` is the number of elements.
22+
23+
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`.
24+
25+
26+
### Syntax
27+
28+
`result = cov(array, dim [, mask [, corrected]])`
29+
30+
### Arguments
31+
32+
`array`: Shall be a rank-1 or a rank-2 array of type `integer`, `real`, or `complex`.
33+
34+
`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`.
35+
36+
`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
37+
38+
`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`.
39+
40+
### Return value
41+
42+
If `array` is of rank 1 and of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
43+
If `array` is of rank 2 and of type `real` or `complex`, the result is of the same type as `array`.
44+
If `array` is of type `integer`, the result is of type `real(dp)`.
45+
46+
If `array` is of rank 1, a scalar with the covariance (that is the variance) of all elements in `array` is returned.
47+
If `array` is of rank 2, a rank-2 array is returned.
48+
49+
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`.
50+
51+
### Example
52+
53+
```fortran
54+
program demo_cov
55+
use stdlib_experimental_stats, only: cov
56+
implicit none
57+
real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
58+
real :: y(1:2, 1:3) = reshape([ 1., 2., 3., 4., 5., 6. ], [ 2, 3])
59+
print *, cov(x, 1) !returns 3.5
60+
print *, cov(x, 1, corrected = .false.) !returns 2.9167
61+
print *, cov(y, 1) !returns a square matrix of size 3 with all elements equal to 0.5
62+
end program demo_cov
63+
```
64+
765
## `mean` - mean of array elements
866

967
### Description
@@ -20,16 +78,16 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a
2078

2179
`array`: Shall be an array of type `integer`, `real`, or `complex`.
2280

23-
`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`.
81+
`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`.
2482

25-
`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
83+
`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
2684

2785
### Return value
2886

2987
If `array` is of type `real` or `complex`, the result is of the same type as `array`.
3088
If `array` is of type `integer`, the result is of type `real(dp)`.
3189

32-
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.
90+
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.
3391

3492
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`.
3593

@@ -63,7 +121,7 @@ The _k_-th order central moment is defined as :
63121
moment(array) = 1/n sum_i (array(i) - mean(array))^k
64122
```
65123

66-
where n is the number of elements.
124+
where `n` is the number of elements.
67125

68126
The _k_-th order moment about `center` is defined as :
69127

@@ -83,18 +141,18 @@ The _k_-th order moment about `center` is defined as :
83141

84142
`order`: Shall be an scalar of type `integer`.
85143

86-
`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`.
144+
`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`.
87145

88146
`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`.
89147

90-
`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
148+
`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
91149

92150
### Return value
93151

94152
If `array` is of type `real` or `complex`, the result is of the same type as `array`.
95153
If `array` is of type `integer`, the result is of type `real(dp)`.
96154

97-
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.
155+
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.
98156

99157
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`.
100158

@@ -127,7 +185,7 @@ Per default, the variance is defined as the best unbiased estimator and is compu
127185
var(array) = 1/(n-1) sum_i (array(i) - mean(array))^2
128186
```
129187

130-
where n is the number of elements.
188+
where `n` is the number of elements.
131189

132190
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`.
133191

@@ -142,9 +200,9 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin
142200

143201
`array`: Shall be an array of type `integer`, `real`, or `complex`.
144202

145-
`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`.
203+
`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`.
146204

147-
`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
205+
`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
148206

149207
`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`.
150208

@@ -153,7 +211,7 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin
153211
If `array` is of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
154212
If `array` is of type `integer`, the result is of type `real(dp)`.
155213

156-
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.
214+
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.
157215

158216
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`.
159217

0 commit comments

Comments
 (0)