Welch's t-test: Difference between revisions
Content added Content deleted
(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: | Line 1,241: | ||
def calculate_Pvalue (array1, array2): |
def calculate_Pvalue (array1, array2): |
||
ARRAY1_SIZE = len(array1) |
|||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
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 (True): |
|||
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 |
#end of function |
||
Line 1,377: | Line 1,377: | ||
print 'the cumulative error is %g' % error</lang> |
print 'the cumulative error is %g' % error</lang> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |