Skip to content

Commit 50367ff

Browse files
jvdp1jvdp1
committed
Addition of corrected in API of var
1 parent 1e0e460 commit 50367ff

File tree

6 files changed

+237
-86
lines changed

6 files changed

+237
-86
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: 35 additions & 18 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
@@ -30,9 +31,11 @@ contains
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._${k1}$, 0._${k1}$,&
35+
(optval(corrected, .true.) .and. n > 0._${k1}$)))
3436
#:else
35-
res = sum(abs(x - mean)**2) / (n - 1._${k1}$)
37+
res = sum(abs(x - mean)**2) / (n - merge(1._${k1}$, 0._${k1}$,&
38+
(optval(corrected, .true.) .and. n > 0._${k1}$)))
3639
#:endif
3740

3841
end function ${RName}$
@@ -43,9 +46,10 @@ contains
4346
#:for k1, t1 in INT_KINDS_TYPES
4447
#:for rank in RANKS
4548
#:set RName = rname("var_all",rank, t1, k1, 'dp')
46-
module function ${RName}$(x, mask) result(res)
49+
module function ${RName}$(x, mask, corrected) result(res)
4750
${t1}$, intent(in) :: x${ranksuffix(rank)}$
4851
logical, intent(in), optional :: mask
52+
logical, intent(in), optional :: corrected
4953
real(dp) :: res
5054

5155
real(dp) :: n, mean
@@ -58,7 +62,8 @@ contains
5862
n = real(size(x, kind = int64), dp)
5963
mean = sum(real(x, dp)) / n
6064

61-
res = sum((real(x, dp) - mean)**2) / (n - 1._dp)
65+
res = sum((real(x, dp) - mean)**2) / (n - merge(1._dp, 0._dp,&
66+
(optval(corrected, .true.) .and. n > 0._dp)))
6267

6368
end function ${RName}$
6469
#:endfor
@@ -68,10 +73,11 @@ contains
6873
#:for k1, t1 in RC_KINDS_TYPES
6974
#:for rank in RANKS
7075
#:set RName = rname("var",rank, t1, k1)
71-
module function ${RName}$(x, dim, mask) result(res)
76+
module function ${RName}$(x, dim, mask, corrected) result(res)
7277
${t1}$, intent(in) :: x${ranksuffix(rank)}$
7378
integer, intent(in) :: dim
7479
logical, intent(in), optional :: mask
80+
logical, intent(in), optional :: corrected
7581
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
7682

7783
integer :: i
@@ -100,7 +106,8 @@ contains
100106
case default
101107
call error_stop("ERROR (mean): wrong dimension")
102108
end select
103-
res = res / (n - 1._${k1}$)
109+
res = res / (n - merge(1._${k1}$, 0._${k1}$,&
110+
(optval(corrected, .true.) .and. n > 0._${k1}$)))
104111

105112
end function ${RName}$
106113
#:endfor
@@ -110,10 +117,11 @@ contains
110117
#:for k1, t1 in INT_KINDS_TYPES
111118
#:for rank in RANKS
112119
#:set RName = rname("var",rank, t1, k1, 'dp')
113-
module function ${RName}$(x, dim, mask) result(res)
120+
module function ${RName}$(x, dim, mask, corrected) result(res)
114121
${t1}$, intent(in) :: x${ranksuffix(rank)}$
115122
integer, intent(in) :: dim
116123
logical, intent(in), optional :: mask
124+
logical, intent(in), optional :: corrected
117125
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
118126

119127
integer :: i
@@ -138,7 +146,8 @@ contains
138146
case default
139147
call error_stop("ERROR (mean): wrong dimension")
140148
end select
141-
res = res / (n - 1._dp)
149+
res = res / (n - merge(1._dp, 0._dp,&
150+
(optval(corrected, .true.) .and. n > 0._dp)))
142151

143152
end function ${RName}$
144153
#:endfor
@@ -148,9 +157,10 @@ contains
148157
#:for k1, t1 in RC_KINDS_TYPES
149158
#:for rank in RANKS
150159
#:set RName = rname("var_mask_all",rank, t1, k1)
151-
module function ${RName}$(x, mask) result(res)
160+
module function ${RName}$(x, mask, corrected) result(res)
152161
${t1}$, intent(in) :: x${ranksuffix(rank)}$
153162
logical, intent(in) :: mask${ranksuffix(rank)}$
163+
logical, intent(in), optional :: corrected
154164
real(${k1}$) :: res
155165

