Welch's t-test: Difference between revisions

Added 11l
(Added 11l)
Line 65:
 
The <math>\ln(\Gamma(x))</math>, or <code>lgammal(x)</code> function is necessary for the program to work with large <code>a</code> values, as [http://rosettacode.org/wiki/Gamma_function Gamma functions] can often return values larger than can be handled by <code>double</code> or <code>long double</code> data types. The <code>lgammal(x)</code> function is standard in <code>math.h</code> with C99 and C11 standards.
 
=={{header|11l}}==
{{trans|Python}}
 
<lang 11l>F betain(x, p, q)
I p <= 0 | q <= 0 | x < 0 | x > 1
X ValueError(0)
 
I x == 0 | x == 1
R x
 
V acu = 1e-15
V lnbeta = lgamma(p) + lgamma(q) - lgamma(p + q)
 
V xx = x
V cx = 1 - x
V pp = p
V qq = q
V indx = 0B
V psq = p + q
I p < psq * x
xx = 1 - x
cx = x
pp = q
qq = p
indx = 1B
 
V term = 1.0
V ai = 1.0
V value = 1.0
V ns = floor(qq + cx * psq)
V rx = xx / cx
V temp = qq - ai
I ns == 0
rx = xx
 
L
term *= temp * rx / (pp + ai)
value += term
temp = abs(term)
 
I temp <= acu & temp <= acu * value
value *= exp(pp * log(xx) + (qq - 1) * log(cx) - lnbeta) / pp
R I indx {1 - value} E value
 
ai++
I --ns >= 0
temp = qq - ai
I ns == 0
rx = xx
E
temp = psq
psq++
 
F welch_ttest(a1, a2)
V n1 = a1.len
V n2 = a2.len
I n1 <= 1 | n2 <= 1
X ValueError(0)
 
V mean1 = sum(a1) / n1
V mean2 = sum(a2) / n2
 
V var1 = sum(a1.map(x -> (x - @mean1) ^ 2)) / (n1 - 1)
V var2 = sum(a2.map(x -> (x - @mean2) ^ 2)) / (n2 - 1)
 
V t = (mean1 - mean2) / sqrt(var1 / n1 + var2 / n2)
V df = (var1 / n1 + var2 / n2) ^ 2 / (var1 ^ 2 / (n1 ^ 2 * (n1 - 1)) + var2 ^ 2 / (n2 ^ 2 * (n2 - 1)))
V p = betain(df / (t ^ 2 + df), df / 2, 1 / 2)
 
R (t, df, p)
 
V a1 = [Float(3), 4, 1, 2.1]
V a2 = [Float(490.2), 340, 433.9]
print(welch_ttest(a1, a2))</lang>
 
{{out}}
<pre>
(-9.5595, 2.00085, 0.0107516)
</pre>
 
=={{header|C}}==
1,481

edits