Anonymous user
Welch's t-test: Difference between revisions
added Ruby
m (→{{header|Perl 6}}: added header for other section) |
(added Ruby) |
||
Line 1,829:
{{out}}
<pre>(0.021378001462867013 0.14884169660532798 0.035972271029796624 0.09077332428567102 0.01075139991904718)</pre>
=={{header|Ruby}}==
{{trans|Perl}}
<lang Ruby>
def calculate_Pvalue (array1, array2)
return 1.0 if array1.size <= 1
return 1.0 if array2.size <= 1
mean1 = array1.sum / array1.size
mean2 = array2.sum / array2.size
return 1.0 if mean1 == mean2
variance1 = 0.0
variance2 = 0.0
array1.each do |x|
variance1 += (mean1 - x)**2
end
array2.each do |x|
variance2 += (mean2 - x)**2
end
return 1.0 if variance1 == 0.0 && variance2 == 0.0
variance1 = variance1 / (array1.size - 1)
variance2 = variance2 / (array2.size - 1)
welch_t_statistic = (mean1 - mean2) / Math.sqrt(variance1 / array1.size + variance2/array2.size)
degrees_of_freedom = ((variance1/array1.size+variance2/array2.size)**2) / (
(variance1*variance1)/(array1.size*array1.size*(array1.size-1))+
(variance2*variance2)/(array2.size*array2.size*(array2.size-1)));
a = degrees_of_freedom / 2
value = degrees_of_freedom / (welch_t_statistic**2 + degrees_of_freedom)
beta = Math.lgamma(a)[0] + 0.57236494292470009 - Math.lgamma(a + 0.5)[0]
acu = 10**(-15)
return value if a <= 0
return value if value < 0.0 || 1.0 < value
return value if value == 0 or value == 1.0
psq = a + 0.5
cx = 1.0 - value
if a < psq * value
xx = cx
cx = value
pp = 0.5
qq = a
indx = 1
else
xx = value
pp = a
qq = 0.5
indx = 0
end
term = 1.0;
ai = 1.0;
value = 1.0;
ns = (qq + cx * psq).to_i;
#Soper reduction formula
rx = xx / cx
temp = qq - ai
while 1
term = term * temp * rx / ( pp + ai );
value = value + term;
temp = term.abs;
if temp <= acu && temp <= acu * value
#p "pp = #{pp}"
#p "xx = #{xx}"
#p "qq = #{qq}"
#p "cx = #{cx}"
#p "beta = #{beta}"
value = value * Math.exp(pp * Math.log(xx) + (qq - 1.0) * Math.log(cx) - beta) / pp;
#p "value = #{value}"
value = 1.0 - value if indx
if indx == 0
value = 1.0 - value
end
break
end
ai = ai + 1.0
ns = ns - 1
#p "ai = #{ai}"
#p "ns = #{ns}"
if 0 <= ns
temp = qq - ai
rx = xx if ns == 0
else
temp = psq
psq = psq + 1.0
end
end
return value
end
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]
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]
d3 = [17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8]
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]
d5 = [19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0]
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]
d7 = [30.02,29.99,30.11,29.97,30.01,29.99]
d8 = [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]
s1 = [1.0/15,10.0/62.0]
s2 = [1.0/10,2/50.0]
v1 = [0.010268,0.000167,0.000167]
v2 = [0.159258,0.136278,0.122389]
z1 = [9/23.0,21/45.0,0/38.0]
z2 = [0/44.0,42/94.0,0/22.0]
CORRECT_ANSWERS = [0.021378001462867, 0.148841696605327, 0.0359722710297968,
0.090773324285671, 0.0107515611497845, 0.00339907162713746, 0.52726574965384, 0.545266866977794]
pvalue = calculate_Pvalue(d1,d2)
error = (pvalue - CORRECT_ANSWERS[0]).abs
printf("Test sets 1 p-value = %.14g\n",pvalue)
pvalue = calculate_Pvalue(d3,d4)
error += (pvalue - CORRECT_ANSWERS[1]).abs
printf("Test sets 2 p-value = %.14g\n",pvalue)
pvalue = calculate_Pvalue(d5,d6)
error += (pvalue - CORRECT_ANSWERS[2]).abs
printf("Test sets 3 p-value = %.14g\n",pvalue)
pvalue = calculate_Pvalue(d7,d8)
error += (pvalue - CORRECT_ANSWERS[3]).abs
printf("Test sets 4 p-value = %.14g\n",pvalue)
pvalue = calculate_Pvalue(x,y)
error += (pvalue - CORRECT_ANSWERS[4]).abs
printf("Test sets 5 p-value = %.14g\n",pvalue)
pvalue = calculate_Pvalue(v1,v2)
error += (pvalue - CORRECT_ANSWERS[5]).abs
printf("Test sets 6 p-value = %.14g\n",pvalue)
pvalue = calculate_Pvalue(s1,s2)
error += (pvalue - CORRECT_ANSWERS[6]).abs
printf("Test sets 7 p-value = %.14g\n",pvalue)
pvalue = calculate_Pvalue(z1,z2)
error += (pvalue - CORRECT_ANSWERS[7]).abs
printf("Test sets z p-value = %.14g\n",pvalue)
printf("the cumulative error is %g\n", error)
</lang>
{{out}}
<pre>
Test sets 1 p-value = 0.021378001462867
Test sets 2 p-value = 0.14884169660533
Test sets 3 p-value = 0.035972271029797
Test sets 4 p-value = 0.090773324285671
Test sets 5 p-value = 0.010751561149784
Test sets 6 p-value = 0.0033990716271375
Test sets 7 p-value = 0.52726574965384
Test sets z p-value = 0.54526686697779
the cumulative error is 1.34961e-15
</pre>
=={{header|SAS}}==
| |||