156166
real(${k1}$) :: n
@@ -160,9 +170,11 @@ contains
160170
mean = sum(x, mask) / n
161171

162172
#:if t1[0] == 'r'
163-
res = sum((x - mean)**2, mask) / (n - 1._${k1}$)
173+
res = sum((x - mean)**2, mask) / (n -&
174+
merge(1._${k1}$, 0._${k1}$, (optval(corrected, .true.)) .and. n > 0._${k1}$))
164175
#:else
165-
res = sum(abs(x - mean)**2, mask) / (n - 1._${k1}$)
176+
res = sum(abs(x - mean)**2, mask) / (n - merge(1._${k1}$, 0._${k1}$,&
177+
(optval(corrected, .true.) .and. n > 0._${k1}$)))
166178
#:endif
167179

168180
end function ${RName}$
@@ -173,17 +185,19 @@ contains
173185
#:for k1, t1 in INT_KINDS_TYPES
174186
#:for rank in RANKS
175187
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
176-
module function ${RName}$(x, mask) result(res)
188+
module function ${RName}$(x, mask, corrected) result(res)
177189
${t1}$, intent(in) :: x${ranksuffix(rank)}$
178190
logical, intent(in) :: mask${ranksuffix(rank)}$
191+
logical, intent(in), optional :: corrected
179192
real(dp) :: res
180193

181194
real(dp) :: n, mean
182195

183196
n = real(count(mask, kind = int64), dp)
184197
mean = sum(real(x, dp), mask) / n
185198

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

188202
end function ${RName}$
189203
#:endfor
@@ -193,10 +207,11 @@ contains
193207
#:for k1, t1 in RC_KINDS_TYPES
194208
#:for rank in RANKS
195209
#:set RName = rname("var_mask",rank, t1, k1)
196-
module function ${RName}$(x, dim, mask) result(res)
210+
module function ${RName}$(x, dim, mask, corrected) result(res)
197211
${t1}$, intent(in) :: x${ranksuffix(rank)}$
198212
integer, intent(in) :: dim
199213
logical, intent(in) :: mask${ranksuffix(rank)}$
214+
logical, intent(in), optional :: corrected
200215
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
201216

202217
integer :: i
@@ -222,7 +237,8 @@ contains
222237
case default
223238
call error_stop("ERROR (mean): wrong dimension")
224239
end select
225-
res = res / (n - 1._${k1}$)
240+
res = res / (n - merge(1._${k1}$, 0._${k1}$,&
241+
(optval(corrected, .true.) .and. n > 0._${k1}$)))
226242

227243
end function ${RName}$
228244
#:endfor
@@ -232,10 +248,11 @@ contains
232248
#:for k1, t1 in INT_KINDS_TYPES
233249
#:for rank in RANKS
234250
#:set RName = rname("var_mask",rank, t1, k1, 'dp')
235-
module function ${RName}$(x, dim, mask) result(res)
251+
module function ${RName}$(x, dim, mask, corrected) result(res)
236252
${t1}$, intent(in) :: x${ranksuffix(rank)}$
237253
integer, intent(in) :: dim
238254
logical, intent(in) :: mask${ranksuffix(rank)}$
255+
logical, intent(in), optional :: corrected
239256
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
240257

241258
integer :: i
@@ -256,7 +273,7 @@ contains
256273
case default
257274
call error_stop("ERROR (mean): wrong dimension")
258275
end select
259-
res = res / (n - 1._dp)
276+
res = res / (n - merge(1._dp, 0._dp, (optval(corrected, .true.) .and. n > 0._dp)))
260277

261278
end function ${RName}$
262279
#:endfor

‎src/tests/stats/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
ADDTEST(mean)
22
ADDTEST(var)
3+
ADDTEST(varn)
34

45
if(DEFINED CMAKE_MAXIMUM_RANK)
56
if(${CMAKE_MAXIMUM_RANK} GREATER 7)

0 commit comments

Comments
 (0)