Welch's t-test: Difference between revisions

Added FreeBASIC
(Added Wren)
(Added FreeBASIC)
 
(13 intermediate revisions by 10 users not shown)
Line 12:
Your task is to discern whether or not the difference in means between the two sets is statistically significant and worth further investigation. P-values are significance tests to gauge the probability that the difference in means between two data sets is significant, or due to chance. A threshold level, alpha, is usually chosen, 0.01 or 0.05, where p-values below alpha are worth further investigation and p-values above alpha are considered not significant. The p-value is not considered a final test of significance, [http://www.nature.com/news/scientific-method-statistical-errors-1.14700 only whether the given variable should be given further consideration].
 
There is more than onone way of calculating the [[wp:Student's_t-test|t-statistic]], and you must choose which method is appropriate for you. Here we use [[wp:Welch's_t_test|Welch's t-test]], which assumes that the variances between the two sets <code>x</code> and <code>y</code> are not equal. Welch's t-test statistic can be computed:
 
<math>t \quad = \quad {\; \overline{X}_1 - \overline{X}_2 \; \over \sqrt{ \; {s_1^2 \over N_1} \; + \; {s_2^2 \over N_2} \quad }} </math>
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}}
 
<syntaxhighlight lang="11l">F betain(x, p, q)
I p <= 0 | q <= 0 | x < 0 | x > 1
X.throw 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.throw 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))</syntaxhighlight>
 
{{out}}
<pre>
(-9.5595, 2.00085, 0.0107516)
</pre>
 
=={{header|C}}==
Line 81 ⟶ 161:
This shows how pvalue can be calculated from any two arrays, using Welch's 2-sided t-test, which doesn't assume equal variance.
This is the equivalent of R's<code>t.test(vector1,vector2, alternative="two.sided", var.equal=FALSE)</code> and as such, it is compared against R's pvalues with the same vectors/arrays to show that the differences are very small (here 10^-14).
<langsyntaxhighlight Clang="c">#include <stdio.h>
#include <math.h>
#include <stdlib.h>
Line 348 ⟶ 428:
return 0;
}
</syntaxhighlight>
</lang>
 
{{out}}
Line 363 ⟶ 443:
'''If''' your computer does not have <code>lgammal</code>, add the following function before <code>main</code> and replace <code>lgammal</code> with <code>lngammal</code> in the <code>calculate_Pvalue</code> function:
 
<langsyntaxhighlight Clang="c">#include <stdio.h>
#include <math.h>
 
Line 383 ⟶ 463:
}
 
</syntaxhighlight>
</lang>
 
=={{header|Fortran}}==
Line 390 ⟶ 470:
Alternatively, the program shows the p-value computed using the IMSL '''BETAI''' function.
 
<langsyntaxhighlight lang="fortran">subroutine welch_ttest(n1, x1, n2, x2, t, df, p)
use tdf_int
implicit none
Line 415 ⟶ 495:
print *, t, df, p
print *, betai(df / (t**2 + df), 0.5d0 * df, 0.5d0)
end program</langsyntaxhighlight>
 
'''Output'''
Line 425 ⟶ 505:
With Absoft Pro Fortran, compile with <code>af90 -m64 pvalue.f90 %SLATEC_LINK%</code>.
 
<langsyntaxhighlight lang="fortran">subroutine welch_ttest(n1, x1, n2, x2, t, df, p)
implicit none
integer :: n1, n2
Line 450 ⟶ 530:
call welch_ttest(4, x, 3, y, t, df, p)
print *, t, df, p
end program</langsyntaxhighlight>
 
'''Output'''
 
<pre> -9.55949772193266 2.00085234885628 1.075156114978449E-002</pre>
 
=== Using GSL ===
{{works with|Fortran|95}}
 
{{libheader|GNU Scientific Library}}
 
Instead of implementing the t-distribution by ourselves, we bind to GNU Scientific Library:
<syntaxhighlight lang="fortran">module t_test_m
 
use, intrinsic :: iso_c_binding
use, intrinsic :: iso_fortran_env, only: wp => real64
implicit none
private
 
public :: t_test, wp
 
interface
function gsl_cdf_tdist_p(x, nu) bind(c, name='gsl_cdf_tdist_P')
import
real(c_double), value :: x
real(c_double), value :: nu
real(c_double) :: gsl_cdf_tdist_p
end function gsl_cdf_tdist_p
end interface
 
contains
 
!> Welch T test
impure subroutine t_test(x, y, p, t, df)
real(wp), intent(in) :: x(:), y(:)
real(wp), intent(out) :: p !! p-value
real(wp), intent(out) :: t !! T value
real(wp), intent(out) :: df !! degrees of freedom
integer :: n1, n2
real(wp) :: m1, m2, v1, v2
 
n1 = size(x)
n2 = size(y)
m1 = sum(x)/n1
m2 = sum(y)/n2
v1 = sum((x - m1)**2)/(n1 - 1)
v2 = sum((y - m2)**2)/(n2 - 1)
 
t = (m1 - m2)/sqrt(v1/n1 + v2/n2)
df = (v1/n1 + v2/n2)**2/(v1**2/(n1**2*(n1 - 1)) + v2**2/(n2**2*(n2 - 1)))
p = 2*gsl_cdf_tdist_p(-abs(t), df)
 
end subroutine t_test
 
end module t_test_m
 
program main
use t_test_m, only: t_test, wp
implicit none
real(wp) :: x(4) = [3.0_wp, 4.0_wp, 1.0_wp, 2.1_wp]
real(wp) :: y(3) = [490.2_wp, 340.0_wp, 433.9_wp]
real(wp) :: t, df, p
 
call t_test(x, y, p, t, df)
print *, t, df, p
 
end program main</syntaxhighlight>
 
'''Output'''
 
<pre> -9.5594977219326580 2.0008523488562844 1.0751561149784494E-002</pre>
 
=={{header|FreeBASIC}}==
===Using Betain===
{{trans|11l}}
<syntaxhighlight lang="vbnet">#include "crt\math.bi"
 
