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 double N = 3000, h = x/N;
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));
}
}
return ((h / 6.0) * ((pow(x,a-1))/(sqrt(1-x)) + 4.0 * sum1 + 2.0 * sum2))/(expl(lgammal(a)+lgammal(0.5)-lgammal(a+0.5)));
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;
}
}
}
//-------------------
//-------------------