Welch's t-test: Difference between revisions

{{trans|Python}}

I p <= 0 | q <= 0 | x < 0 | x > 1

X ValueError(0)

V a1 = [Float(3), 4, 1, 2.1]

V a2 = [Float(490.2), 340, 433.9]

{{out}}

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).

#include <math.h>

#include <stdlib.h>

return 0;

}

{{out}}

'''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:

#include <math.h>

}

=={{header|Fortran}}==

Alternatively, the program shows the p-value computed using the IMSL '''BETAI''' function.

use tdf_int

implicit none

print *, t, df, p

print *, betai(df / (t**2 + df), 0.5d0 * df, 0.5d0)

'''Output'''

With Absoft Pro Fortran, compile with <code>af90 -m64 pvalue.f90 %SLATEC_LINK%</code>.

implicit none

integer :: n1, n2

call welch_ttest(4, x, 3, y, t, df, p)

print *, t, df, p

'''Output'''

Instead of implementing the t-distribution by ourselves, we bind to GNU Scientific Library:

use, intrinsic :: iso_c_binding

print *, t, df, p

'''Output'''

=={{header|Go}}==

import (

func(r float64) float64 { return math.Pow(r, ν/2-1) / math.Sqrt(1-r) }) /

math.Exp(g1+g2-g3)

{{out}}

<pre>

Implementation:

'a b steps'=. 3{.y,128

size=. (b - a)%steps

hi=. v%(t^2)+v

(F f. simpson integrate lo,hi) % 0.5 B v%2

<code>integrate</code> and <code>simpson</code> are from the [[Numerical_integration#J|Numerical integration]] task.

Data for task examples:

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

d8=: 29.89 29.93 29.72 29.98 30.02 29.98

d9=: 3 4 1 2.1

Task examples:

0.021378

d3 p2_tail d4

0.0907733

d9 p2_tail da

=={{header|Java}}==

Using the '''[http://commons.apache.org/proper/commons-math/ Apache Commons Mathematics Library]'''.

public class WelchTTest {

System.out.println("p = " + res[2]);

}

'''Result'''

jq's `lgamma` returns the natural logarithm of the absolute value of the gamma function of x.

# Sample variance using division by (length-1)

pValue(d5; d6),

pValue(d7; d8),

{{out}}

<pre>

{{works with|Julia|0.6}}

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]

ttest = UnequalVarianceTTest(y1, y2)

println("\nData:\n y1 = $y1\n y2 = $y2\nP-value for unequal variance TTest: ", round(pvalue(ttest), 4))

{{out}}

=={{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:

typealias Func = (Double) -> Double

println(f.format(p2Tail(da7, da8)))

println(f.format(p2Tail(x, y)))

{{out}}

=={{header|Nim}}==

{{trans|Kotlin}}

func sqr(f: float): float = f * f

echo p2Tail(Da5, Da6).formatFloat(ffDecimal, 6)

echo p2Tail(Da7, Da8).formatFloat(ffDecimal, 6)

{{out}}

=={{header|Maple}}==

uses Statistics;

local n1:=nops(x),n2:=nops(y),

y:=[490.2,340,433.9]:

WelschTTest(x,y);

=={{header|Octave}}==

{{trans|Stata}}

y = [490.2,340.0,433.9];

n1 = length(x);

ans =

=={{header|PARI/GP}}==

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);

{{out}}

<pre>%1 = 0.010751561149784496723954539777213062928</pre>

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

use List::Util qw(sum);

use Math::AnyNum qw(gamma pi);

my ($left, $right) = splice(@tests, 0, 2);

print p_value($left, $right), "\n";

{{out}}

<pre>

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

use warnings;

use List::Util 'sum';

printf("p-value = %.14g\n",$pvalue);

}

{{out}}

<pre>p-value = 0.021378001462867

{{trans|Go}}

{{trans|Kotlin}}

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

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

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

