s_tst_ranksum_core Subroutine

private pure subroutine s_tst_ranksum_core(x1, x2, u, p, h1)

The ranks sum test (Wilcoxon rank-sum test or Mann–Whitney U test).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: x1(:)

x1 vector (samples)
real(kind=wp), intent(in) :: x2(:)

x2 vector (samples)
real(kind=wp), intent(out) :: u

U statistic
real(kind=wp), intent(out) :: p

p-value
character(len=*), intent(in) :: h1

option: "two" (default), "lt", or "gt"

Calls

proc~~s_tst_ranksum_core~~CallsGraph proc~s_tst_ranksum_core s_tst_ranksum_core proc~f_dst_norm_cdf_core f_dst_norm_cdf_core proc~s_tst_ranksum_core->proc~f_dst_norm_cdf_core proc~s_utl_rank s_utl_rank proc~s_tst_ranksum_core->proc~s_utl_rank

Called by

proc~~s_tst_ranksum_core~~CalledByGraph proc~s_tst_ranksum_core s_tst_ranksum_core proc~s_tst_ranksum s_tst_ranksum proc~s_tst_ranksum->proc~s_tst_ranksum_core interface~fsml_ranksum fsml_ranksum interface~fsml_ranksum->proc~s_tst_ranksum

Source Code

pure subroutine s_tst_ranksum_core(x1, x2, u, p, h1)

! ==== Description
!! The ranks sum test (Wilcoxon rank-sum test or Mann–Whitney U test).

! ==== Declarations
  real(wp)         , intent(in)  :: x1(:)    !! x1 vector (samples)
  real(wp)         , intent(in)  :: x2(:)    !! x2 vector (samples)
  real(wp)         , intent(out) :: u        !! U statistic
  real(wp)         , intent(out) :: p        !! p-value
  character(len=*) , intent(in)  :: h1       !! \( H_{1} \) option: "two" (default), "lt", or "gt"
  integer(i4)                    :: n1       !! sample size for x1
  integer(i4)                    :: n2       !! sample size for x2
  real(wp)                       :: n1r      !! float sample size for x1
  real(wp)                       :: n2r      !! float sample size for x2
  real(wp)         , allocatable :: x(:)     !! combined x1 and x2
  real(wp)         , allocatable :: ranks(:) !! stores ranks
  real(wp)                       :: rx1      !! rank sum for x1
  real(wp)                       :: rx2      !! rank sum for x2
  real(wp)                       :: u1       !! U statistic 1
  real(wp)                       :: u2       !! U statistic 2
  real(wp)                       :: mu       !! expected value of U (under H0)
  real(wp)                       :: s        !! standard deviation of U (under H0)
  real(wp)                       :: z        !! z statistic
  integer(i4)                    :: i

! ==== Instructions

  ! get dims and combine and rank all values
  n1 = size(x1)
  n2 = size(x2)
  n1r = real(n1, kind=wp)
  n2r = real(n2, kind=wp)
  allocate(x(n1+n2))
  allocate(ranks(n1+n2))
  x(1:n1)  = x1
  x(n1+1:) = x2

  ! assigns ranks, use stdlib procedure
  call s_utl_rank(x, ranks)

  ! sum ranks of sample x
  rx1 = sum( ranks(1:n1) )
  rx2 = sum( ranks(n1+1:) )

  ! compute U statistics
  u1 = n1r * n2r + (n1r * (n1r + 1.0_wp) / 2.0_wp) - rx1
  u2 = n1r * n2r + (n2r * (n2r + 1.0_wp) / 2.0_wp) - rx2
  u  = min(u1, u2)

  ! expected value and standard deviation of U under H0
  mu = n1r * n2r / 2.0_wp
  s  = sqrt(n1r * n2r * (n1r + n2r + 1.0_wp) / 12.0_wp)

  ! z statistic (U is approximately normal for large samples)
  z = (u - mu) / s

  ! get p-value
  select case(h1)
     ! less than
     case("lt")
        p = f_dst_norm_cdf_core(z, mu=0.0_wp, sigma=1.0_wp, tail="left")
     ! greater than
     case("gt")
        p = f_dst_norm_cdf_core(z, mu=0.0_wp, sigma=1.0_wp, tail="right")
     ! two-sided
     case("two")
        p = f_dst_norm_cdf_core(z, mu=0.0_wp, sigma=1.0_wp, tail="two")
  end select

  ! deallocate
  deallocate(x)
  deallocate(ranks)

end subroutine s_tst_ranksum_core