s_tst_anova_1w_core Subroutine

pure subroutine s_tst_anova_1w_core(x, f, df_b, df_w, p)

! ==== Description
!! One-way ANOVA.

! ==== Declarations
  real(wp), intent(in)  :: x(:,:)    !! 2D array, each column is a group
  real(wp), intent(out) :: f         !! F-statistic
  real(wp), intent(out) :: p         !! p-value from F distribution
  real(wp), intent(out) :: df_b      !! degrees of freedom between groups
  real(wp), intent(out) :: df_w      !! degrees of freedom within groups
  integer(i4)           :: ni        !! number of elements in groups
  integer(i4)           :: n         !! number of elements in total
  integer(i4)           :: k         !! number of groups
  real(wp)              :: mu_t      !! grand/total mean
  real(wp)              :: mu_g      !! group mean
  real(wp)              :: ss_b      !! ss between groups
  real(wp)              :: ss_w      !! ss within groups
  real(wp), allocatable :: x_flat(:) !! flattened x
  integer(i4)           :: i, j      !! flexible integers

! ==== Instructions

  ! flatten all elements to compute grand mean and total n
  n = size(x)
  allocate(x_flat(n))
  x_flat = reshape(x, [n])

  ! get grand mean
  mu_t = f_sts_mean_core(x_flat)

  ! initialise sums
  ss_b = 0.0_wp
  ss_w = 0.0_wp

  ! get number of groups
  k = size(x, 2)

  ni = size(x, 1)
  do j = 1, k
     mu_g = f_sts_mean_core(x(:,j))
     ss_b = ss_b + real(ni, kind=wp) * (mu_g - mu_t) ** 2
     ss_w = ss_w + sum( (x(:,j) - mu_g) ** 2 )
  enddo

  ! degrees of freedom
  df_b = real(k - 1, kind=wp)
  df_w = real(n - k, kind=wp)

  ! calculate f statistics
  f = (ss_b / df_b) / (ss_w / df_w)

  ! get right tail p value with F distribution CDF procedure; use default loc, scale
  p = f_dst_f_cdf_core(f, df_b, df_w, 0.0_wp, 1.0_wp, "right")

end subroutine s_tst_anova_1w_core