Welch's t-test: Difference between revisions

</pre>

{{trans|Python}}

The above was a bit off on the fifth test, so I also tried this.<br>

using gamma() from [[Gamma_function#Phix]] (the one from above is probably also fine, but I didn't test that)

<lang Phix>function lgamma(atom d)

return log(gamma(d))

end function

function betain(atom x, p, q)

if p<=0 or q<=0 or x<0 or x>1 then ?9/0 end if

if x == 0 or x == 1 then return x end if

atom acu = 1e-15,

lnbeta = lgamma(p) + lgamma(q) - lgamma(p + q),

psq = p + q, cx = 1-x

bool indx = (p<psq*x)

if indx then

{cx,x,p,q} = {x,1-x,q,p}

end if

atom term = 1,

ai = 1,

val = 1,

ns = floor(q + cx*psq),

rx = iff(ns=0?x:x/cx),

temp = q - ai

while true do

term *= temp * rx / (p + ai)

val += term

temp = abs(term)

if temp<=acu and temp<=acu*val then

val *= exp(p*log(x) + (q-1)*log(cx) - lnbeta) / p

return iff(indx?1-val:val)

end if

ai += 1

ns -= 1

if ns>=0 then

temp = q - ai

if ns == 0 then

rx = x

end if

else

temp = psq

psq += 1

end if

end while

end function

function welch_ttest(sequence ab)

sequence {a, b} = ab

integer la = length(a),

lb = length(b)

atom ma = sum(a)/la,

mb = sum(b)/lb,

va = sum(sq_power(sq_sub(a,ma),2))/(la-1),

vb = sum(sq_power(sq_sub(b,mb),2))/(lb-1),

n = va/la + vb/lb,

t = (ma-mb)/sqrt(n),

df = (n*n) / (va*va/(la*la*(la-1)) + vb*vb/(lb*lb*(lb-1)))

return betain(df/(t*t+df), df/2, 1/2)

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

{{0.010268,0.000167,0.000167},

{0.159258,0.136278,0.122389}},

{{1.0/15,10.0/62.0},

{1.0/10,2/50.0}},

{{9/23.0,21/45.0,0/38.0},

{0/44.0,42/94.0,0/22.0}}},

correct = {0.021378001462867,

0.148841696605327,

0.0359722710297968,

0.090773324285671,

0.0107515611497845,

0.00339907162713746,

0.52726574965384,

0.545266866977794}

atom cerr = 0

for i=1 to length(tests) do

atom r = welch_ttest(tests[i])

?r

cerr += abs(r-correct[i])

end for

?{"cumulative error",cerr}</lang>

{{out}}

<pre>

0.02137800146

0.1488416966

0.03597227103

0.09077332429

0.01075156115

0.003399071627

0.5272657497

0.545266867

{"cumulative error",1.989380882e-14} -- (32 bit)

{"cumulative error",4.915115776e-15} -- (64-bit)