stdlib_stats_var.fypp Source File


Source Code

#:include "common.fypp"
#:set RANKS = range(1, MAXRANK + 1)
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
submodule (stdlib_stats) stdlib_stats_var

  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
  use stdlib_error, only: error_stop
  use stdlib_optval, only: optval
  implicit none

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, 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
        ${t1}$ :: mean

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

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

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

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


  #: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, 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

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

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

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

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


  #:for k1, t1 in RC_KINDS_TYPES
    #:for rank in RANKS
      #:set RName = rname("var",rank, t1, k1)
      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
        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

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

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


  #: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, 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
        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

        res = 0._dp
        select case(dim)
          #:for fi in range(1, rank+1)
          case(${fi}$)
            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
            end do
          #:endfor
          case default
            call error_stop("ERROR (var): wrong dimension")
        end select
        res = res / (n - merge(1, 0, optval(corrected, .true.)))

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


  #: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, 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}$)
        mean = sum(x, mask) / n

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

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


  #: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, 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)
        mean = sum(real(x, dp), mask) / n

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

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


  #: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, 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
        real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
        ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$

        res = 0._${k1}$
        select case(dim)
          #:for fi in range(1, rank+1)
          case(${fi}$)
            n = count(mask, dim)
            mean = sum(x, dim, mask) / n
            do i = 1, size(x, dim)
              #:if t1[0] == 'r'
                res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2,&
              #:else
                res = res + merge( abs(x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2,&
              #:endif
                                  0._${k1}$,&
                                  mask${select_subarray(rank, [(fi, 'i')])}$)
            end do
          #:endfor
          case default
            call error_stop("ERROR (var): wrong dimension")
        end select
        res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0)))

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


  #: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, 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
        real(dp) :: n${reduced_shape('x', rank, 'dim')}$
        real(dp) :: mean${reduced_shape('x', rank, 'dim')}$

        res = 0._dp
        select case(dim)
          #:for fi in range(1, rank+1)
          case(${fi}$)
            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,&
                                  0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
            end do
          #:endfor
          case default
            call error_stop("ERROR (var): wrong dimension")
        end select
        res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0)))

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

end submodule