Skip to content

Addition of corrected in API of var (following #149) #151

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.

Sign up for GitHub

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 4 commits into from
Feb 24, 2020
Merged

Addition of corrected in API of var (following #149) #151

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
50367ff
Addition of corrected in API of var
jvdp1 Feb 22, 2020
fdfd7f0
var_corrected_dev:remove a unnecessary condition
jvdp1 Feb 22, 2020
21f68dd
var_corrected_dev: review following comments in #151
jvdp1 Feb 23, 2020
9af45ef
var_corrected_dev: addition of tests for dp, int, and complex(dp) for…
jvdp1 Feb 23, 2020
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
24 changes: 16 additions & 8 deletions src/stdlib_experimental_stats.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,9 +110,10 @@ module stdlib_experimental_stats
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(${k1}$) :: res
end function ${RName}$
#:endfor
Expand All @@ -121,9 +122,10 @@ module stdlib_experimental_stats
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(dp) :: res
end function ${RName}$
#:endfor
Expand All @@ -132,10 +134,11 @@ module stdlib_experimental_stats
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
Expand All @@ -144,10 +147,11 @@ module stdlib_experimental_stats
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
Expand All @@ -156,9 +160,10 @@ module stdlib_experimental_stats
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(${k1}$) :: res
end function ${RName}$
#:endfor
Expand All @@ -167,9 +172,10 @@ module stdlib_experimental_stats
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(dp) :: res
end function ${RName}$
#:endfor
Expand All @@ -178,10 +184,11 @@ module stdlib_experimental_stats
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
Expand All @@ -190,10 +197,11 @@ module stdlib_experimental_stats
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
end function ${RName}$
#:endfor
Expand Down
28 changes: 20 additions & 8 deletions src/stdlib_experimental_stats.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,17 +55,22 @@ end program demo_mean

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

The variance is defined as the best unbiased estimator and is computed as:
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
```

where n is the number of elements.

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


### Syntax

`result = var(array [, mask])`
`result = var(array [, mask [, corrected]])`

`result = var(array, dim [, mask])`
`result = var(array, dim [, mask [, corrected]])`

### Arguments

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

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

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

### 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 `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
If `array` is of type `integer`, the result is of type `real(dp)`.

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
ray` with dimension `dim` dropped is returned.
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.

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

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

### Example

```fortran
Expand All @@ -93,10 +101,14 @@ program demo_var
implicit none
real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
print *, var(x) !returns 3.5
print *, var(x, corrected = .false.) !returns 2.9167
print *, var( reshape(x, [ 2, 3 ] )) !returns 3.5
print *, var( reshape(x, [ 2, 3 ] ), 1) !returns [0.5, 0.5, 0.5]
print *, var( reshape(x, [ 2, 3 ] ), 1,&
reshape(x, [ 2, 3 ] ) > 3.) !returns [0., NaN, 0.5]
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]
end program demo_var
```

62 changes: 36 additions & 26 deletions src/stdlib_experimental_stats_var.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,10 @@ contains
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(${k1}$) :: res

real(${k1}$) :: n
Expand All @@ -26,13 +27,13 @@ contains
return
end if

n = real(size(x, kind = int64), ${k1}$)
n = size(x, kind = int64)
mean = sum(x) / n

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

end function ${RName}$
Expand All @@ -43,9 +44,10 @@ contains
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(dp) :: res

real(dp) :: n, mean
Expand All @@ -55,10 +57,10 @@ contains
return
end if

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

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

end function ${RName}$
#:endfor
Expand All @@ -68,10 +70,11 @@ contains
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$

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

end function ${RName}$
#:endfor
Expand All @@ -110,10 +113,11 @@ contains
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
logical, intent(in), optional :: corrected
real(dp) :: res${reduced_shape('x', rank, 'dim')}$

integer :: i
Expand All @@ -129,7 +133,7 @@ contains
select case(dim)
#:for fi in range(1, rank+1)
case(${fi}$)
n = real(size(x, dim), dp)
n = size(x, dim)
mean = sum(real(x, dp), dim) / n
do i = 1, size(x, dim)
res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2
Expand All @@ -138,7 +142,7 @@ contains
case default
call error_stop("ERROR (mean): wrong dimension")
end select
res = res / (n - 1._dp)
res = res / (n - merge(1, 0, optval(corrected, .true.)))

end function ${RName}$
#:endfor
Expand All @@ -148,22 +152,24 @@ contains
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(${k1}$) :: res

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

n = real(count(mask, kind = int64), ${k1}$)
n = count(mask, kind = int64)
mean = sum(x, mask) / n

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

end function ${RName}$
#:endfor
Expand All @@ -173,17 +179,19 @@ contains
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
module function ${RName}$(x, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(dp) :: res

real(dp) :: n, mean

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

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

end function ${RName}$
#:endfor
Expand All @@ -193,10 +201,11 @@ contains
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$

integer :: i
Expand All @@ -207,7 +216,7 @@ contains
select case(dim)
#:for fi in range(1, rank+1)
case(${fi}$)
n = real(count(mask, dim), ${k1}$)
n = count(mask, dim)
mean = sum(x, dim, mask) / n
do i = 1, size(x, dim)
#:if t1[0] == 'r'
Expand All @@ -222,7 +231,7 @@ contains
case default
call error_stop("ERROR (mean): wrong dimension")
end select
res = res / (n - 1._${k1}$)
res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0)))

end function ${RName}$
#:endfor
Expand All @@ -232,10 +241,11 @@ contains
#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
#:set RName = rname("var_mask",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
module function ${RName}$(x, dim, mask, corrected) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
logical, intent(in), optional :: corrected
real(dp) :: res${reduced_shape('x', rank, 'dim')}$

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

end function ${RName}$
#:endfor
Expand Down
Loading