Skip to content

Commit 5f35833

Browse files
milancurcicmilancurcic
authored
Merge pull request #151 from jvdp1/var_corrected_dev
Addition of corrected in API of var (following #149)
2 parents 1e0e460 + 9af45ef commit 5f35833

File tree

6 files changed

+467
-94
lines changed

6 files changed

+467
-94
lines changed

‎src/stdlib_experimental_stats.fypp

Lines changed: 16 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -110,9 +110,10 @@ module stdlib_experimental_stats
110110
#:for k1, t1 in RC_KINDS_TYPES
111111
#:for rank in RANKS
112112
#:set RName = rname("var_all",rank, t1, k1)
113-
module function ${RName}$(x, mask) result(res)
113+
module function ${RName}$(x, mask, corrected) result(res)
114114
${t1}$, intent(in) :: x${ranksuffix(rank)}$
115115
logical, intent(in), optional :: mask
116+
logical, intent(in), optional :: corrected
116117
real(${k1}$) :: res
117118
end function ${RName}$
118119
#:endfor
@@ -121,9 +122,10 @@ module stdlib_experimental_stats
121122
#:for k1, t1 in INT_KINDS_TYPES
122123
#:for rank in RANKS
123124
#:set RName = rname("var_all",rank, t1, k1, 'dp')
124-
module function ${RName}$(x, mask) result(res)
125+
module function ${RName}$(x, mask, corrected) result(res)
125126
${t1}$, intent(in) :: x${ranksuffix(rank)}$
126127
logical, intent(in), optional :: mask
128+
logical, intent(in), optional :: corrected
127129
real(dp) :: res
128130
end function ${RName}$
129131
#:endfor
@@ -132,10 +134,11 @@ module stdlib_experimental_stats
132134
#:for k1, t1 in RC_KINDS_TYPES
133135
#:for rank in RANKS
134136
#:set RName = rname("var",rank, t1, k1)
135-
module function ${RName}$(x, dim, mask) result(res)
137+
module function ${RName}$(x, dim, mask, corrected) result(res)
136138
${t1}$, intent(in) :: x${ranksuffix(rank)}$
137139
integer, intent(in) :: dim
138140
logical, intent(in), optional :: mask
141+
logical, intent(in), optional :: corrected
139142
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
140143
end function ${RName}$
141144
#:endfor
@@ -144,10 +147,11 @@ module stdlib_experimental_stats
144147
#:for k1, t1 in INT_KINDS_TYPES
145148
#:for rank in RANKS
146149
#:set RName = rname("var",rank, t1, k1, 'dp')
147-
module function ${RName}$(x, dim, mask) result(res)
150+
module function ${RName}$(x, dim, mask, corrected) result(res)
148151
${t1}$, intent(in) :: x${ranksuffix(rank)}$
149152
integer, intent(in) :: dim
150153
logical, intent(in), optional :: mask
154+
logical, intent(in), optional :: corrected
151155
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
152156
end function ${RName}$
153157
#:endfor
@@ -156,9 +160,10 @@ module stdlib_experimental_stats
156160
#:for k1, t1 in RC_KINDS_TYPES
157161
#:for rank in RANKS
158162
#:set RName = rname("var_mask_all",rank, t1, k1)
159-
module function ${RName}$(x, mask) result(res)
163+
module function ${RName}$(x, mask, corrected) result(res)
160164
${t1}$, intent(in) :: x${ranksuffix(rank)}$
161165
logical, intent(in) :: mask${ranksuffix(rank)}$
166+
logical, intent(in), optional :: corrected
162167
real(${k1}$) :: res
163168
end function ${RName}$
164169
#:endfor
@@ -167,9 +172,10 @@ module stdlib_experimental_stats
167172
#:for k1, t1 in INT_KINDS_TYPES
168173
#:for rank in RANKS
169174
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
170-
module function ${RName}$(x, mask) result(res)
175+
module function ${RName}$(x, mask, corrected) result(res)
171176
${t1}$, intent(in) :: x${ranksuffix(rank)}$
172177
logical, intent(in) :: mask${ranksuffix(rank)}$
178+
logical, intent(in), optional :: corrected
173179
real(dp) :: res
174180
end function ${RName}$
175181
#:endfor
@@ -178,10 +184,11 @@ module stdlib_experimental_stats
178184
#:for k1, t1 in RC_KINDS_TYPES
179185
#:for rank in RANKS
180186
#:set RName = rname("var_mask",rank, t1, k1)
181-
module function ${RName}$(x, dim, mask) result(res)
187+
module function ${RName}$(x, dim, mask, corrected) result(res)
182188
${t1}$, intent(in) :: x${ranksuffix(rank)}$
183189
integer, intent(in) :: dim
184190
logical, intent(in) :: mask${ranksuffix(rank)}$
191+
logical, intent(in), optional :: corrected
185192
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
186193
end function ${RName}$
187194
#:endfor
@@ -190,10 +197,11 @@ module stdlib_experimental_stats
190197
#:for k1, t1 in INT_KINDS_TYPES
191198
#:for rank in RANKS
192199
#:set RName = rname("var_mask",rank, t1, k1, 'dp')
193-
module function ${RName}$(x, dim, mask) result(res)
200+
module function ${RName}$(x, dim, mask, corrected) result(res)
194201
${t1}$, intent(in) :: x${ranksuffix(rank)}$
195202
integer, intent(in) :: dim
196203
logical, intent(in) :: mask${ranksuffix(rank)}$
204+
logical, intent(in), optional :: corrected
197205
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
198206
end function ${RName}$
199207
#:endfor

‎src/stdlib_experimental_stats.md

Lines changed: 20 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -55,17 +55,22 @@ end program demo_mean
5555

5656
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`.
5757

58-
The variance is defined as the best unbiased estimator and is computed as:
58+
Per default, the variance is defined as the best unbiased estimator and is computed as:
5959

6060
```
6161
var(x) = 1/(n-1) sum_i (array(i) - mean(array))^2
6262
```
6363

64+
where n is the number of elements.
65+
66+
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`.
67+
68+
6469
### Syntax
6570

66-
`result = var(array [, mask])`
71+
`result = var(array [, mask [, corrected]])`
6772

68-
`result = var(array, dim [, mask])`
73+
`result = var(array, dim [, mask [, corrected]])`
6974

7075
### Arguments
7176

@@ -75,16 +80,19 @@ The variance is defined as the best unbiased estimator and is computed as:
7580

7681
`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
7782

83+
`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`.
84+
7885
### Return value
7986

80-
If `array` is of type `real` or `complex`, the result is of the same type as `array`.
81-
If `array` is of type `integer`, the result is of type `double precision`.
87+
If `array` is of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
88+
If `array` is of type `integer`, the result is of type `real(dp)`.
8289

83-
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
84-
ray` with dimension `dim` dropped is returned.
90+
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.
8591

8692
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`.
8793

94+
If the variance is computed with only one single element, then the result is IEEE `NaN` if `corrected` is `.true.` and is `0.` if `corrected` is `.false.`.
95+
8896
### Example
8997

9098
```fortran
@@ -93,10 +101,14 @@ program demo_var
93101
implicit none
94102
real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
95103
print *, var(x) !returns 3.5
104+
print *, var(x, corrected = .false.) !returns 2.9167
96105
print *, var( reshape(x, [ 2, 3 ] )) !returns 3.5
97106
print *, var( reshape(x, [ 2, 3 ] ), 1) !returns [0.5, 0.5, 0.5]
98107
print *, var( reshape(x, [ 2, 3 ] ), 1,&
99-
reshape(x, [ 2, 3 ] ) > 3.) !returns [0., NaN, 0.5]
108+
reshape(x, [ 2, 3 ] ) > 3.) !returns [NaN, NaN, 0.5]
109+
print *, var( reshape(x, [ 2, 3 ] ), 1,&
110+
reshape(x, [ 2, 3 ] ) > 3.,&
111+
corrected=.false.) !returns [NaN, 0., 0.5]
100112
end program demo_var
101113
```
102114

‎src/stdlib_experimental_stats_var.fypp

Lines changed: 36 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,10 @@ contains
1313
#:for k1, t1 in RC_KINDS_TYPES
1414
#:for rank in RANKS
1515
#:set RName = rname("var_all",rank, t1, k1)
16-
module function ${RName}$(x, mask) result(res)
16+
module function ${RName}$(x, mask, corrected) result(res)
1717
${t1}$, intent(in) :: x${ranksuffix(rank)}$
1818
logical, intent(in), optional :: mask
19+
logical, intent(in), optional :: corrected
1920
real(${k1}$) :: res
2021

2122
real(${k1}$) :: n
@@ -26,13 +27,13 @@ contains
2627
return
2728
end if
2829

29-
n = real(size(x, kind = int64), ${k1}$)
30+
n = size(x, kind = int64)
3031
mean = sum(x) / n
3132

3233
#:if t1[0] == 'r'
33-
res = sum((x - mean)**2) / (n - 1._${k1}$)
34+
res = sum((x - mean)**2) / (n - merge(1, 0 , optval(corrected, .true.)))
3435
#:else
35-
res = sum(abs(x - mean)**2) / (n - 1._${k1}$)
36+
res = sum(abs(x - mean)**2) / (n - merge(1, 0, optval(corrected, .true.)))
3637
#:endif
3738

3839
end function ${RName}$
@@ -43,9 +44,10 @@ contains
4344
#:for k1, t1 in INT_KINDS_TYPES
4445
#:for rank in RANKS
4546
#:set RName = rname("var_all",rank, t1, k1, 'dp')
46-
module function ${RName}$(x, mask) result(res)
47+
module function ${RName}$(x, mask, corrected) result(res)
4748
${t1}$, intent(in) :: x${ranksuffix(rank)}$
4849
logical, intent(in), optional :: mask
50+
logical, intent(in), optional :: corrected
4951
real(dp) :: res
5052

5153
real(dp) :: n, mean
@@ -55,10 +57,10 @@ contains
5557
return
5658
end if
5759

58-
n = real(size(x, kind = int64), dp)
60+
n = size(x, kind = int64)
5961
mean = sum(real(x, dp)) / n
6062

61-
res = sum((real(x, dp) - mean)**2) / (n - 1._dp)
63+
res = sum((real(x, dp) - mean)**2) / (n - merge(1, 0, optval(corrected, .true.)))
6264

6365
end function ${RName}$
6466
#:endfor
@@ -68,10 +70,11 @@ contains
6870
#:for k1, t1 in RC_KINDS_TYPES
6971
#:for rank in RANKS
7072
#:set RName = rname("var",rank, t1, k1)
71-
module function ${RName}$(x, dim, mask) result(res)
73+
module function ${RName}$(x, dim, mask, corrected) result(res)
7274
${t1}$, intent(in) :: x${ranksuffix(rank)}$
7375
integer, intent(in) :: dim
7476
logical, intent(in), optional :: mask
77+
logical, intent(in), optional :: corrected
7578
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
7679

7780
integer :: i
@@ -87,7 +90,7 @@ contains
8790
select case(dim)
8891
#:for fi in range(1, rank+1)
8992
case(${fi}$)
90-
n = real(size(x, dim), ${k1}$)
93+
n = size(x, dim)
9194
mean = sum(x, dim) / n
9295
do i = 1, size(x, dim)
9396
#:if t1[0] == 'r'
@@ -100,7 +103,7 @@ contains
100103
case default
101104
call error_stop("ERROR (mean): wrong dimension")
102105
end select
103-
res = res / (n - 1._${k1}$)
106+
res = res / (n - merge(1, 0, optval(corrected, .true.)))
104107

105108
end function ${RName}$
106109
#:endfor
@@ -110,10 +113,11 @@ contains
110113
#:for k1, t1 in INT_KINDS_TYPES
111114
#:for rank in RANKS
112115
#:set RName = rname("var",rank, t1, k1, 'dp')
113-
module function ${RName}$(x, dim, mask) result(res)
116+
module function ${RName}$(x, dim, mask, corrected) result(res)
114117
${t1}$, intent(in) :: x${ranksuffix(rank)}$
115118
integer, intent(in) :: dim
116119
logical, intent(in), optional :: mask
120+
logical, intent(in), optional :: corrected
117121
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
118122

119123
integer :: i
@@ -129,7 +133,7 @@ contains
129133
select case(dim)
130134
#:for fi in range(1, rank+1)
131135
case(${fi}$)
132-
n = real(size(x, dim), dp)
136+
n = size(x, dim)
133137
mean = sum(real(x, dp), dim) / n
134138
do i = 1, size(x, dim)
135139
res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2
@@ -138,7 +142,7 @@ contains
138142
case default
139143
call error_stop("ERROR (mean): wrong dimension")
140144
end select
141-
res = res / (n - 1._dp)
145+
res = res / (n - merge(1, 0, optval(corrected, .true.)))
142146

143147
end function ${RName}$
144148
#:endfor
@@ -148,22 +152,24 @@ contains
148152
#:for k1, t1 in RC_KINDS_TYPES
149153
#:for rank in RANKS
150154
#:set RName = rname("var_mask_all",rank, t1, k1)
151-
module function ${RName}$(x, mask) result(res)
155+
module function ${RName}$(x, mask, corrected) result(res)
152156
${t1}$, intent(in) :: x${ranksuffix(rank)}$
153157
logical, intent(in) :: mask${ranksuffix(rank)}$
158+
logical, intent(in), optional :: corrected
154159
real(${k1}$) :: res
155160

156161
real(${k1}$) :: n
157162
${t1}$ :: mean
158163

159-
n = real(count(mask, kind = int64), ${k1}$)
164+
n = count(mask, kind = int64)
160165
mean = sum(x, mask) / n
161166

162167
#:if t1[0] == 'r'
163-
res = sum((x - mean)**2, mask) / (n - 1._${k1}$)
168+
res = sum((x - mean)**2, mask) / (n -&
164169
#:else
165-
res = sum(abs(x - mean)**2, mask) / (n - 1._${k1}$)
170+
res = sum(abs(x - mean)**2, mask) / (n -&
166171
#:endif
172+
merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
167173

168174
end function ${RName}$
169175
#:endfor
@@ -173,17 +179,19 @@ contains
173179
#:for k1, t1 in INT_KINDS_TYPES
174180
#:for rank in RANKS
175181
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
176-
module function ${RName}$(x, mask) result(res)
182+
module function ${RName}$(x, mask, corrected) result(res)
177183
${t1}$, intent(in) :: x${ranksuffix(rank)}$
178184
logical, intent(in) :: mask${ranksuffix(rank)}$
185+
logical, intent(in), optional :: corrected
179186
real(dp) :: res
180187

181188
real(dp) :: n, mean
182189

183-
n = real(count(mask, kind = int64), dp)
190+
n = count(mask, kind = int64)
184191
mean = sum(real(x, dp), mask) / n
185192

186-
res = sum((real(x, dp) - mean)**2, mask) / (n - 1._dp)
193+
res = sum((real(x, dp) - mean)**2, mask) / (n -&
194+
merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
187195

188196
end function ${RName}$
189197
#:endfor
@@ -193,10 +201,11 @@ contains
193201
#:for k1, t1 in RC_KINDS_TYPES
194202
#:for rank in RANKS
195203
#:set RName = rname("var_mask",rank, t1, k1)
196-
module function ${RName}$(x, dim, mask) result(res)
204+
module function ${RName}$(x, dim, mask, corrected) result(res)
197205
${t1}$, intent(in) :: x${ranksuffix(rank)}$
198206
integer, intent(in) :: dim
199207
logical, intent(in) :: mask${ranksuffix(rank)}$
208+
logical, intent(in), optional :: corrected
200209
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
201210

202211
integer :: i
@@ -207,7 +216,7 @@ contains
207216
select case(dim)
208217
#:for fi in range(1, rank+1)
209218
case(${fi}$)
210-
n = real(count(mask, dim), ${k1}$)
219+
n = count(mask, dim)
211220
mean = sum(x, dim, mask) / n
212221
do i = 1, size(x, dim)
213222
#:if t1[0] == 'r'
@@ -222,7 +231,7 @@ contains
222231
case default
223232
call error_stop("ERROR (mean): wrong dimension")
224233
end select
225-
res = res / (n - 1._${k1}$)
234+
res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
226235

227236
end function ${RName}$
228237
#:endfor
@@ -232,10 +241,11 @@ contains
232241
#:for k1, t1 in INT_KINDS_TYPES
233242
#:for rank in RANKS
234243
#:set RName = rname("var_mask",rank, t1, k1, 'dp')
235-
module function ${RName}$(x, dim, mask) result(res)
244+
module function ${RName}$(x, dim, mask, corrected) result(res)
236245
${t1}$, intent(in) :: x${ranksuffix(rank)}$
237246
integer, intent(in) :: dim
238247
logical, intent(in) :: mask${ranksuffix(rank)}$
248+
logical, intent(in), optional :: corrected
239249
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
240250

241251
integer :: i
@@ -246,7 +256,7 @@ contains
246256
select case(dim)
247257
#:for fi in range(1, rank+1)
248258
case(${fi}$)
249-
n = real(count(mask, dim), dp)
259+
n = count(mask, dim)
250260
mean = sum(real(x, dp), dim, mask) / n
251261
do i = 1, size(x, dim)
252262
res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2,&
@@ -256,7 +266,7 @@ contains
256266
case default
257267
call error_stop("ERROR (mean): wrong dimension")
258268
end select
259-
res = res / (n - 1._dp)
269+
res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0)))
260270

261271
end function ${RName}$
262272
#:endfor

0 commit comments

Comments
 (0)