Skip to content

Addition of a function to compute the k-th order central moment of an array #153

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Mar 14, 2020
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -6,6 +6,7 @@ set(fppFiles
stdlib_experimental_optval.fypp
stdlib_experimental_stats.fypp
stdlib_experimental_stats_mean.fypp
stdlib_experimental_stats_moment.fypp
stdlib_experimental_stats_var.fypp
stdlib_experimental_quadrature.fypp
stdlib_experimental_quadrature_trapz.fypp
6 changes: 6 additions & 0 deletions src/Makefile.manual
Original file line number Diff line number Diff line change
@@ -8,6 +8,7 @@ SRC = f18estop.f90 \
stdlib_experimental_quadrature_trapz.f90 \
stdlib_experimental_stats.f90 \
stdlib_experimental_stats_mean.f90 \
stdlib_experimental_stats_moment.f90 \
stdlib_experimental_stats_var.f90

LIB = libstdlib.a
@@ -46,6 +47,10 @@ stdlib_experimental_stats_mean.o: \
stdlib_experimental_optval.o \
stdlib_experimental_kinds.o \
stdlib_experimental_stats.o
stdlib_experimental_stats_moment.o: \
stdlib_experimental_optval.o \
stdlib_experimental_kinds.o \
stdlib_experimental_stats.o
stdlib_experimental_stats_var.o: \
stdlib_experimental_optval.o \
stdlib_experimental_kinds.o \
@@ -56,4 +61,5 @@ stdlib_experimental_io.f90: stdlib_experimental_io.fypp
stdlib_experimental_quadrature.f90: stdlib_experimental_quadrature.fypp
stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp
stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp
stdlib_experimental_stats_moment.f90: stdlib_experimental_stats_moment.fypp
stdlib_experimental_stats_var.f90: stdlib_experimental_stats_var.fypp
103 changes: 102 additions & 1 deletion src/stdlib_experimental_stats.fypp
Original file line number Diff line number Diff line change
@@ -7,7 +7,7 @@ module stdlib_experimental_stats
implicit none
private
! Public API
public :: mean, var
public :: mean, moment, var

interface mean
#:for k1, t1 in RC_KINDS_TYPES
@@ -209,5 +209,106 @@ module stdlib_experimental_stats

end interface var

interface moment
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_all",rank, t1, k1)
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in), optional :: mask
${t1}$ :: res
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_all",rank, t1, k1, 'dp')
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in), optional :: mask
real(dp) :: res
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment",rank, t1, k1)
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in), optional :: mask
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment",rank, t1, k1, 'dp')
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in), optional :: mask
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask_all",rank, t1, k1)
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in) :: mask${ranksuffix(rank)}$
${t1}$ :: res
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask_all",rank, t1, k1, 'dp')
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in) :: mask${ranksuffix(rank)}$
real(dp) :: res
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask",rank, t1, k1)
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask",rank, t1, k1, 'dp')
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
#:endfor
end interface moment

end module stdlib_experimental_stats
86 changes: 80 additions & 6 deletions src/stdlib_experimental_stats.md
Original file line number Diff line number Diff line change
@@ -1,9 +1,29 @@
# Descriptive statistics

## Implemented

* `mean`
* `var`
## Implemented
<!-- vim-markdown-toc GFM -->

