Welch's t-test

func p_value (A, B) {
    [A.len, B.len].all { _ > 1 } || return 1

    var x̄_a = Math.avg(A...)
    var x̄_b = Math.avg(B...)

    var a_var = (A.map {|n| (x̄_a - n)**2 }.sum / A.end)
    var b_var = (B.map {|n| (x̄_b - n)**2 }.sum / B.end)

    (a_var && b_var) || return 1

    var Welsh_𝒕_statistic = ((x̄_a - x̄_b) / √(a_var/A.len + b_var/B.len))

    var DoF = ((a_var/A.len + b_var/B.len)**2 /
              ((a_var**2 / (A.len**3 - A.len**2)) + (b_var**2 / (B.len**3 - B.len**2))))

    var sa = (DoF/2 - 1)
    var x  = (DoF/(Welsh_𝒕_statistic**2 + DoF))
    var N  = 65355
    var h  = x/N

    var (sum1=0, sum2=0)

    ^N -> lazy.map { _ * h }.each { |i|
        sum1 += (((i + h/2) ** sa) / √(1 - (i + h/2)))
        sum2 += (( i        ** sa) / √(1 - (i      )))
    }

    (h/6 * (x**sa / √(1-x) + 4*sum1 + 2*sum2)) /
        (gamma(sa + 1) * √(Num.pi) / gamma(sa + 1.5))
}

# Testing
var tests = [
  %n<27.5 21.0 19.0 23.6 17.0 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19.0 21.7 21.4>,
  %n<27.1 22.0 20.8 23.4 23.4 23.5 25.8 22.0 24.8 20.2 21.9 22.1 22.9 20.5 24.4>,

  %n<17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8>,
  %n<21.5 22.8 21.0 23.0 21.6 23.6 22.5 20.7 23.4 21.8 20.7 21.7 21.5 22.5 23.6 21.5 22.5 23.5 21.5 21.8>,

  %n<19.8 20.4 19.6 17.8 18.5 18.9 18.3 18.9 19.5 22.0>,
  %n<28.2 26.6 20.1 23.3 25.2 22.1 17.7 27.6 20.6 13.7 23.2 17.5 20.6 18.0 23.9 21.6 24.3 20.4 24.0 13.2>,

  %n<30.02 29.99 30.11 29.97 30.01 29.99>,
  %n<29.89 29.93 29.72 29.98 30.02 29.98>,

  %n<3.0 4.0 1.0 2.1>,
  %n<490.2 340.0 433.9>
]

tests.each_slice(2, {|left, right|
    say p_value(left, right)
})

Output:

0.0213780014628670325061113281387220205111519317756
0.148841696605327985083613019511085971435711697961
0.0359722710297967180871367618538977446933248150651
0.0907733242856668878840956275523536083406692525656
0.0107515340333929755465323718028856669932912031012