Merge pull request fortran-lang#192 from jvdp1/cov_correction

Addition of tests and clarification for cov

Author: milancurcic

src/stdlib_experimental_stats_cov.fypp

@@ -178,52 +178,58 @@ contains
                           , merge(size(x, 1), size(x, 2), mask = 1<dim))
 
       integer :: i, j, n
-      ${t1}$ :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
-      ${t1}$ :: center(size(x, 1),size(x, 2))
+      ${t1}$ :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      ${t1}$ :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))
 
-      mean_ = mean(x, dim, mask = mask)
       select case(dim)
         case(1)
-          do i = 1, size(x, 1)
-            center(i, :) = merge( x(i, :) - mean_,&
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
+             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
               #:if t1[0] == 'r'
                 0._${k1}$,&
               #:else
                 cmplx(0,0,kind=${k1}$),&
               #:endif
-                mask(i, :))
-          end do
-          #:if t1[0] == 'r'
-            res = matmul( transpose(center), center)
-          #:else
-            res = matmul( transpose(conjg(center)), center)
-          #:endif
-          do j = 1, size(res, 2)
-            do i = 1, size(res, 1)
-              n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
-              res(i, j) = res(i, j) / (n - merge(1, 0,&
+                mask_)
+             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask_)
+
+              n = count(mask_)
+              res(j, i) = dot_product( centerj_, centeri_)&
+                           / (n - merge(1, 0,&
                             optval(corrected, .true.) .and. n > 0))
             end do
           end do
         case(2)
-          do i = 1, size(x, 2)
-            center(:, i) = merge( x(:, i) - mean_,&
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
+             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
               #:if t1[0] == 'r'
                 0._${k1}$,&
               #:else
                 cmplx(0,0,kind=${k1}$),&
               #:endif
-                mask(:, i))
-          end do
-          #:if t1[0] == 'r'
-            res = matmul( center, transpose(center))
-          #:else
-            res = matmul( center, transpose(conjg(center)))
-          #:endif
-          do j = 1, size(res, 2)
-            do i = 1, size(res, 1)
-              n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
-              res(i, j) = res(i, j) / (n - merge(1, 0,&
+                mask_)
+             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask_)
+
+              n = count(mask_)
+              res(j, i) = dot_product( centeri_, centerj_)&
+                           / (n - merge(1, 0,&
                             optval(corrected, .true.) .and. n > 0))
             end do
           end do
@@ -246,36 +252,38 @@ contains
                           , merge(size(x, 1), size(x, 2), mask = 1<dim))
 
       integer :: i, j, n
-      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
-      real(dp) :: center(size(x, 1),size(x, 2))
+      real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))
 
-      mean_ = mean(x, dim, mask = mask)
       select case(dim)
         case(1)
-          do i = 1, size(x, 1)
-            center(i, :) = merge( x(i, :) - mean_,&
-                0._dp,&
-                mask(i, :))
-          end do
-          res = matmul( transpose(center), center)
-          do j = 1, size(res, 2)
-            do i = 1, size(res, 1)
-              n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
-              res(i, j) = res(i, j) / (n - merge(1, 0,&
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+              mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
+              centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
+                 0._dp, mask_)
+              centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
+                 0._dp, mask_)
+
+              n = count(mask_)
+              res(j, i) = dot_product( centerj_, centeri_)&
+                           / (n - merge(1, 0,&
                             optval(corrected, .true.) .and. n > 0))
             end do
           end do
         case(2)
-          do i = 1, size(x, 2)
-            center(:, i) = merge( x(:, i) - mean_,&
-                0._dp,&
-                mask(:, i))
-          end do
-            res = matmul( center, transpose(center))
-          do j = 1, size(res, 2)
-            do i = 1, size(res, 1)
-              n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
-              res(i, j) = res(i, j) / (n - merge(1, 0,&
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+              mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
+              centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
+                 0._dp, mask_)
+              centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
+                 0._dp, mask_)
+
+              n = count(mask_)
+              res(j, i) = dot_product( centeri_, centerj_)&
+                           / (n - merge(1, 0,&
                             optval(corrected, .true.) .and. n > 0))
             end do
           end do

src/tests/stats/test_cov.f90

@@ -126,6 +126,14 @@ subroutine test_sp(x, x2)
             ) < sptol)&
             , 'sp check 16')
 