Function betain(x As Double, p As Double, q As Double) As Double
If p <= 0 Or q <= 0 Or x < 0 Or x > 1 Then
Print "ValueError"
End
End If
If x = 0 Or x = 1 Then Return x
Dim As Double acu = 1e-15
'Dim As Double lnbeta = LogGamma(p) + LogGamma(q) - LogGamma(p + q)
Dim As Double lnbeta = lGamma(p) + lGamma(q) - lGamma(p + q)
Dim As Double xx = x
Dim As Double cx = 1 - x
Dim As Double pp = p
Dim As Double qq = q
Dim As Integer indx = 0
Dim As Double psq = p + q
If p < psq * x Then
xx = 1 - x
cx = x
pp = q
qq = p
indx = 1
End If
Dim As Double term = 1.0
Dim As Double ai = 1.0
Dim As Double value = 1.0
Dim As Integer ns = Int(qq + cx * psq)
Dim As Double rx = xx / cx
Dim As Double temp = qq - ai
If ns = 0 Then rx = xx
Do
term *= temp * rx / (pp + ai)
value += term
temp = Abs(term)
If temp <= acu And temp <= acu * value Then
value *= Exp(pp * Log(xx) + (qq - 1) * Log(cx) - lnbeta) / pp
Return Iif(indx, 1 - value, value)
End If
ai += 1
If ns > 0 Then
ns -= 1
temp = qq - ai
If ns = 0 Then
rx = xx
Else
temp = psq
psq += 1
End If
End If
Loop
End Function
 
Sub welch_ttest(a1() As Double, a2() As Double, Byref t As Double, Byref df As Double, Byref p As Double)
Dim As Integer n1 = Ubound(a1) + 1
Dim As Integer n2 = Ubound(a2) + 1
If n1 <= 1 Or n2 <= 1 Then
Print "ValueError"
End
End If
 
Dim As Double mean1 = 0
For i As Integer = 0 To n1 - 1
mean1 += a1(i)
Next i
mean1 /= n1
 
Dim As Double mean2 = 0
For i As Integer = 0 To n2 - 1
mean2 += a2(i)
Next i
mean2 /= n2
 
Dim As Double var1 = 0
For i As Integer = 0 To n1 - 1
var1 += (a1(i) - mean1) ^ 2
Next i
var1 /= (n1 - 1)
 
Dim As Double var2 = 0
For i As Integer = 0 To n2 - 1
var2 += (a2(i) - mean2) ^ 2
Next i
var2 /= (n2 - 1)
 
t = (mean1 - mean2) / Sqr(var1 / n1 + var2 / n2)
df = (var1 / n1 + var2 / n2) ^ 2 / (var1 ^ 2 / (n1 ^ 2 * (n1 - 1)) + var2 ^ 2 / (n2 ^ 2 * (n2 - 1)))
p = betain(df / (t ^ 2 + df), df / 2, 1 / 2)
End Sub
 
Dim As Double a1(3) = {3, 4, 1, 2.1}
Dim As Double a2(2) = {490.2, 340, 433.9}
Dim As Double t, df, p
welch_ttest(a1(), a2(), t, df, p)
Print " t: "; t
Print "df: "; df
Print " p: "; p
 
Sleep</syntaxhighlight>
{{out}}
<pre> t: -9.559497721932658
df: 2.000852348856284
p: 0.01075155600241868</pre>
 
=={{header|Go}}==
<langsyntaxhighlight lang="go">package main
 
import (
Line 548 ⟶ 808:
func(r float64) float64 { return math.Pow(r, ν/2-1) / math.Sqrt(1-r) }) /
math.Exp(g1+g2-g3)
}</langsyntaxhighlight>
{{out}}
<pre>
Line 562 ⟶ 822:
Implementation:
 
<langsyntaxhighlight Jlang="j">integrate=: adverb define
'a b steps'=. 3{.y,128
size=. (b - a)%steps
Line 585 ⟶ 845:
hi=. v%(t^2)+v
(F f. simpson integrate lo,hi) % 0.5 B v%2
)</langsyntaxhighlight>
 
<code>integrate</code> and <code>simpson</code> are from the [[Numerical_integration#J|Numerical integration]] task.
Line 602 ⟶ 862:
 
Data for task examples:
<langsyntaxhighlight Jlang="j">d1=: 27.5 21 19 23.6 17 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19 21.7 21.4
d2=: 27.1 22 20.8 23.4 23.4 23.5 25.8 22 24.8 20.2 21.9 22.1 22.9 20.5 24.4
d3=: 17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8
Line 611 ⟶ 871:
d8=: 29.89 29.93 29.72 29.98 30.02 29.98
d9=: 3 4 1 2.1
da=: 490.2 340 433.9</langsyntaxhighlight>
 
Task examples:
<langsyntaxhighlight Jlang="j"> d1 p2_tail d2
0.021378
d3 p2_tail d4
Line 623 ⟶ 883:
0.0907733
d9 p2_tail da
0.0107377</langsyntaxhighlight>
 
=={{header|Java}}==
Using the '''[http://commons.apache.org/proper/commons-math/ Apache Commons Mathematics Library]'''.
<langsyntaxhighlight lang="java">import org.apache.commons.math3.distribution.TDistribution;
 
public class WelchTTest {
Line 676 ⟶ 936:
System.out.println("p = " + res[2]);
}
}</langsyntaxhighlight>
 
'''Result'''
Line 685 ⟶ 945:
df = 2.0008523488562844
p = 0.010751561149784485</pre>
 
=={{header|jq}}==
# {{trans|Wren}}
{{works with|jq}}
'''Works with gojq, the Go implementation of jq'''
 
Notice how jq supports the closure, f, in the same way as Wren.
 
jq's `lgamma` returns the natural logarithm of the absolute value of the gamma function of x.
<syntaxhighlight lang="jq">def mean: add / length;
 
# Sample variance using division by (length-1)
def variance:
mean as $m
| (reduce .[] as $x (0; . + (($x - $m) | .*.))) / (length-1) ;
 
def welch(a; b):
((a|mean) - (b|mean)) /
(((a|variance/length) + (b|variance/length)) | sqrt) ;
 
def dof(a; b):
(a|variance) as $sva
| (b|variance) as $svb
| (a|length) as $la
| (b|length) as $lb
| ($sva/$la + $svb/$lb) as $n
| $n * $n / ($sva*$sva/($la*$la*($la-1)) + $svb*$svb/($lb*$lb*($lb-1))) ;
 
def simpson0(nf; upper; filter):
(upper/nf) as $dx0
| {sum: (( (0|filter) + ($dx0 * 0.5|filter) * 4) * $dx0),
x0: $dx0 }
| reduce range(1; nf) as $i (.;
( ($i + 1) * upper / nf ) as $x1
| ((.x0 + $x1) * 0.5) as $xmid
| ($x1 - .x0) as $dx
| .sum = .sum + ((.x0|filter)*2 + ($xmid|filter)*4) * $dx
| .x0 = $x1)
| (.sum + (upper|filter)*$dx0) / 6 ;
 
def pValue(a; b):
dof(a; b) as $nu
| def f:
. as $r
| pow($r; ($nu/2) - 1) / ((1 - $r)|sqrt);
 
welch(a; b) as $t
| (($nu/2)|lgamma) as $g1
| (0.5|lgamma) as $g2
| (($nu/2 + 0.5)|lgamma) as $g3
| simpson0(2000; $nu/($t*$t + $nu); f) / (($g1 + $g2 - $g3)|exp) ;
def d1: [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];
def d2: [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];
def d3: [17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8];
def d4: [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];
def d5: [19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0];
def d6: [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];
def d7: [30.02, 29.99, 30.11, 29.97, 30.01, 29.99];
def d8: [29.89, 29.93, 29.72, 29.98, 30.02, 29.98];
def x : [3.0, 4.0, 1.0, 2.1];
def y : [490.2, 340.0, 433.9];
 
pValue(d1; d2),
pValue(d3; d4),
pValue(d5; d6),
pValue(d7; d8),
pValue(x; y)</syntaxhighlight>
{{out}}
<pre>
0.02137800146288292
0.1488416966053347
0.03597227102982764
0.09077332428566065
0.010750673736239608
</pre>
 
 
=={{header|Julia}}==
{{works with|Julia|0.6}}
 
<langsyntaxhighlight lang="julia">using HypothesisTests
 
d1 = [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]
Line 709 ⟶ 1,048:
ttest = UnequalVarianceTTest(y1, y2)
println("\nData:\n y1 = $y1\n y2 = $y2\nP-value for unequal variance TTest: ", round(pvalue(ttest), 4))
end</langsyntaxhighlight>
 
{{out}}
Line 741 ⟶ 1,080:
=={{header|Kotlin}}==
This program brings in code from other tasks for gamma functions and integration by Simpson's rule as Kotlin doesn't have these built-in:
<langsyntaxhighlight lang="scala">// version 1.1.4-3
 
typealias Func = (Double) -> Double
Line 846 ⟶ 1,185:
println(f.format(p2Tail(da7, da8)))
println(f.format(p2Tail(x, y)))
}</langsyntaxhighlight>
 