{{out}}

<pre>

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)

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #000080;font-style:italic;">--&lt;copy of gamma from Gamma_function#Phix&gt;</span>

<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>

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

{{out}}

<pre>

=== Using NumPy & SciPy ===

import scipy as sp

import scipy.stats

welch_ttest(np.array([3.0, 4.0, 1.0, 2.1]), np.array([490.2, 340.0, 433.9]))

=== 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.

def betain(x, p, q):

else:

temp = psq

The Python code is then straightforward:

def welch_ttest(a1, a2):

p = betain(df / (t**2 + df), df / 2, 1 / 2)

'''Example'''

a2 = [490.2, 340, 433.9]

'''Output'''

=={{header|R}}==

printf <- function(...) cat(sprintf(...))

results <- t.test(z1,z2, alternative="two.sided", var.equal=FALSE)

printf("%.15g\n", results$p.value);

The output here is used to compare against C's output above.

{{trans|C}}

(require math/statistics math/special-functions)

(p-value (list 3.0 4.0 1.0 2.1)

{{out}}

Perhaps "inspired by C example" may be more accurate. Gamma subroutine from [[Gamma_function#Raku|Gamma function task]].

constant g = 9;

z < .5 ?? π / sin(π × z) / Γ(1 - z) !!

[<3.0 4.0 1.0 2.1>],

[<490.2 340.0 433.9>]

{{out}}

<pre>0.0213780014628669

This uses the Soper reduction formula to evaluate the integral, which converges much more quickly than Simpson's formula.

use NativeCall;

sub lgamma (num64 --> num64) is native {}

$error += abs($p-value - shift @answers);

}

{{out}}

<pre>p-value = 0.021378001462867

=={{header|Ruby}}==

{{trans|Perl}}

return 1.0 if array1.size <= 1

return 1.0 if array2.size <= 1

printf("the cumulative error is %g\n", error)

{{out}}

<pre>

=={{header|SAS}}==

{{trans|Stata}}

input value group @@;

cards;

class group;

var value;

'''Output'''

Implementation in IML:

use tbl;

read all var {value} into x where(group=1);

p = 2*probt(-abs(t), df);

print t df p;

'''Output'''

=={{header|Scala}}==

object WelchTTest extends App {

println(s"\nSuccessfully completed without errors. [total ${scala.compat.Platform.currentTime - executionStart} ms]")

=={{header|Scilab}}==

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.

y = [490.2,340.0,433.9];

n1 = length(x);

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);

'''Output'''

=={{header|Sidef}}==

{{trans|Raku}}

[A.len, B.len].all { _ > 1 } || return 1

tests.each_slice(2, {|left, right|

say p_value(left, right)

{{out}}

<pre>

Notice that here we use the option '''unequal''' of the '''ttest''' command, and not '''welch''', so that Stata uses the Welch-Satterthwaite approximation.

clear

svmat double a

di r(p)

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).

x = select(a[., 1], a[., 2] :== 1)

y = select(a[., 1], a[., 2] :== 2)

+----------------------------------------------+

1 | -9.559497722 2.000852349 .0107515611 |

=={{header|Tcl}}==

This is not particularly idiomatic Tcl, but perhaps illustrates some of the language's relationship with the Lisp family.

package require math::statistics

puts [pValue $left $right]

}

{{out}}

{{libheader|Wren-math}}

{{libheader|Wren-fmt}}

import "/fmt" for Fmt

Fmt.print("$0.6f", pValue.call(d5, d6))

Fmt.print("$0.6f", pValue.call(d7, d8))

{{out}}

=={{header|zkl}}==

{{trans|C}}

if (array1.len()<=1 or array2.len()<=1) return(1.0);

foreach x in (cof){ ser+=(x/(y+=1)); }

return((2.5066282746310005 * ser / x).log() - tmp);

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

foreach x,y in (testSets)

{{out}}

<pre>