Welch's t-test: Difference between revisions

</pre>

=={{header|Phix}}==

{{trans|Go}}

{{trans|Kotlin}}

<lang Phix>function mean(sequence a)

return sum(a) / length(a)

end function

function sv(sequence a)

integer la = length(a)

atom m := mean(a),

tot := 0

for i=1 to la do

atom d = a[i] - m

tot += d * d

end for

return tot / (la-1)

end function

function welch(sequence a, b)

integer la = length(a),

lb = length(b)

return (mean(a) - mean(b)) / sqrt(sv(a)/la+sv(b)/lb)

end function

function dof(sequence a, b)

integer la = length(a),

lb = length(b)

atom sva := sv(a),

svb := sv(b),

n := sva/la + svb/lb

return n * n / (sva*sva/(la*la*(la-1)) +

svb*svb/(lb*lb*(lb-1)))

end function

function f(atom r, v)

return power(r, v/2-1) / sqrt(1-r)

end function

function simpson0(integer n, atom high, v)

atom tot := 0,

dx0 := high / n,

x0 := dx0, x1, xmid, dx

tot += f(0,v) * dx0

tot += f(dx0*.5,v) * dx0 * 4

for i=1 to n-1 do

x1 := (i+1) * high / n

xmid := (x0 + x1) * .5

dx := x1 - x0

tot += f(x0,v) * dx * 2

tot += f(xmid,v) * dx * 4

x0 = x1

end for

return (tot + f(high,v)*dx0) / 6

end function

constant p = {

0.99999999999980993,

676.5203681218851,

-1259.1392167224028,

771.32342877765313,

-176.61502916214059,

12.507343278686905,

-0.13857109526572012,

9.9843695780195716e-6,

1.5056327351493116e-7

}

function gamma(atom d)

atom dd = d,

g = 7

if dd<0.5 then

return PI / (sin(PI*dd) * gamma(1-dd))

end if

dd -= 1

atom a = p[1],

t = dd + g + 0.5

for i=2 to length(p) do a += p[i] / (dd + i - 1) end for

return sqrt(2*PI) * power(t, dd + 0.5) * exp(-t) * a

end function

function lGamma(atom d)

return log(gamma(d))

end function

function pValue(sequence ab)

sequence {a, b} = ab

atom v := dof(a, b),

t := welch(a, b),

g1 := lGamma(v / 2),

g2 := lGamma(.5),

g3 := lGamma(v/2 + .5)

return simpson0(2000, v/(t*t+v), v) / exp(g1+g2-g3)

end function

constant tests = {{{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},

{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}},

{{17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8},

{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}},

{{19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0},

{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}},

{{30.02, 29.99, 30.11, 29.97, 30.01, 29.99},

{29.89, 29.93, 29.72, 29.98, 30.02, 29.98}},

{{3.0, 4.0, 1.0, 2.1},

{490.2, 340.0, 433.9}}

}

for i=1 to length(tests) do

?pValue(tests[i])

end for</lang>

{{out}}

<pre>

0.02137800146

0.1488416966

0.03597227103

0.09077332429

0.01075067374