Welch's t-test: Difference between revisions
Content added Content deleted
No edit summary |
|||
Line 45: | Line 45: | ||
<lang C>#include <stdio.h> |
<lang C>#include <stdio.h> |
||
#include <math.h> |
#include <math.h> |
||
double calculate_Pvalue (const double *array1, const size_t array1_size, const double *array2, const size_t array2_size) { |
double calculate_Pvalue (const double *array1, const size_t array1_size, const double *array2, const size_t array2_size) { |
||
if (array1_size <= 1) { |
if (array1_size <= 1) { |
||
Line 85: | Line 84: | ||
); |
); |
||
const double a = DEGREES_OF_FREEDOM/2, x = DEGREES_OF_FREEDOM/(WELCH_T_STATISTIC*WELCH_T_STATISTIC+DEGREES_OF_FREEDOM); |
const double a = DEGREES_OF_FREEDOM/2, x = DEGREES_OF_FREEDOM/(WELCH_T_STATISTIC*WELCH_T_STATISTIC+DEGREES_OF_FREEDOM); |
||
const |
const unsigned short int N = 65535; |
||
const double h = x/N; |
|||
double sum1 = 0.0, sum2 = 0.0; |
double sum1 = 0.0, sum2 = 0.0; |
||
for(unsigned int i = 0;i < N; i++) { |
for(unsigned short int i = 0;i < N; i++) { |
||
sum1 += (pow(h * i + h / 2.0,a-1))/(sqrt(1-(h * i + h / 2.0))); |
sum1 += (pow(h * i + h / 2.0,a-1))/(sqrt(1-(h * i + h / 2.0))); |
||
⚫ | |||
for(unsigned int i = 1;i < N;i++) { |
|||
sum2 += (pow(h * i,a-1))/(sqrt(1-h * i)); |
sum2 += (pow(h * i,a-1))/(sqrt(1-h * i)); |
||
} |
} |
||
double return_value = ((h / 6.0) * ((pow(x,a-1))/(sqrt(1-x)) + 4.0 * sum1 + 2.0 * sum2))/(expl(lgammal(a)+0.57236494292470009-lgammal(a+0.5)));//0.5723... is lgammal(0.5) |
|||
if ((isinf(return_value) != 0) || (return_value > 1.0)) { |
|||
return 1.0; |
|||
} else { |
|||
return return_value; |
|||
⚫ | |||
} |
} |
||
//------------------- |
//------------------- |