Skip to content

Commit cd2c16c

Browse files
committed
mv_to_central_moment: changed default to central moment
1 parent cb69df8 commit cd2c16c

File tree

3 files changed

+720
-676
lines changed

3 files changed

+720
-676
lines changed

src/stdlib_experimental_stats_moment.fypp

Lines changed: 20 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ contains
3232
if (present(center)) then
3333
res = sum((x - center)**order) / n
3434
else
35-
res = sum((x)**order) / n
35+
res = sum((x - mean(x))**order) / n
3636
end if
3737

3838
end function ${RName}$
@@ -62,7 +62,7 @@ contains
6262
if (present(center)) then
6363
res = sum((real(x, dp) - center)**order) / n
6464
else
65-
res = sum((real(x, dp))**order) / n
65+
res = sum((real(x, dp) - mean(x))**order) / n
6666
end if
6767

6868
end function ${RName}$
@@ -83,6 +83,7 @@ contains
8383

8484
integer :: i
8585
real(${k1}$) :: n
86+
${t1}$, allocatable :: mean_${ranksuffix(rank-1)}$
8687

8788
if (.not.optval(mask, .true.)) then
8889
res = ieee_value(1._${k1}$, ieee_quiet_nan)
@@ -100,9 +101,11 @@ contains
100101
res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - center)**order
101102
end do
102103
else
104+
allocate(mean_, source = mean(x, ${fi}$))
103105
do i = 1, size(x, ${fi}$)
104-
res = res + (x${select_subarray(rank, [(fi, 'i')])}$)**order
106+
res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean_)**order
105107
end do
108+
deallocate(mean_)
106109
end if
107110
#:endfor
108111
case default
@@ -128,6 +131,7 @@ contains
128131

129132
integer :: i
130133
real(dp) :: n
134+
real(dp), allocatable :: mean_${ranksuffix(rank-1)}$
131135

132136
if (.not.optval(mask, .true.)) then
133137
res = ieee_value(1._dp, ieee_quiet_nan)
@@ -146,9 +150,11 @@ contains
146150
center)**order
147151
end do
148152
else
153+
allocate(mean_, source = mean(x, ${fi}$))
149154
do i = 1, size(x, ${fi}$)
150-
res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp))**order
155+
res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean_)**order
151156
end do
157+
deallocate(mean_)
152158
end if
153159
#:endfor
154160
case default
@@ -178,7 +184,7 @@ contains
178184
if (present(center)) then
179185
res = sum((x - center)**order, mask) / n
180186
else
181-
res = sum((x)**order, mask) / n
187+
res = sum((x - mean(x, mask))**order, mask) / n
182188
end if
183189

184190
end function ${RName}$
@@ -203,7 +209,7 @@ contains
203209
if (present(center)) then
204210
res = sum((real(x, dp) - center)**order, mask) / n
205211
else
206-
res = sum((real(x, dp))**order, mask) / n
212+
res = sum((real(x, dp) - mean(x,mask))**order, mask) / n
207213
end if
208214

209215
end function ${RName}$
@@ -224,6 +230,7 @@ contains
224230

225231
integer :: i
226232
real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
233+
${t1}$, allocatable :: mean_${ranksuffix(rank-1)}$
227234

228235
n = count(mask, dim)
229236

@@ -243,15 +250,17 @@ contains
243250
mask${select_subarray(rank, [(fi, 'i')])}$)
244251
end do
245252
else
253+
allocate(mean_, source = mean(x, ${fi}$, mask))
246254
do i = 1, size(x, ${fi}$)
247-
res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$)**order,&
255+
res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean_)**order,&
248256
#:if t1[0] == 'r'
249257
0._${k1}$,&
250258
#:else
251259
cmplx(0,0,kind=${k1}$),&
252260
#:endif
253261
mask${select_subarray(rank, [(fi, 'i')])}$)
254262
end do
263+
deallocate(mean_)
255264
end if
256265
#:endfor
257266
case default
@@ -277,6 +286,7 @@ contains
277286

278287
integer :: i
279288
real(dp) :: n${reduced_shape('x', rank, 'dim')}$
289+
real(dp), allocatable :: mean_${ranksuffix(rank-1)}$
280290

281291
n = count(mask, dim)
282292

@@ -291,11 +301,13 @@ contains
291301
0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
292302
end do
293303
else
304+
allocate(mean_, source = mean(x, ${fi}$, mask))
294305
do i = 1, size(x, ${fi}$)
295-
res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp))&
306+
res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean_)&
296307
**order,&
297308
0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
298309
end do
310+
deallocate(mean_)
299311
end if
300312
#:endfor
301313
case default

0 commit comments

Comments
 (0)