Welch's t-test: Difference between revisions

{{trans|C}}

<lang zkl>fcn calculate_Pvalue(array1,array2){

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

mean1,mean2 := array1.sum(0.0),array2.sum(0.0);

if(mean1==mean2) return(1.0);

mean1/=array1.len();

mean2/=array2.len();

variance1:=array1.reduce('wrap(sum,x){ sum + (x-mean1).pow(2) },0.0);

variance2:=array2.reduce('wrap(sum,x){ sum + (x-mean2).pow(2) },0.0);

variance1/=(array1.len() - 1);

variance2/=(array2.len() - 1);

WELCH_T_STATISTIC:=(mean1-mean2)/

(variance1/array1.len() + variance2/array2.len()).sqrt();

DEGREES_OF_FREEDOM:=

( variance1/array1.len() + variance2/array2.len() ).pow(2) // numerator

/ (

(variance1*variance1)/(array1.len().pow(2)*(array1.len() - 1)) +

(variance2*variance2)/(array2.len().pow(2)*(array2.len() - 1))

);

a:=DEGREES_OF_FREEDOM/2;

x:=DEGREES_OF_FREEDOM/( WELCH_T_STATISTIC.pow(2) + DEGREES_OF_FREEDOM );

N,h := 65535, x/N;

sum1,sum2 := 0.0, 0.0;

foreach i in (N){

sum1+=((h*i + h/2.0).pow(a - 1))/(1.0 - (h*i + h/2.0)).sqrt();

sum2+=((h*i).pow(a - 1))/(1.0 - h*i).sqrt();

}

return_value:=((h/6.0)*( x.pow(a - 1)/(1.0 - x).sqrt() +

4.0*sum1 + 2.0*sum2) ) /

((0.0).e.pow(lngammal(a) + 0.57236494292470009 - lngammal(a + 0.5)));

if(return_value > 1.0) return(1.0); // or return_value is infinite, throws

return_value;

}

fcn lngammal(xx){

var [const] cof=List( // static

76.18009172947146, -86.50532032941677,

24.01409824083091, -1.231739572450155,

0.1208650973866179e-2,-0.5395239384953e-5

);

y:=x:=xx;

tmp:=x + 5.5 - (x + 0.5) * (x + 5.5).log();

ser:=1.000000000190015;

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

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

}</lang>

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

T(T(17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8),

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

T(T(19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0),

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

T(T(30.02,29.99,30.11,29.97,30.01,29.99),

T(29.89,29.93,29.72,29.98,30.02,29.98)),

T(T(3.0,4.0,1.0,2.1),T(490.2,340.0,433.9)) );

foreach x,y in (testSets)

{ println("Test set 1 p-value = %f".fmt(calculate_Pvalue(x,y))); }</lang>

Test set 1 p-value = 0.021378

Test set 1 p-value = 0.148842

Test set 1 p-value = 0.035972

Test set 1 p-value = 0.090773

Test set 1 p-value = 0.010752