Skip to content

Commit 340309e

Browse files
jvdp1jvdp1
committed
covariance: init
1 parent 25965bc commit 340309e

File tree

5 files changed

+233
-1
lines changed

5 files changed

+233
-1
lines changed

‎src/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ set(fppFiles
55
stdlib_experimental_io.fypp
66
stdlib_experimental_optval.fypp
77
stdlib_experimental_stats.fypp
8+
stdlib_experimental_stats_cov.fypp
89
stdlib_experimental_stats_mean.fypp
910
stdlib_experimental_stats_moment.fypp
1011
stdlib_experimental_stats_var.fypp

‎src/stdlib_experimental_stats.fypp‎

Lines changed: 27 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ module stdlib_experimental_stats
88
implicit none
99
private
1010
! Public API
11-
public :: mean, moment, var
11+
public :: cov, mean, moment, var
1212

1313
interface mean
1414
#:for k1, t1 in RC_KINDS_TYPES
@@ -376,4 +376,30 @@ module stdlib_experimental_stats
376376
#:endfor
377377
end interface moment
378378

379+
interface cov
380+
#:for k1, t1 in RC_KINDS_TYPES
381+
#:set RName = rname("cov",2, t1, k1)
382+
module function ${RName}$(x, dim, mask, corrected) result(res)
383+
${t1}$, intent(in) :: x(:, :)
384+
integer, intent(in) :: dim
385+
logical, intent(in), optional :: mask
386+
logical, intent(in), optional :: corrected
387+
real(${k1}$) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
388+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
389+
end function ${RName}$
390+
#:endfor
391+
392+
#:for k1, t1 in INT_KINDS_TYPES
393+
#:set RName = rname("cov",2, t1, k1, 'dp')
394+
module function ${RName}$(x, dim, mask, corrected) result(res)
395+
${t1}$, intent(in) :: x(:, :)
396+
integer, intent(in) :: dim
397+
logical, intent(in), optional :: mask
398+
logical, intent(in), optional :: corrected
399+
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
400+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
401+
end function ${RName}$
402+
#:endfor
403+
end interface cov
404+
379405
end module stdlib_experimental_stats

‎src/stdlib_experimental_stats_cov.fypp‎

Lines changed: 99 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,99 @@
1+
#:include "common.fypp"
2+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
submodule (stdlib_experimental_stats) stdlib_experimental_stats_cov
4+
5+
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
6+
use stdlib_experimental_error, only: error_stop
7+
use stdlib_experimental_optval, only: optval
8+
implicit none
9+
10+
contains
11+
12+
#:for k1, t1 in RC_KINDS_TYPES
13+
#:set RName = rname("cov",2, t1, k1)
14+
module function ${RName}$(x, dim, mask, corrected) result(res)
15+
${t1}$, intent(in) :: x(:, :)
16+
integer, intent(in) :: dim
17+
logical, intent(in), optional :: mask
18+
logical, intent(in), optional :: corrected
19+
real(${k1}$) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
20+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
21+
22+
integer :: i
23+
real(${k1}$) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
24+
real(${k1}$) :: center(size(x, 1),size(x, 2))
25+
26+
if (.not.optval(mask, .true.)) then
27+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
28+
return
29+
end if
30+
31+
mean_ = mean(x, dim)
32+
select case(dim)
33+
case(1)
34+
do i = 1, size(x, 1)
35+
center(i, :) = x(i, :) - mean_
36+
end do
37+
#:if t1[0] == 'r'
38+
res = matmul( transpose(center), center)
39+
#:else
40+
error stop
41+
#:endif
42+
case(2)
43+
do i = 1, size(x, 2)
44+
center(:, i) = x(:, i) - mean_
45+
end do
46+
#:if t1[0] == 'r'
47+
res = matmul( center, transpose(center))
48+
#:else
49+
error stop
50+
#:endif
51+
case default
52+
call error_stop("ERROR (mean): wrong dimension")
53+
end select
54+
res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
55+
56+
end function ${RName}$
57+
#:endfor
58+
59+
60+
#:for k1, t1 in INT_KINDS_TYPES
61+
#:set RName = rname("cov",2, t1, k1, 'dp')
62+
module function ${RName}$(x, dim, mask, corrected) result(res)
63+
${t1}$, intent(in) :: x(:, :)
64+
integer, intent(in) :: dim
65+
logical, intent(in), optional :: mask
66+
logical, intent(in), optional :: corrected
67+
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
68+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
69+
70+
integer :: i
71+
real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
72+
real(dp) :: center(size(x, 1),size(x, 2))
73+
74+
if (.not.optval(mask, .true.)) then
75+
res = ieee_value(1._dp, ieee_quiet_nan)
76+
return
77+
end if
78+
79+
mean_ = mean(x, dim)
80+
select case(dim)
81+
case(1)
82+
do i = 1, size(x, 1)
83+
center(i, :) = real(x(i, :), dp) - mean_
84+
end do
85+
res = matmul( transpose(center), center)
86+
case(2)
87+
do i = 1, size(x, 2)
88+
center(:, i) = real(x(:, i), dp) - mean_
89+
end do
90+
res = matmul( center, transpose(center))
91+
case default
92+
call error_stop("ERROR (mean): wrong dimension")
93+
end select
94+
res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
95+
96+
end function ${RName}$
97+
#:endfor
98+
99+
end submodule

‎src/tests/stats/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
ADDTEST(cov)
12
ADDTEST(mean)
23
ADDTEST(moment)
34
ADDTEST(rawmoment)

‎src/tests/stats/test_cov.f90‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
program test_moment
2+
use stdlib_experimental_error, only: check
3+
use stdlib_experimental_kinds, only: sp, dp, int32, int64
4+
use stdlib_experimental_stats, only: cov, var
5+
use,intrinsic :: ieee_arithmetic, only : ieee_is_nan
6+
implicit none
7+
8+
9+
real(sp), parameter :: sptol = 1000 * epsilon(1._sp)
10+
real(dp), parameter :: dptol = 1000 * epsilon(1._dp)
11+
12+
real(dp) :: d(4, 3) = reshape([1._dp, 3._dp, 5._dp, 7._dp,&
13+
2._dp, 4._dp, 6._dp, 8._dp,&
14+
9._dp, 10._dp, 11._dp, 12._dp], [4, 3])
15+
16+
17+
call test_dp(d)
18+
call test_int32(int(d, int32))
19+
20+
contains
21+
22+
subroutine test_dp(x2)
23+
real(dp), intent(in) :: x2(:, :)
24+
25+
call check( any(ieee_is_nan(cov(x2, 1, mask = .false.)))&
26+
, 'dp check 1')
27+
call check( any(ieee_is_nan(cov(x2, 2, mask = .false.)))&
28+
, 'dp check 2')
29+
30+
call check( all( abs( cov(x2, 1) - reshape([&
31+
60._dp/9, 60._dp/9, 30._dp/9&
32+
,60._dp/9, 60._dp/9, 30._dp/9&
33+
,30._dp/9, 30._dp/9, 15._dp/9]&
34+
,[ size(x2, 2), size(x2, 2)])&
35+
) < dptol)&
36+
, 'dp check 3')
37+
call check( all( abs( cov(x2, 2) - reshape([&
38+
19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
39+
,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
40+
, 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
41+
,[ size(x2, 1), size(x2, 1)])&
42+
) < dptol)&
43+
, 'dp check 4')
44+
45+
call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
46+
60._dp/9, 60._dp/9, 30._dp/9&
47+
,60._dp/9, 60._dp/9, 30._dp/9&
48+
,30._dp/9, 30._dp/9, 15._dp/9]&
49+
*(size(x2, 1)-1._dp)/size(x2, 1)&
50+
,[ size(x2, 2), size(x2, 2)])&
51+
) < dptol)&
52+
, 'dp check 5')
53+
call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
54+
19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
55+
,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
56+
, 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
57+
*(size(x2, 2)-1._dp)/size(x2, 2)&
58+
,[ size(x2, 1), size(x2, 1)])&
59+
) < dptol)&
60+
, 'dp check 6')
61+
62+
end subroutine test_dp
63+
64+
subroutine test_int32(x2)
65+
integer(int32), intent(in) :: x2(:, :)
66+
67+
call check( any(ieee_is_nan(cov(x2, 1, mask = .false.)))&
68+
, 'int32 check 1')
69+
call check( any(ieee_is_nan(cov(x2, 2, mask = .false.)))&
70+
, 'int32 check 2')
71+
72+
call check( all( abs( cov(x2, 1) - reshape([&
73+
60._dp/9, 60._dp/9, 30._dp/9&
74+
,60._dp/9, 60._dp/9, 30._dp/9&
75+
,30._dp/9, 30._dp/9, 15._dp/9]&
76+
,[ size(x2, 2), size(x2, 2)])&
77+
) < dptol)&
78+
, 'int32 check 3')
79+
call check( all( abs( cov(x2, 2) - reshape([&
80+
19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
81+
,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
82+
, 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
83+
,[ size(x2, 1), size(x2, 1)])&
84+
) < dptol)&
85+
, 'int32 check 4')
86+
87+
call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
88+
60._dp/9, 60._dp/9, 30._dp/9&
89+
,60._dp/9, 60._dp/9, 30._dp/9&
90+
,30._dp/9, 30._dp/9, 15._dp/9]&
91+
*(size(x2, 1)-1._dp)/size(x2, 1)&
92+
,[ size(x2, 2), size(x2, 2)])&
93+
) < dptol)&
94+
, 'int32 check 5')
95+
call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
96+
19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
97+
,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
98+
, 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
99+
*(size(x2, 2)-1._dp)/size(x2, 2)&
100+
,[ size(x2, 1), size(x2, 1)])&
101+
) < dptol)&
102+
, 'int32 check 6')
103+
104+
end subroutine test_int32
105+
end program test_moment

0 commit comments

Comments
 (0)