{{out}}
Line 856 ⟶ 1,195:
0.010751
</pre>
 
=={{header|Nim}}==
{{trans|Kotlin}}
<syntaxhighlight lang="nim">import math, stats, strutils, sugar
 
func sqr(f: float): float = f * f
 
func degreesFreedom(da1, da2: openArray[float]): float =
let s1 = varianceS(da1)
let s2 = varianceS(da2)
let n1 = da1.len.toFloat
let n2 = da2.len.toFloat
let n = sqr(s1 / n1 + s2 / n2)
let d = sqr(s1) / (n1 * n1 * (n1 - 1)) + sqr(s2) / (n2 * n2 * (n2 - 1))
result = n / d
 
func welch(da1, da2: openArray[float]): float =
let f = varianceS(da1) / da1.len.toFloat + varianceS(da2) / da2.len.toFloat
result = (mean(da1) - mean(da2)) / sqrt(f)
 
func simpson(a, b: float; n: int; f: float -> float): float =
let h = (b - a) / n.toFloat
var sum = 0.0
for i in 0..<n:
let x = a + i.toFloat * h
sum += (f(x) + 4 * f(x + h / 2) + f(x + h)) / 6
result = sum * h
 
func p2Tail(da1, da2: openArray[float]): float =
let ν = degreesFreedom(da1, da2)
let t = welch(da1, da2)
let g = exp(lGamma(ν / 2) + lGamma(0.5) - lGamma(ν / 2 + 0.5))
let b = ν / (t * t + ν)
proc f(r: float): float = pow(r, ν / 2 - 1) / sqrt(1 - r)
result = simpson(0, b, 10000, f) / g # n = 10000 seems more than enough here.
 
 
when isMainModule:
 
const
Da1 = [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]
Da2 = [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]
Da3 = [17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8]
Da4 = [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]
Da5 = [19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0]
Da6 = [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]
Da7 = [30.02, 29.99, 30.11, 29.97, 30.01, 29.99]
Da8 = [29.89, 29.93, 29.72, 29.98, 30.02, 29.98]
 
X = [3.0, 4.0, 1.0, 2.1]
Y = [490.2, 340.0, 433.9]
 
echo p2Tail(Da1, Da2).formatFloat(ffDecimal, 6)
echo p2Tail(Da3, Da4).formatFloat(ffDecimal, 6)
echo p2Tail(Da5, Da6).formatFloat(ffDecimal, 6)
echo p2Tail(Da7, Da8).formatFloat(ffDecimal, 6)
echo p2Tail(X, Y).formatFloat(ffDecimal, 6)</syntaxhighlight>
 
{{out}}
<pre>0.021378
0.148842
0.035972
0.090773
0.010751</pre>
 
=={{header|Maple}}==
 
<syntaxhighlight lang="maple">WelschTTest:=proc(x::list(numeric),y::list(numeric))
uses Statistics;
local n1:=nops(x),n2:=nops(y),
m1:=Mean(x),m2:=Mean(y),
v1:=Variance(x),v2:=Variance(y),
t,nu,p;
t:=(m1-m2)/sqrt(v1/n1+v2/n2);
nu:=(v1/n1+v2/n2)^2/(v1^2/(n1^2*(n1-1))+v2^2/(n2^2*(n2-1)));
p:=2*CDF(StudentTDistribution(nu),-abs(t));
t,nu,p
end proc:
 
x:=[3,4,1,2.1]:
y:=[490.2,340,433.9]:
WelschTTest(x,y);
# -9.55949772193266, 2.00085234885628, 0.0107515611497845</syntaxhighlight>
 
=={{header|Octave}}==
{{trans|Stata}}
<langsyntaxhighlight lang="octave">x = [3.0,4.0,1.0,2.1];
y = [490.2,340.0,433.9];
n1 = length(x);
Line 872 ⟶ 1,298:
ans =
 
-9.559498 2.000852 0.010752</langsyntaxhighlight>
 
