s_tst_ttest_2s_core Subroutine

pure subroutine s_tst_ttest_2s_core(x1, x2, t, df, p, eq_var, h1)

! ==== Description
!! The 2-sample t-test.

! ==== Declarations
  real(wp)        , intent(in)  :: x1(:)    !! x1 vector (samples)
  real(wp)        , intent(in)  :: x2(:)    !! x2 vector (samples)
  character(len=*), intent(in)  :: h1       !! \( H_{1} \) option: two (default), le, ge
  logical         , intent(in)  :: eq_var   !! true if equal variances assumed
  real(wp)        , intent(out) :: t        !! test statistic
  real(wp)        , intent(out) :: df       !! degrees of freedom
  real(wp)        , intent(out) :: p        !! p-value
  real(wp)                      :: x1bar    !! sample mean for x1
  real(wp)                      :: x2bar    !! sample mean for x2
  real(wp)                      :: s1       !! sample 1 standard deviation
  real(wp)                      :: s2       !! sample 2 standard deviation
  integer(i4)                   :: n1       !! sample size for x1
  integer(i4)                   :: n2       !! sample size for x2
  real(wp)                      :: s1n, s2n !! equation terms

! ==== Instructions

! ---- conduct test

  ! get means, sample sizes, and sample standard deviations (using n-1)
  x1bar = f_sts_mean_core(x1)
  x2bar = f_sts_mean_core(x2)
  n1 = size(x1)
  n2 = size(x2)
  s1 = sqrt( dot_product( (x1 - x1bar), (x1 - x1bar) ) / &
     & real( (n1-1), kind=wp ) )
  s2 = sqrt( dot_product( (x2 - x2bar), (x2 - x2bar) ) / &
     & real( (n2-1), kind=wp ) )

  ! get test statistic
  if (eq_var) then
     t = f_tst_ttest_2s_pooled_t(x1bar, x2bar, s1, s2, n1, n2)
  else
     t = f_tst_ttest_2s_welch_t(x1bar, x2bar, s1, s2, n1, n2)
  endif

  ! get degrees of freedom
  if (eq_var) then
     df = real( (n1 + n2 - 2), kind=wp )
  else
     s1n = (s1 * s1) / real( (n1), kind=wp )
     s2n = (s2 * s2) / real( (n2), kind=wp )
     df = ( (s1n + s2n) * (s1n + s2n) ) / ( &
        & (s1n * s1n) / real( (n1-1), kind=wp ) + &
        & (s2n * s2n) / real( (n2-1), kind=wp ) )
  endif

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

  contains

  ! --------------------------------------------------------------- !
  pure function f_tst_ttest_2s_welch_t(x1bar, x2bar, s1, s2, n1, n2) result(t)

     ! ==== Description
     !! Calculates the test statstic \( t \) for 2 sample t-test for unequal variances.
     ! TODO: Think about making elemental and public for batch processing

     ! ==== Declarations
     real(wp)   , intent(in) :: x1bar !! sample mean for x1
     real(wp)   , intent(in) :: x2bar !! sample mean for x2
     real(wp)   , intent(in) :: s1    !! sample 1 standard deviation
     real(wp)   , intent(in) :: s2    !! sample 2 standard deviation
     integer(i4), intent(in) :: n1    !! sample size for x1
     integer(i4), intent(in) :: n2    !! sample size for x2
     real(wp)                :: t     !! test statistic

     ! ==== Instructions
     t = (x1bar - x2bar) / ( sqrt( &
       & ( (s1 * s1) / real(n1, kind=wp)) + &
       & ( (s2 * s2) / real(n2, kind=wp)) ) )

  end function f_tst_ttest_2s_welch_t

  ! --------------------------------------------------------------- !
  pure function f_tst_ttest_2s_pooled_t(x1bar, x2bar, s1, s2, n1, n2) result(t)

     ! ==== Description
     !! Calculates the test statistic \( t \) for a two-sample t-test for equal variances.
     !! The function uses the pooled standard deviation.
     ! TODO: Think about making elemental and public for batch processing

     ! ==== Declarations
     real(wp)   , intent(in) :: x1bar !! sample mean for x1
     real(wp)   , intent(in) :: x2bar !! sample mean for x2
     real(wp)   , intent(in) :: s1    !! sample 1 standard deviation
     real(wp)   , intent(in) :: s2    !! sample 2 standard deviation
     integer(i4), intent(in) :: n1    !! sample size for x1
     integer(i4), intent(in) :: n2    !! sample size for x2
     real(wp)                :: t     !! test statistic
     real(wp)                :: sp    !! pooled variance

     ! ==== Instructions

     ! pooled standard deviation
     sp = sqrt( ( &
        & real(n1 - 1, kind=wp) * (s1 * s1) + &
        & real(n2 - 1, kind=wp) * (s2 * s2) ) / real(n1 + n2 - 2, kind=wp) )

     ! test statistic
     t = (x1bar - x2bar) / ( sp * sqrt( &
       & 1.0_wp / real(n1, kind=wp) + 1.0_wp / real(n2, kind=wp) ) )

  end function f_tst_ttest_2s_pooled_t

end subroutine s_tst_ttest_2s_core