s_tst_signedrank_1s_core Subroutine

pure subroutine s_tst_signedrank_1s_core(x, mu0, w, p, h1)

! ==== Description
!! The 1-sample Wilcoxon signed rank test.

! ==== Declarations
  real(wp)         , intent(in)  :: x(:)     !! x vector (samples)
  real(wp)         , intent(in)  :: mu0      !! population mean (null hypothesis expected value)
  real(wp)         , intent(out) :: w        !! W statistic (sum of signed ranks)
  real(wp)         , intent(out) :: p        !! p-value
  character(len=*) , intent(in)  :: h1       !! \( H_{1} \): "two" (default), "lt", "gt"
  integer(i4)                    :: n        !! sample size
  integer(i4)                    :: m        !! non-zero sample size
  real(wp)                       :: mr       !! float sample size for x
  real(wp)         , allocatable :: d(:)     !! differences
  real(wp)         , allocatable :: dm(:)    !! absolute (modulus of) differences
  real(wp)         , allocatable :: ranks(:) !! ranks of absolute differences
  integer(i4)      , allocatable :: idx(:)   !! indeces for x > 0
  real(wp)                       :: rpos     !! sum of positive ranks
  real(wp)                       :: rneg     !! sum of negative ranks
  real(wp)                       :: mu       !! expected W under H0
  real(wp)                       :: s        !! standard deviation under H0
  real(wp)                       :: z        !! z statistic
  integer(i4)                    :: i, j

! ==== Instructions

  ! get dims and allocate
  n = size(x)
  allocate(d(n))

  ! compute differences and absolute values
  d  = x - mu0

  ! filter out zero differences
  m = count(abs(d) .gt. 0.0_wp)
  if (m .eq. 0) then
     w = 0.0_wp
     p = 1.0_wp
     return
  endif

  ! allocation based on new vector
  allocate(dm(m))
  allocate(ranks(m))
  allocate(idx(m))
  mr = real(m, kind=wp)

  ! extract non-zero entries and their signs
  i = 0
  do j = 1, n
     if (abs(d(j)) .gt. 0.0_wp) then
        i = i + 1
        dm(i)  = abs(d(j))
        idx(i) = j
     endif
  enddo

  ! rank non-zero differences
  call s_utl_rank(dm, ranks)

  ! compute rank sums
  rpos = 0.0_wp
  rneg = 0.0_wp
  do i = 1, m
     if (d(idx(i)) .gt. 0.0_wp) then
        rpos = rpos + ranks(i)
     else
        rneg = rneg + ranks(i)
     endif
  enddo

  ! W is the smaller of the rank sums
  w = min(rpos, rneg)

  ! expected value and standard deviation under H0
  mu = mr * (mr + 1.0_wp) / 4.0_wp
  s  = sqrt(mr * (mr + 1.0_wp) * (2.0_wp * mr + 1.0_wp) / 24.0_wp)

  ! z statistic
  z = (w - 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(d)
  deallocate(dm)
  deallocate(ranks)
  deallocate(idx)

end subroutine s_tst_signedrank_1s_core