=={{header|PARI/GP}}==
<langsyntaxhighlight lang="parigp">B2(x,y)=exp(lngamma(x)+lngamma(y)-lngamma(x+y))
B3(x,a,b)=a--;b--;intnum(r=0,x,r^a*(1-r)^b)
Welch2(u,v)=my(m1=vecsum(u)/#u, m2=vecsum(v)/#v, v1=var(u,m1), v2=var(v,m2), s=v1/#u+v2/#v, t=(m1-m2)/sqrt(s), nu=s^2/(v1^2/#u^2/(#u-1)+v2^2/#v^2/(#v-1))); B3(nu/(t^2+nu),nu/2,1/2)/B2(nu/2,1/2);
Welch2([3,4,1,2.1], [490.2,340,433.9])</langsyntaxhighlight>
{{out}}
<pre>%1 = 0.010751561149784496723954539777213062928</pre>
Line 886 ⟶ 1,312:
Uses Math::AnyNum for gamma and pi. It is possible to use some other modules (e.g. Math::Cephes) if Math::AnyNum has problematic dependencies.
{{trans|Sidef}}
<langsyntaxhighlight lang="perl">use utf8;
use List::Util qw(sum);
use Math::AnyNum qw(gamma pi);
 
sub p_value :prototype($$) {
my ($A, $B) = @_;
 
Line 946 ⟶ 1,372:
my ($left, $right) = splice(@tests, 0, 2);
print p_value($left, $right), "\n";
}</langsyntaxhighlight>
{{out}}
<pre>
Line 959 ⟶ 1,385:
We use a slightly more accurate lgamma than the C code. Note that Perl can be compiled with different underlying floating point representations -- double, long double, or quad double.
{{trans|C}}
<langsyntaxhighlight lang="perl">use strict;
use warnings;
use List::Util 'sum';
Line 1,086 ⟶ 1,512:
printf("p-value = %.14g\n",$pvalue);
}
printf("cumulative error is %g\n", $error);</langsyntaxhighlight>
{{out}}
<pre>p-value = 0.021378001462867
Line 1,101 ⟶ 1,527:
{{trans|Go}}
{{trans|Kotlin}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<lang Phix>function mean(sequence a)
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
return sum(a) / length(a)
<span style="color: #008080;">function</span> <span style="color: #000000;">mean</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function sv(sequence a)
integer la = length(a)
<span style="color: #008080;">function</span> <span style="color: #000000;">sv</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
atom m := mean(a),
<span style="color: #004080;">integer</span> <span style="color: #000000;">la</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
tot := 0
<span style="color: #004080;">atom</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">mean</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span>
for i=1 to la do
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">0</span>
atom d = a[i] - m
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">la</span> <span style="color: #008080;">do</span>
tot += d * d
<span style="color: #004080;">atom</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">m</span>
end for
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">d</span>
return tot / (la-1)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">la</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function welch(sequence a, b)
integer la = length(a),
<span style="color: #008080;">function</span> <span style="color: #000000;">welch</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>
lb = length(b)
<span style="color: #004080;">integer</span> <span style="color: #000000;">la</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span>
return (mean(a) - mean(b)) / sqrt(sv(a)/la+sv(b)/lb)
<span style="color: #000000;">lb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">return</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">mean</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">mean</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">))</span> <span style="color: #0000FF;">/</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">la</span><span style="color: #0000FF;">+</span><span style="color: #000000;">sv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function dof(sequence a, b)
integer la = length(a),
<span style="color: #008080;">function</span> <span style="color: #000000;">dof</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>
lb = length(b)
<span style="color: #004080;">integer</span> <span style="color: #000000;">la</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span>
atom sva := sv(a),
<span style="color: #000000;">lb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>
svb := sv(b),
<span style="color: #004080;">atom</span> <span style="color: #000000;">sva</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">sv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span>
n := sva/la + svb/lb
<span style="color: #000000;">svb</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">sv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">),</span>
return n * n / (sva*sva/(la*la*(la-1)) +
<span style="color: #000000;">n</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">sva</span><span style="color: #0000FF;">/</span><span style="color: #000000;">la</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">svb</span><span style="color: #0000FF;">/</span><span style="color: #000000;">lb</span>
svb*svb/(lb*lb*(lb-1)))
<span style="color: #008080;">return</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">sva</span><span style="color: #0000FF;">*</span><span style="color: #000000;">sva</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">la</span><span style="color: #0000FF;">*</span><span style="color: #000000;">la</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">la</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span> <span style="color: #0000FF;">+</span>
end function
<span style="color: #000000;">svb</span><span style="color: #0000FF;">*</span><span style="color: #000000;">svb</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function f(atom r, v)
return power(r, v/2-1) / sqrt(1-r)
<span style="color: #008080;">function</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">return</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function simpson0(integer n, atom high, v)
atom tot := 0,
<span style="color: #008080;">function</span> <span style="color: #000000;">simpson0</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">high</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
dx0 := high / n,
<span style="color: #004080;">atom</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span>
x0 := dx0, x1, xmid, dx
<span style="color: #000000;">dx0</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">high</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span>
tot += f(0,v) * dx0
<span style="color: #000000;">x0</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">dx0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">xmid</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dx</span>
tot += f(dx0*.5,v) * dx0 * 4
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">dx0</span>
for i=1 to n-1 do
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dx0</span><span style="color: #0000FF;">*.</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">dx0</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">4</span>
x1 := (i+1) * high / n
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
xmid := (x0 + x1) * .5
<span style="color: #000000;">x1</span> <span style="color: #0000FF;">:=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">high</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">n</span>
dx := x1 - x0
<span style="color: #000000;">xmid</span> <span style="color: #0000FF;">:=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">x0</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">x1</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #0000FF;">.</span><span style="color: #000000;">5</span>
tot += f(x0,v) * dx * 2
<span style="color: #000000;">dx</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">x1</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">x0</span>
tot += f(xmid,v) * dx * 4
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">dx</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">2</span>
x0 = x1
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xmid</span><span style="color: #0000FF;">,</span><span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">dx</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">4</span>
end for
<span style="color: #000000;">x0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x1</span>
return (tot + f(high,v)*dx0) / 6
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">tot</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">high</span><span style="color: #0000FF;">,</span><span style="color: #000000;">v</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">dx0</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">6</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
constant p = {
0.99999999999980993,
<span style="color: #008080;">constant</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span>
676.5203681218851,
<span style="color: #000000;">0.99999999999980993</span><span style="color: #0000FF;">,</span>
-1259.1392167224028,
<span style="color: #000000;">676.5203681218851</span><span style="color: #0000FF;">,</span>
771.32342877765313,
<span style="color: #0000FF;">-</span><span style="color: #000000;">1259.1392167224028</span><span style="color: #0000FF;">,</span>
-176.61502916214059,
<span style="color: #000000;">771.32342877765313</span><span style="color: #0000FF;">,</span>
12.507343278686905,
<span style="color: #0000FF;">-</span><span style="color: #000000;">176.61502916214059</span><span style="color: #0000FF;">,</span>
-0.13857109526572012,
<span style="color: #000000;">12.507343278686905</span><span style="color: #0000FF;">,</span>
9.9843695780195716e-6,
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.13857109526572012</span><span style="color: #0000FF;">,</span>
1.5056327351493116e-7
<span style="color: #000000;">9.9843695780195716e-6</span><span style="color: #0000FF;">,</span>
}
<span style="color: #000000;">1.5056327351493116e-7</span>
 