* [`mean` - mean of array elements](#mean---mean-of-array-elements)
* [Description](#description)
* [Syntax](#syntax)
* [Arguments](#arguments)
* [Return value](#return-value)
* [Example](#example)
* [`moment` - central moment of array elements](#moment---central-moment-of-array-elements)
* [Description](#description-1)
* [Syntax](#syntax-1)
* [Arguments](#arguments-1)
* [Return value](#return-value-1)
* [Example](#example-1)
* [`var` - variance of array elements](#var---variance-of-array-elements)
* [Description](#description-2)
* [Syntax](#syntax-2)
* [Arguments](#arguments-2)
* [Return value](#return-value-2)
* [Example](#example-2)

<!-- vim-markdown-toc -->

## `mean` - mean of array elements

@@ -28,7 +48,7 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a
### 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 `double precision`.
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.

@@ -49,6 +69,60 @@ program demo_mean
end program demo_mean
```

## `moment` - central moment of array elements

### Description

Returns the _k_-th order central moment 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`.

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.

### Syntax

`result = moment(array, order [, mask])`

`result = moment(array, order, dim [, mask])`

### Arguments

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

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

`mask` (optional): Shall be of type `logical` and either by 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 `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`.

### Example

```fortran
program demo_moment
use stdlib_experimental_stats, only: moment
implicit none
real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
print *, moment(x, 2) !returns 2.9167
print *, moment( reshape(x, [ 2, 3 ] ), 2) !returns 2.9167
print *, moment( reshape(x, [ 2, 3 ] ), 2, 1) !returns [0.25, 0.25, 0.25]
print *, moment( reshape(x, [ 2, 3 ] ), 2, 1,&
reshape(x, [ 2, 3 ] ) > 3.) !returns [NaN, 0., 0.25]
end program demo_moment
```

## `var` - variance of array elements

### Description
@@ -58,7 +132,7 @@ Returns the variance of all the elements of `array`, or of the elements of `arra
Per default, the variance is defined as the best unbiased estimator and is computed as:

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

where n is the number of elements.
@@ -108,7 +182,7 @@ program demo_var
reshape(x, [ 2, 3 ] ) > 3.) !returns [NaN, NaN, 0.5]
print *, var( reshape(x, [ 2, 3 ] ), 1,&
reshape(x, [ 2, 3 ] ) > 3.,&
corrected=.false.) !returns [NaN, 0., 0.5]
corrected=.false.) !returns [NaN, 0., 0.25]
end program demo_var
```

265 changes: 265 additions & 0 deletions src/stdlib_experimental_stats_moment.fypp
Original file line number Diff line number Diff line change
@@ -0,0 +1,265 @@
#:include "common.fypp"
#:set RANKS = range(1, MAXRANK + 1)
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
submodule (stdlib_experimental_stats) stdlib_experimental_stats_moment

use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
use stdlib_experimental_error, only: error_stop
use stdlib_experimental_optval, only: optval
implicit none

contains

#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_all",rank, t1, k1)
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in), optional :: mask
${t1}$ :: res

real(${k1}$) :: n
${t1}$ :: mean

if (.not.optval(mask, .true.)) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if

n = size(x, kind = int64)
mean = sum(x) / n

res = sum((x - mean)**order) / n
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not

res = mean((x - mean(x))**order)

here? Is there a concern about the overhead from calling mean()?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was indeed to avoid the overhead of calling mean (especially when mean is an array).


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


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_all",rank, t1, k1, 'dp')
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in), optional :: mask
real(dp) :: res

real(dp) :: n, mean

if (.not.optval(mask, .true.)) then
res = ieee_value(1._dp, ieee_quiet_nan)
return
end if

n = size(x, kind = int64)
mean = sum(real(x, dp)) / n

res = sum((real(x, dp) - mean)**order) / n

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


#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment",rank, t1, k1)
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in), optional :: mask
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$

integer :: i
real(${k1}$) :: n
${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$

if (.not.optval(mask, .true.)) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if

n = size(x, dim)
mean = sum(x, dim) / n

res = 0
select case(dim)
#:for fi in range(1, rank+1)
case(${fi}$)
do i = 1, size(x, dim)
res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**order
end do
#:endfor
case default
call error_stop("ERROR (moment): wrong dimension")
end select
res = res / n

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


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment",rank, t1, k1, 'dp')
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in), optional :: mask
real(dp) :: res${reduced_shape('x', rank, 'dim')}$

integer :: i
real(dp) :: n
real(dp) :: mean${reduced_shape('x', rank, 'dim')}$

if (.not.optval(mask, .true.)) then
res = ieee_value(1._dp, ieee_quiet_nan)
return
end if

n = size(x, dim)
mean = sum(real(x, dp), dim) / n

res = 0
select case(dim)
#:for fi in range(1, rank+1)
case(${fi}$)
do i = 1, size(x, dim)
res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**order
end do
#:endfor
case default
call error_stop("ERROR (moment): wrong dimension")
end select
res = res / n

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


#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask_all",rank, t1, k1)
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in) :: mask${ranksuffix(rank)}$
${t1}$ :: res

real(${k1}$) :: n
${t1}$ :: mean

n = count(mask, kind = int64)
mean = sum(x, mask) / n

res = sum((x - mean)**order, mask) / n

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


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask_all",rank, t1, k1, 'dp')
module function ${RName}$(x, order, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
logical, intent(in) :: mask${ranksuffix(rank)}$
real(dp) :: res

real(dp) :: n, mean

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

res = sum((real(x, dp) - mean)**order, mask) / n

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


#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask",rank, t1, k1)
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$

integer :: i
real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$

n = count(mask, dim)
mean = sum(x, dim, mask) / n

res = 0
select case(dim)
#:for fi in range(1, rank+1)
case(${fi}$)
do i = 1, size(x, dim)
res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**order,&
#:if t1[0] == 'r'
0._${k1}$,&
#:else
cmplx(0,0,kind=${k1}$),&
#:endif
mask${select_subarray(rank, [(fi, 'i')])}$)
end do
#:endfor
case default
call error_stop("ERROR (moment): wrong dimension")
end select
res = res / n

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


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("moment_mask",rank, t1, k1, 'dp')
module function ${RName}$(x, order, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: order
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
real(dp) :: res${reduced_shape('x', rank, 'dim')}$

integer :: i
real(dp) :: n${reduced_shape('x', rank, 'dim')}$
real(dp) :: mean${reduced_shape('x', rank, 'dim')}$

n = count(mask, dim)
mean = sum(real(x, dp), dim, mask) / n

res = 0
select case(dim)
#:for fi in range(1, rank+1)
case(${fi}$)
do i = 1, size(x, dim)
res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**order,&
0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
end do
#:endfor
case default
call error_stop("ERROR (moment): wrong dimension")
end select
res = res / n

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

end submodule
1 change: 1 addition & 0 deletions src/tests/stats/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
ADDTEST(mean)
ADDTEST(moment)
ADDTEST(var)
ADDTEST(varn)

2 changes: 1 addition & 1 deletion src/tests/stats/Makefile.manual
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
PROGS_SRC = test_mean.f90 test_var.f90
PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90

include ../Makefile.manual.test.mk
706 changes: 706 additions & 0 deletions src/tests/stats/test_moment.f90

Large diffs are not rendered by default.