Anonymous user
Welch's t-test: Difference between revisions
added lgamma C code in case user's computer does not have lgamma
(answered critiques of Rdm and Sonia) |
(added lgamma C code in case user's computer does not have lgamma) |
||
Line 163:
}
</lang>
{{out}}
<pre>Test sets 1 p-value = 0.021378
Line 169 ⟶ 170:
Test sets 4 p-value = 0.090773
Test sets 5 p-value = 0.010751</pre>
'''If''' your computer does not have <code>lgamma</code>, use the following program:
<lang C>
#include <stdio.h>
#include <math.h>
double LnGamma(const double xx) {
unsigned int j;
double x,y,tmp,ser;
const double cof[6] = {
76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2,-0.5395239384953e-5
};
y = x = xx;
tmp = x + 5.5 - (x + 0.5) * logl(x + 5.5);
ser = 1.000000000190015;
for (j=0;j<=5;j++)
ser += (cof[j] / ++y);
return(log(2.5066282746310005 * ser / x) - tmp);
}
double calculate_Pvalue (const double *array1, const size_t array1_size, const double *array2, const size_t array2_size) {
if (array1_size <= 1) {
return 1.0;
}
if (array2_size <= 1) {
return 1.0;
}
double mean1 = 0.0, mean2 = 0.0;
for (size_t x = 0; x < array1_size; x++) {
mean1 += array1[x];
}
for (size_t x = 0; x < array2_size; x++) {
mean2 += array2[x];
}
if (mean1 == mean2) {
return 1.0;
}
mean1 /= array1_size;
mean2 /= array2_size;
double variance1 = 0.0, variance2 = 0.0;
for (size_t x = 0; x < array1_size; x++) {
variance1 += (array1[x]-mean1)*(array1[x]-mean1);
}
for (size_t x = 0; x < array2_size; x++) {
variance2 += (array2[x]-mean2)*(array2[x]-mean2);
}
if ((variance1 == 0.0) && (variance2 == 0.0)) {
return 1.0;
}
variance1 = variance1/(array1_size-1);
variance2 = variance2/(array2_size-1);
const double WELCH_T_STATISTIC = (mean1-mean2)/sqrt(variance1/array1_size+variance2/array2_size);
const double DEGREES_OF_FREEDOM = pow((variance1/array1_size+variance2/array2_size),2.0)
/
(
(variance1*variance1)/(array1_size*array1_size*(array1_size-1))+
(variance2*variance2)/(array2_size*array2_size*(array2_size-1))
);
const double a = DEGREES_OF_FREEDOM/2, x = DEGREES_OF_FREEDOM/(WELCH_T_STATISTIC*WELCH_T_STATISTIC+DEGREES_OF_FREEDOM);
const unsigned short int N = 65535;
const double h = x/N;
double sum1 = 0.0, sum2 = 0.0;
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)));
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(LnGamma(a)+0.57236494292470009-LnGamma(a+0.5)));
if ((isfinite(return_value) == 0) || (return_value > 1.0)) {
return 1.0;
} else {
return return_value;
}
}
//-------------------
int main(void) {
const double 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};
const double 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};
const double d3[] = {17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8};
const double 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};
const double d5[] = {19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0};
const double 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};
const double d7[] = {30.02,29.99,30.11,29.97,30.01,29.99};
const double d8[] = {29.89,29.93,29.72,29.98,30.02,29.98};
const double x[] = {3.0,4.0,1.0,2.1};
const double y[] = {490.2,340.0,433.9};
printf("Test sets 1 p-value = %lf\n",calculate_Pvalue(d1,sizeof(d1)/sizeof(*d1),d2,sizeof(d2)/sizeof(*d2)));
printf("Test sets 2 p-value = %lf\n",calculate_Pvalue(d3,sizeof(d3)/sizeof(*d3),d4,sizeof(d4)/sizeof(*d4)));
printf("Test sets 3 p-value = %lf\n",calculate_Pvalue(d5,sizeof(d5)/sizeof(*d5),d6,sizeof(d6)/sizeof(*d6)));
printf("Test sets 4 p-value = %lf\n",calculate_Pvalue(d7,sizeof(d7)/sizeof(*d7),d8,sizeof(d8)/sizeof(*d8)));
printf("Test sets 5 p-value = %lf\n",calculate_Pvalue(x,sizeof(x)/sizeof(*x),y,sizeof(y)/sizeof(*y)));
return 0;
}
</lang>
{{out}}
<pre>Test sets 1 p-value = 0.021378
Test sets 2 p-value = 0.148842
Test sets 3 p-value = 0.035972
Test sets 4 p-value = 0.090773
Test sets 5 p-value = 0.010752</pre>
=={{header|Go}}==
| |||