<span style="color: #0000FF;">}</span>
function gamma(atom d)
atom dd = d,
<span style="color: #008080;">function</span> <span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
g = 7
<span style="color: #004080;">atom</span> <span style="color: #000000;">dd</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">,</span>
if dd<0.5 then
<span style="color: #000000;">g</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">7</span>
return PI / (sin(PI*dd) * gamma(1-dd))
<span style="color: #008080;">if</span> <span style="color: #000000;">dd</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0.5</span> <span style="color: #008080;">then</span>
end if
<span style="color: #008080;">return</span> <span style="color: #004600;">PI</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">*</span><span style="color: #000000;">dd</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">dd</span><span style="color: #0000FF;">))</span>
dd -= 1
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
atom a = p[1],
<span style="color: #000000;">dd</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">1</span>
t = dd + g + 0.5
<span style="color: #004080;">atom</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span>
for i=2 to length(p) do a += p[i] / (dd + i - 1) end for
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">dd</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">g</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">0.5</span>
return sqrt(2*PI) * power(t, dd + 0.5) * exp(-t) * a
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">dd</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dd</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">a</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function lGamma(atom d)
return log(gamma(d))
<span style="color: #008080;">function</span> <span style="color: #000000;">lGamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">return</span> <span style="color: #7060A8;">log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function pValue(sequence ab)
sequence {a, b} = ab
<span style="color: #008080;">function</span> <span style="color: #000000;">pValue</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">ab</span><span style="color: #0000FF;">)</span>
atom v := dof(a, b),
<span style="color: #004080;">sequence</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ab</span>
t := welch(a, b),
<span style="color: #004080;">atom</span> <span style="color: #000000;">v</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">dof</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">),</span>
g1 := lGamma(v / 2),
<span style="color: #000000;">t</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">welch</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">),</span>
g2 := lGamma(.5),
<span style="color: #000000;">g1</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">lGamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">v</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">),</span>
g3 := lGamma(v/2 + .5)
<span style="color: #000000;">g2</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">lGamma</span><span style="color: #0000FF;">(.</span><span style="color: #000000;">5</span><span style="color: #0000FF;">),</span>
return simpson0(2000, v/(t*t+v), v) / exp(g1+g2-g3)
<span style="color: #000000;">g3</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">lGamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">v</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span> <span style="color: #0000FF;">+</span> <span style="color: #0000FF;">.</span><span style="color: #000000;">5</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">simpson0</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">*</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">v</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">g2</span><span style="color: #0000FF;">-</span><span style="color: #000000;">g3</span><span style="color: #0000FF;">)</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
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}},
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span><span style="color: #000000;">27.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">16.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.4</span><span style="color: #0000FF;">},</span>
{{17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8},
<span style="color: #0000FF;">{</span><span style="color: #000000;">27.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">25.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.4</span><span style="color: #0000FF;">}},</span>
{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}},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">17.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">14.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.8</span><span style="color: #0000FF;">},</span>
{{19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0},
<span style="color: #0000FF;">{</span><span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.8</span><span style="color: #0000FF;">}},</span>
{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}},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">19.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.0</span><span style="color: #0000FF;">},</span>
{{30.02, 29.99, 30.11, 29.97, 30.01, 29.99},
<span style="color: #0000FF;">{</span><span style="color: #000000;">28.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">26.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">25.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">27.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">13.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">13.2</span><span style="color: #0000FF;">}},</span>
{29.89, 29.93, 29.72, 29.98, 30.02, 29.98}},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">30.02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.99</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">30.11</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.97</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">30.01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.99</span><span style="color: #0000FF;">},</span>
{{3.0, 4.0, 1.0, 2.1},
<span style="color: #0000FF;">{</span><span style="color: #000000;">29.89</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.93</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.72</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.98</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">30.02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.98</span><span style="color: #0000FF;">}},</span>
{490.2, 340.0, 433.9}}
<span style="color: #0000FF;">{{</span><span style="color: #000000;">3.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.1</span><span style="color: #0000FF;">},</span>
}
<span style="color: #0000FF;">{</span><span style="color: #000000;">490.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">340.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">433.9</span><span style="color: #0000FF;">}}</span>
<span style="color: #0000FF;">}</span>
for i=1 to length(tests) do
?pValue(tests[i])
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end for</lang>
<span style="color: #0000FF;">?</span><span style="color: #000000;">pValue</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
0.0213780015
0.02137800146
0.1488416966
0.035972271
0.03597227103
0.0907733243
0.09077332429
0.0107506737
0.01075067374
</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)
<!--<syntaxhighlight lang="phix">(phixonline)-->
<lang Phix>function lgamma(atom d)
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
return log(gamma(d))
<span style="color: #000080;font-style:italic;">--&lt;copy of gamma from Gamma_function#Phix&gt;</span>
end function
<span style="color: #004080;">sequence</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">12</span><span style="color: #0000FF;">)</span>
 
