s_tst_kruskalwallis_core Subroutine

pure subroutine s_tst_kruskalwallis_core(x, h, df, p)

! ==== Description
!! Kruskal-Wallis H-test for independent samples. No tie correction.

! ==== Declarations
  real(wp)   , intent(in)  :: x(:,:)     !! 2D array, each column is a group
  real(wp)   , intent(out) :: h          !! Kruskal-Wallis H-statistic
  real(wp)   , intent(out) :: df         !! degrees of freedom (k - 1)
  real(wp)   , intent(out) :: p          !! p-value from chi-squared distribution
  integer(i4)              :: k          !! number of groups
  integer(i4)              :: n          !! total number of samples
  integer(i4)              :: ni         !! group size (assumes balanced)
  real(wp)   , allocatable :: x_flat(:)  !! flattened x array
  real(wp)   , allocatable :: ranks(:)   !! real ranks
  real(wp)   , allocatable :: r_sum(:)   !! sum of ranks per group
  integer(i4)              :: idx        !! group indexing
  integer(i4)              :: i, j       !! flexible integers

! ==== Instructions

  ! get sizes
  k = size(x, 2)
  ni = size(x, 1)
  n = k * ni

  ! allocate ranks
  allocate(x_flat(n))
  allocate(ranks(n))
  allocate(r_sum(k))

  ! flatten x
  x_flat = reshape(x, [n])

  ! compute ranks
  call s_utl_rank(x_flat, ranks)

  ! allocate and compute rank sums for each group
  r_sum = 0.0_wp
  idx = 1
  do j = 1, k
     r_sum(j) = sum( ranks(idx:idx + ni - 1) )
     idx = idx + ni
  end do

  ! compute H-statistic
  h = 0.0_wp
  do j = 1, k
     h = h + ( r_sum(j) ** 2 ) / real(ni, kind=wp)
  end do
  h = 12.0_wp * h / real(n * (n + 1), wp) - 3.0_wp * real(n + 1, kind=wp)

  ! degrees of freedom and p-value
  df = real(k - 1, kind=wp)
  p = f_dst_chi2_cdf_core(h, df, 0.0_wp, 1.0_wp, "right")

  ! deallocate
  deallocate(x_flat)
  deallocate(ranks)
  deallocate(r_sum)

end subroutine s_tst_kruskalwallis_core