Welch's t-test: Difference between revisions

m
→‎Using Burkardt's betain: while (1) changed to while(True), tab indent changed to 4 spaces, array2_size accounted for as well
(Scala solution added)
m (→‎Using Burkardt's betain: while (1) changed to while(True), tab indent changed to 4 spaces, array2_size accounted for as well)
Line 1,241:
 
def calculate_Pvalue (array1, array2):
ARRAY1_SIZE = len(array1)
if ARRAY1_SIZE <= 1:
ARRAY2_SIZE = len(array2)
return 1.0
if ARRAY1_SIZE <= 1:
ARRAY2_SIZE = len(array2)
return 1.0
elif ARRAY1_SIZE if ARRAY2_SIZE <= 1:
return 1.0
fmean1 = sum(array1)/ARRAY1_SIZE
fmean2 = sum(array2)/ARRAY2_SIZE
if fmean1 == fmean2:
return 1.0
unbiased_sample_variance1 = 0.0
unbiased_sample_variance2 = 0.0
for v in range(ARRAY1_SIZE):
unbiased_sample_variance1 += (array1[v] - fmean1) * (array1[v] - fmean1)
for v in range(ARRAY2_SIZE):
unbiased_sample_variance2 += (array2[v] - fmean2) * (array2[v] - fmean2)
unbiased_sample_variance1 = unbiased_sample_variance1/(ARRAY1_SIZE-1)
unbiased_sample_variance2 = unbiased_sample_variance2/(ARRAY2_SIZE-1)
WELCH_T_STATISTIC = (fmean1-fmean2)/math.sqrt(unbiased_sample_variance1/ARRAY1_SIZE+unbiased_sample_variance2/ARRAY2_SIZE)
DEGREES_OF_FREEDOM = pow((unbiased_sample_variance1/ARRAY1_SIZE+unbiased_sample_variance2/ARRAY2_SIZE),2.0) / ((unbiased_sample_variance1*unbiased_sample_variance1)/(ARRAY1_SIZE*ARRAY1_SIZE*(ARRAY1_SIZE-1))+ (unbiased_sample_variance2*unbiased_sample_variance2)/(ARRAY2_SIZE*ARRAY2_SIZE*(ARRAY2_SIZE-1))
)
)
a = DEGREES_OF_FREEDOM/2
value = DEGREES_OF_FREEDOM/(WELCH_T_STATISTIC*WELCH_T_STATISTIC+DEGREES_OF_FREEDOM)
beta = math.lgamma(a) + 0.57236494292470009-math.lgamma(a+0.5)
acu = 0.1E-14
if a <= 0.0:
return value
if value < 0.0 or 1.0 < value:
return value
if value == 0 or value == 1.0:
return value
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
term = 1.0
ai = 1.0
value = 1.0
ns = int(qq + cx*psq)
rx = xx/cx
temp = qq - ai
if ns == 0:
rx = xx
while (1True):
term = term * temp * rx / ( pp + ai )
value = value + term
temp = abs ( term )
if ((temp <= acu) and (temp <= acu * value)):
value = value * math.exp ( pp * math.log ( xx ) + ( qq - 1.0 ) * math.log ( cx ) - beta ) / pp
if indx:
value = 1.0 - value
break
 
ai = ai + 1.0
ns = ns - 1
 
if 0 <= ns:
temp = qq - ai
if ns == 0:
rx = xx
else:
temp = psq
psq = psq + 1.0
return value
#end of function
 
Line 1,377:
 
print 'the cumulative error is %g' % error</lang>
 
{{out}}
<pre>