function betain(atom x, p, q)
<span style="color: #008080;">function</span> <span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
if p<=0 or q<=0 or x<0 or x>1 then ?9/0 end if
<span style="color: #004080;">atom</span> <span style="color: #000000;">accm</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
if x == 0 or x == 1 then return x end if
<span style="color: #008080;">if</span> <span style="color: #000000;">accm</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">accm</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">accm</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">k1_factrl</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> <span style="color: #000080;font-style:italic;">-- (k - 1)!*(-1)^k with 0!==1</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">12</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">13</span><span style="color: #0000FF;">-</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)*</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">13</span><span style="color: #0000FF;">-</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1.5</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">k1_factrl</span>
<span style="color: #000000;">k1_factrl</span> <span style="color: #0000FF;">*=</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">12</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">accm</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]/(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">accm</span> <span style="color: #0000FF;">*=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(-(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">12</span><span style="color: #0000FF;">))*</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">12</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- Gamma(z+1)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">accm</span><span style="color: #0000FF;">/</span><span style="color: #000000;">z</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000080;font-style:italic;">--&lt;/copy of gamma&gt;</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">lgamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
atom acu = 1e-15,
<span style="color: #008080;">return</span> <span style="color: #7060A8;">log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">))</span>
lnbeta = lgamma(p) + lgamma(q) - lgamma(p + q),
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
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
<span style="color: #008080;">function</span> <span style="color: #000000;">betain</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">)</span>
while true do
<span style="color: #008080;">if</span> <span style="color: #000000;">p</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">q</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">x</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
term *= temp * rx / (p + ai)
<span style="color: #008080;">if</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">==</span> <span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">==</span> <span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
val += term
temp = abs(term)
<span style="color: #004080;">atom</span> <span style="color: #000000;">acu</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e-15</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">lnbeta</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lgamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">lgamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">lgamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">psq</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">x</span>
<span style="color: #004080;">bool</span> <span style="color: #000000;">indx</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;"><</span><span style="color: #000000;">psq</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">indx</span> <span style="color: #008080;">then</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">cx</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">term</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">ai</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">val</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">ns</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">cx</span><span style="color: #0000FF;">*</span><span style="color: #000000;">psq</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">rx</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ns</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">?</span><span style="color: #000000;">x</span><span style="color: #0000FF;">:</span><span style="color: #000000;">x</span><span style="color: #0000FF;">/</span><span style="color: #000000;">cx</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">temp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">q</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">ai</span>
<span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">term</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">temp</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">rx</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">p</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">ai</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">val</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">term</span>
<span style="color: #000000;">temp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">term</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">temp</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">acu</span> <span style="color: #008080;">and</span> <span style="color: #000000;">temp</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">acu</span><span style="color: #0000FF;">*</span><span style="color: #000000;">val</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">val</span> <span style="color: #0000FF;">*=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)*</span><span style="color: #7060A8;">log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cx</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">lnbeta</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">p</span>
<span style="color: #008080;">return</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">indx</span><span style="color: #0000FF;">?</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">val</span><span style="color: #0000FF;">:</span><span style="color: #000000;">val</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">ai</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
if temp<=acu and temp<=acu*val then
<span style="color: #000000;">ns</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">1</span>
val *= exp(p*log(x) + (q-1)*log(cx) - lnbeta) / p
<span style="color: #008080;">if</span> <span style="color: #000000;">ns</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
return iff(indx?1-val:val)
<span style="color: #000000;">temp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">q</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">ai</span>
end if
<span style="color: #008080;">if</span> <span style="color: #000000;">ns</span> <span style="color: #0000FF;">==</span> <span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
 
<span style="color: #000000;">rx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span>
ai += 1
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
ns -= 1
if ns><span style=0"color: then#008080;">else</span>
<span style="color: #000000;">temp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">psq</span>
temp = q - ai
<span style="color: #000000;">psq</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
if ns == 0 then
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
rx = x
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
else
temp = psq
<span style="color: #008080;">function</span> <span style="color: #000000;">welch_ttest</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">ab</span><span style="color: #0000FF;">)</span>
psq += 1
<span style="color: #004080;">sequence</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ab</span>
end if
<span style="color: #004080;">integer</span> <span style="color: #000000;">la</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span>
end while
<span style="color: #000000;">lb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #004080;">atom</span> <span style="color: #000000;">ma</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">la</span><span style="color: #0000FF;">,</span>
 
<span style="color: #000000;">mb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">,</span>
function welch_ttest(sequence ab)
<span style="color: #000000;">va</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ma</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))/(</span><span style="color: #000000;">la</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span>
sequence {a, b} = ab
<span style="color: #000000;">vb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mb</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))/(</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span>
integer la = length(a),
<span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">va</span><span style="color: #0000FF;">/</span><span style="color: #000000;">la</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">vb</span><span style="color: #0000FF;">/</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">,</span>
lb = length(b)
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">ma</span><span style="color: #0000FF;">-</span><span style="color: #000000;">mb</span><span style="color: #0000FF;">)/</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span>
atom ma = sum(a)/la,
<span style="color: #000000;">df</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">*</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">va</span><span style="color: #0000FF;">*</span><span style="color: #000000;">va</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">la</span><span style="color: #0000FF;">*</span><span style="color: #000000;">la</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">la</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">vb</span><span style="color: #0000FF;">*</span><span style="color: #000000;">vb</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">lb</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)))</span>
mb = sum(b)/lb,
<span style="color: #008080;">return</span> <span style="color: #000000;">betain</span><span style="color: #0000FF;">(</span><span style="color: #000000;">df</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">*</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">df</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">df</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
va = sum(sq_power(sq_sub(a,ma),2))/(la-1),
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
vb = sum(sq_power(sq_sub(b,mb),2))/(lb-1),
n = va/la + vb/lb,
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span><span style="color: #000000;">27.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">16.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.4</span><span style="color: #0000FF;">},</span>
t = (ma-mb)/sqrt(n),
<span style="color: #0000FF;">{</span><span style="color: #000000;">27.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">25.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.4</span><span style="color: #0000FF;">}},</span>
df = (n*n) / (va*va/(la*la*(la-1)) + vb*vb/(lb*lb*(lb-1)))
<span style="color: #0000FF;">{{</span><span style="color: #000000;">17.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">14.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.8</span><span style="color: #0000FF;">},</span>
return betain(df/(t*t+df), df/2, 1/2)
<span style="color: #0000FF;">{</span><span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.8</span><span style="color: #0000FF;">}},</span>
end function
<span style="color: #0000FF;">{{</span><span style="color: #000000;">19.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.0</span><span style="color: #0000FF;">},</span>
 
<span style="color: #0000FF;">{</span><span style="color: #000000;">28.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">26.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">25.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">22.1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">27.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">13.7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23.9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21.6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20.4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">24.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">13.2</span><span style="color: #0000FF;">}},</span>
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},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">30.02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.99</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">30.11</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.97</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">30.01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.99</span><span style="color: #0000FF;">},</span>
{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}},
<span style="color: #0000FF;">{</span><span style="color: #000000;">29.89</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.93</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.72</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.98</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">30.02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29.98</span><span style="color: #0000FF;">}},</span>
{{17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">3.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.1</span><span style="color: #0000FF;">},</span>
{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}},
<span style="color: #0000FF;">{</span><span style="color: #000000;">490.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">340.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">433.9</span><span style="color: #0000FF;">}},</span>
{{19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">0.010268</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.000167</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.000167</span><span style="color: #0000FF;">},</span>
{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}},
<span style="color: #0000FF;">{</span><span style="color: #000000;">0.159258</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.136278</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.122389</span><span style="color: #0000FF;">}},</span>
{{30.02, 29.99, 30.11, 29.97, 30.01, 29.99},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">1.0</span><span style="color: #0000FF;">/</span><span style="color: #000000;">15</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10.0</span><span style="color: #0000FF;">/</span><span style="color: #000000;">62.0</span><span style="color: #0000FF;">},</span>
{29.89, 29.93, 29.72, 29.98, 30.02, 29.98}},
<span style="color: #0000FF;">{</span><span style="color: #000000;">1.0</span><span style="color: #0000FF;">/</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">/</span><span style="color: #000000;">50.0</span><span style="color: #0000FF;">}},</span>
{{3.0, 4.0, 1.0, 2.1},
<span style="color: #0000FF;">{{</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">23.0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">21</span><span style="color: #0000FF;">/</span><span style="color: #000000;">45.0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">/</span><span style="color: #000000;">38.0</span><span style="color: #0000FF;">},</span>
{490.2, 340.0, 433.9}},
<span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">/</span><span style="color: #000000;">44.0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">42</span><span style="color: #0000FF;">/</span><span style="color: #000000;">94.0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">/</span><span style="color: #000000;">22.0</span><span style="color: #0000FF;">}}},</span>
{{0.010268,0.000167,0.000167},
<span style="color: #000000;">correct</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0.021378001462867</span><span style="color: #0000FF;">,</span>
{0.159258,0.136278,0.122389}},
<span style="color: #000000;">0.148841696605327</span><span style="color: #0000FF;">,</span>
{{1.0/15,10.0/62.0},
<span style="color: #000000;">0.0359722710297968</span><span style="color: #0000FF;">,</span>
{1.0/10,2/50.0}},
<span style="color: #000000;">0.090773324285671</span><span style="color: #0000FF;">,</span>
{{9/23.0,21/45.0,0/38.0},
<span style="color: #000000;">0.0107515611497845</span><span style="color: #0000FF;">,</span>
{0/44.0,42/94.0,0/22.0}}},
<span style="color: #000000;">0.00339907162713746</span><span style="color: #0000FF;">,</span>
correct = {0.021378001462867,
<span style="color: #000000;">0.14884169660532752726574965384</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.545266866977794</span><span style="color: #0000FF;">}</span>
0.0359722710297968,
0.090773324285671,
<span style="color: #004080;">atom</span> <span style="color: #000000;">cerr</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
0.0107515611497845,
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
0.00339907162713746,
<span style="color: #004080;">atom</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">welch_ttest</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
0.52726574965384,
<span style="color: #0000FF;">?</span><span style="color: #000000;">r</span>
0.545266866977794}
<span style="color: #000000;">cerr</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">-</span><span style="color: #000000;">correct</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
atom cerr = 0
<span style="color: #0000FF;">?{</span><span style="color: #008000;">"cumulative error"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cerr</span><span style="color: #0000FF;">}</span>
for i=1 to length(tests) do
<!--</syntaxhighlight>-->
atom r = welch_ttest(tests[i])
?r
cerr += abs(r-correct[i])
end for
?{"cumulative error",cerr}</lang>
{{out}}
<pre>
Line 1,321 ⟶ 1,775:
0.5272657497
0.545266867
{"cumulative error",1.989380882e-14} -- (32 bit/p2js)
{"cumulative error",4.915115776e-15} -- (64-bit)
</pre>
Line 1,328 ⟶ 1,782:
 
=== Using NumPy & SciPy ===
<langsyntaxhighlight lang="python">import numpy as np
import scipy as sp
import scipy.stats
Line 1,345 ⟶ 1,799:
 
welch_ttest(np.array([3.0, 4.0, 1.0, 2.1]), np.array([490.2, 340.0, 433.9]))
(-9.559497721932658, 2.0008523488562844, 0.01075156114978449)</langsyntaxhighlight>
=== Using betain from AS 63 ===
First, the implementation of betain (translated from the Stata program in the discussion page). The original Fortran code is under copyrighted by the Royal Statistical Society. The C translation is under GPL, written by John Burkardt. The exact statement of the RSS license is unclear.
 
<langsyntaxhighlight lang="python">import math
 
def betain(x, p, q):
Line 1,399 ⟶ 1,853:
else:
temp = psq
psq += 1</langsyntaxhighlight>
 
The Python code is then straightforward:
 
<langsyntaxhighlight lang="python">import math
 
def welch_ttest(a1, a2):
Line 1,421 ⟶ 1,875:
p = betain(df / (t**2 + df), df / 2, 1 / 2)
return t, df, p</langsyntaxhighlight>
 
'''Example'''
 
<langsyntaxhighlight lang="python">a1 = [3, 4, 1, 2.1]
a2 = [490.2, 340, 433.9]
print(welch_ttest(a1, a2))</langsyntaxhighlight>
 
'''Output'''
Line 1,433 ⟶ 1,887:
 
=={{header|R}}==
<langsyntaxhighlight Rlang="r">#!/usr/bin/R
 
printf <- function(...) cat(sprintf(...))
Line 1,470 ⟶ 1,924:
results <- t.test(z1,z2, alternative="two.sided", var.equal=FALSE)
printf("%.15g\n", results$p.value);
</syntaxhighlight>
</lang>
 
The output here is used to compare against C's output above.
Line 1,487 ⟶ 1,941:
{{trans|C}}
 
<langsyntaxhighlight lang="racket">#lang racket
(require math/statistics math/special-functions)
 
Line 1,536 ⟶ 1,990:
(p-value (list 3.0 4.0 1.0 2.1)
(list 490.2 340.0 433.9))))</langsyntaxhighlight>
 
{{out}}
Line 1,550 ⟶ 2,004:
Perhaps "inspired by C example" may be more accurate. Gamma subroutine from [[Gamma_function#Raku|Gamma function task]].
 
<syntaxhighlight lang="raku" perl6line>sub Γ(\z) {
constant g = 9;
z < .5 ?? π / sin(π × z) / Γ(1 - z) !!
Line 1,615 ⟶ 2,069:
[<3.0 4.0 1.0 2.1>],
[<490.2 340.0 433.9>]
) -> @left, @right { say p-value @left, @right }</langsyntaxhighlight>
{{out}}
<pre>0.0213780014628669
Line 1,631 ⟶ 2,085:
This uses the Soper reduction formula to evaluate the integral, which converges much more quickly than Simpson's formula.
 
<syntaxhighlight lang="raku" perl6line>sub lgamma ( Num(Real) \n --> Num ){
use NativeCall;
sub lgamma (num64 --> num64) is native {}
Line 1,735 ⟶ 2,189:
$error += abs($p-value - shift @answers);
}
printf("cumulative error is %g\n", $error);</langsyntaxhighlight>
{{out}}
<pre>p-value = 0.021378001462867
Line 1,749 ⟶ 2,203:
=={{header|Ruby}}==
{{trans|Perl}}
<langsyntaxhighlight Rubylang="ruby">def calculate_p_value(array1, array2)
return 1.0 if array1.size <= 1
return 1.0 if array2.size <= 1
Line 1,874 ⟶ 2,328:
 
printf("the cumulative error is %g\n", error)
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 1,890 ⟶ 2,344:
=={{header|SAS}}==
{{trans|Stata}}
<syntaxhighlight lang="text">data tbl;
input value group @@;
cards;
Line 1,900 ⟶ 2,354:
class group;
var value;
run;</langsyntaxhighlight>
 
'''Output'''
Line 2,061 ⟶ 2,515:
Implementation in IML:
 
<langsyntaxhighlight lang="sas">proc iml;
use tbl;
read all var {value} into x where(group=1);
Line 2,074 ⟶ 2,528:
p = 2*probt(-abs(t), df);
print t df p;
quit;</langsyntaxhighlight>
 
'''Output'''
Line 2,081 ⟶ 2,535:
 
=={{header|Scala}}==
<langsyntaxhighlight Scalalang="scala">import org.apache.commons.math3.distribution.TDistribution
 
object WelchTTest extends App {
Line 2,114 ⟶ 2,568:
println(s"\nSuccessfully completed without errors. [total ${scala.compat.Platform.currentTime - executionStart} ms]")
 
}</langsyntaxhighlight>
 
=={{header|Scilab}}==
Line 2,121 ⟶ 2,575:
Scilab will print a warning because the number of degrees of freedom is not an integer. However, the underlying implementation makes use of the [http://www.netlib.org/random/ dcdflib] Fortran library, which happily accepts a noninteger df.
 
<syntaxhighlight lang="text">x = [3.0,4.0,1.0,2.1];
y = [490.2,340.0,433.9];
n1 = length(x);
Line 2,130 ⟶ 2,584:
df = (v1/n1+v2/n2)^2/(v1^2/(n1^2*(n1-1))+v2^2/(n2^2*(n2-1)));
[p, q] = cdft("PQ", -abs(t), df);
[t df 2*p]</langsyntaxhighlight>
 
'''Output'''
Line 2,140 ⟶ 2,594:
=={{header|Sidef}}==
{{trans|Raku}}
<langsyntaxhighlight lang="ruby">func p_value (A, B) {
[A.len, B.len].all { _ > 1 } || return 1
 
Line 2,192 ⟶ 2,646:
tests.each_slice(2, {|left, right|
say p_value(left, right)
})</langsyntaxhighlight>
{{out}}
<pre>
Line 2,207 ⟶ 2,661:
Notice that here we use the option '''unequal''' of the '''ttest''' command, and not '''welch''', so that Stata uses the Welch-Satterthwaite approximation.
 
<langsyntaxhighlight lang="stata">mat a=(3,4,1,2.1,490.2,340,433.9\1,1,1,1,2,2,2)'
clear
svmat double a
Line 2,234 ⟶ 2,688:
 
di r(p)
.01075156</langsyntaxhighlight>
 
The computation can easily be implemented in Mata. Here is how to compute the t statistic (t), the approximate degrees of freedom (df) and the p-value (p).
 
<langsyntaxhighlight lang="stata">st_view(a=., ., .)
x = select(a[., 1], a[., 2] :== 1)
y = select(a[., 1], a[., 2] :== 2)
Line 2,252 ⟶ 2,706:
+----------------------------------------------+
1 | -9.559497722 2.000852349 .0107515611 |
+----------------------------------------------+</langsyntaxhighlight>
 
=={{header|Tcl}}==
Line 2,262 ⟶ 2,716:
This is not particularly idiomatic Tcl, but perhaps illustrates some of the language's relationship with the Lisp family.
 
<langsyntaxhighlight Tcllang="tcl">#!/usr/bin/tclsh
 
package require math::statistics
Line 2,337 ⟶ 2,791:
puts [pValue $left $right]
}
</syntaxhighlight>
</lang>
 
{{out}}
Line 2,351 ⟶ 2,805:
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
<langsyntaxhighlight ecmascriptlang="wren">import "./math" for Math, Nums
import "./fmt" for Fmt
 
var welch = Fn.new { |a, b|
Line 2,389 ⟶ 2,843:
var g3 = Math.gamma(nu/2 + 0.5).log
var f = Fn.new { |r| r.pow(nu/2-1) / (1 - r).sqrt }
return simpson0.call(2000, nu/(t*t + nu), f) / Math.exp(g1 + g2 - g3).exp
}
 
Line 2,408 ⟶ 2,862:
Fmt.print("$0.6f", pValue.call(d5, d6))
Fmt.print("$0.6f", pValue.call(d7, d8))
Fmt.print("$0.6f", pValue.call(x, y))</langsyntaxhighlight>
 
{{out}}
Line 2,421 ⟶ 2,875:
=={{header|zkl}}==
{{trans|C}}
<langsyntaxhighlight lang="zkl">fcn calculate_Pvalue(array1,array2){
if (array1.len()<=1 or array2.len()<=1) return(1.0);
 
Line 2,471 ⟶ 2,925:
foreach x in (cof){ ser+=(x/(y+=1)); }
return((2.5066282746310005 * ser / x).log() - tmp);
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">testSets:=T(
T(T(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),
T(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)),
Line 2,484 ⟶ 2,938:
 
foreach x,y in (testSets)
{ println("Test set 1 p-value = %f".fmt(calculate_Pvalue(x,y))); }</langsyntaxhighlight>
{{out}}
<pre>
2,122

edits