+        call check( all( abs( cov(x2, 1, mask = x2 < 1000) - cov(x2, 1))&
+            < sptol)&
+            , 'sp check 17')
+
+        call check( all( abs( cov(x2, 2, mask = x2 < 1000) - cov(x2, 2))&
+            < sptol)&
+            , 'sp check 18')
+
     end subroutine test_sp
 
     subroutine test_dp(x, x2)
@@ -213,6 +221,13 @@ subroutine test_dp(x, x2)
             ) < dptol)&
             , 'dp check 16')
 
+        call check( all( abs( cov(x2, 1, mask = x2 < 1000) - cov(x2, 1))&
+            < dptol)&
+            , 'dp check 17')
+
+        call check( all( abs( cov(x2, 2, mask = x2 < 1000) - cov(x2, 2))&
+            < dptol)&
+            , 'dp check 18')
 
     end subroutine test_dp
 
@@ -301,6 +316,14 @@ subroutine test_int32(x, x2)
             ) < dptol)&
             , 'int32 check 16')
 
+        call check( all( abs( cov(x2, 1, mask = x2 < 1000) - cov(x2, 1))&
+            < dptol)&
+            , 'int32 check 17')
+
+        call check( all( abs( cov(x2, 2, mask = x2 < 1000) - cov(x2, 2))&
+            < dptol)&
+            , 'int32 check 18')
+
     end subroutine test_int32
 
     subroutine test_int64(x, x2)
@@ -388,6 +411,14 @@ subroutine test_int64(x, x2)
             ) < dptol)&
             , 'int64 check 16')
 
+        call check( all( abs( cov(x2, 1, mask = x2 < 1000) - cov(x2, 1))&
+            < dptol)&
+            , 'int64 check 17')
+
+        call check( all( abs( cov(x2, 2, mask = x2 < 1000) - cov(x2, 2))&
+            < dptol)&
+            , 'int64 check 18')
+
     end subroutine test_int64
 
     subroutine test_csp(x, x2)
@@ -459,20 +490,29 @@ subroutine test_csp(x, x2)
 !        call check( ieee_is_nan(real(cd(3,3)))&
 !            , 'csp check 10 bis')
 
+
+        call check( all( abs( cov(x2, 1, mask = aimag(x2) < 8) - cov(x2, 1))&
+            < sptol)&
+            , 'csp check 11')
+
+        call check( all( abs( cov(x2, 2, mask = aimag(x2) < 8) - cov(x2, 2))&
+            < sptol)&
+            , 'csp check 12')
+
         call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6) - reshape([&
              (4._sp,0._sp), (0._sp,2._sp)&
              ,(0._sp,-2._sp), (2._sp,0._sp)]&
             ,[ size(x2, 1), size(x2, 1)])&
             ) < sptol)&
-            , 'csp check 11')
+            , 'csp check 13')
 
         call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6, corrected = .false.) -&
             reshape([&
              (2.6666666666666666_sp,0._sp), (0._sp,1._sp)&
              ,(0._sp,-1._sp), (1._sp,0._sp)]&
             ,[ size(x2, 1), size(x2, 1)])&
             ) < sp)&
-            , 'csp check 12')
+            , 'csp check 14')
 
     end subroutine test_csp
 
@@ -546,20 +586,28 @@ subroutine test_cdp(x, x2)
 !        call check( ieee_is_nan(real(cd(3,3)))&
 !            , 'cdp check 10 bis')
 
+        call check( all( abs( cov(x2, 1, mask = aimag(x2) < 8) - cov(x2, 1))&
+            < dptol)&
+            , 'cdp check 11')
+
+        call check( all( abs( cov(x2, 2, mask = aimag(x2) < 8) - cov(x2, 2))&
+            < dptol)&
+            , 'cdp check 12')
+
         call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6) - reshape([&
              (4._dp,0._dp), (0._dp,2._dp)&
              ,(0._dp,-2._dp), (2._dp,0._dp)]&
             ,[ size(x2, 1), size(x2, 1)])&
             ) < dptol)&
-            , 'cdp check 11')
+            , 'cdp check 13')
 
         call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6, corrected = .false.) -&
             reshape([&
              (2.6666666666666666_dp,0._dp), (0._dp,1._dp)&
              ,(0._dp,-1._dp), (1._dp,0._dp)]&
             ,[ size(x2, 1), size(x2, 1)])&
             ) < dptol)&
-            , 'cdp check 12')
+            , 'cdp check 14')
 
     end subroutine test_cdp