Rosetta Code

Welch's t-test: Difference between revisions

Page Discussion

Welch's t-test (view source)

Revision as of 20:41, 30 January 2018

219 bytes removed ,  6 years ago
→‎{{header|C}}: merged betain function inside Pvalue function, should run faster and produce smaller executable
(→‎{{header|Perl}}: Add second version)
(→‎{{header|C}}: merged betain function inside Pvalue function, should run faster and produce smaller executable)
Line 72:
This program, for example, pvalue.c, can be compiled by
 
<code>clang -o pvalue pvalue.c -Wall -pedantic -std=c11 -lm -O2O3</code>
 
or
 
<code>gcc -o pvalue pvalue.c -Wall -pedantic -std=c11 -lm -O2O4</code>.
 
This shows how pvalue can be calculated from any two arrays, using Welch's 2-sided t-test, which doesn't assume equal variance.
Line 84:
#include <stdlib.h>
 
double betainPvalue ( const double x*restrict ARRAY1, const doublesize_t pARRAY1_SIZE, const double q*restrict ARRAY2, const doublesize_t beta,ARRAY2_SIZE) int{//calculate *restricta ifaultp-value )based on an array
if (ARRAY1_SIZE <= 1) {
return 1.0;
} else if (ARRAY2_SIZE <= 1) {
return 1.0;
}
double fmean1 = 0.0, fmean2 = 0.0;
for (size_t x = 0; x < ARRAY1_SIZE; x++) {//get sum of values in ARRAY1
if (isfinite(ARRAY1[x]) == 0) {//check to make sure this is a real numbere
puts("Got a non-finite number in 1st array, can't calculate P-value.");
exit(EXIT_FAILURE);
}
fmean1 += ARRAY1[x];
}
fmean1 /= ARRAY1_SIZE;
for (size_t x = 0; x < ARRAY2_SIZE; x++) {//get sum of values in ARRAY2
if (isfinite(ARRAY2[x]) == 0) {//check to make sure this is a real number
puts("Got a non-finite number in 2nd array, can't calculate P-value.");
exit(EXIT_FAILURE);
}
fmean2 += ARRAY2[x];
}
fmean2 /= ARRAY2_SIZE;
// printf("mean1 = %lf mean2 = %lf\n", fmean1, fmean2);
if (fmean1 == fmean2) {
return 1.0;//if the means are equal, the p-value is 1, leave the function
}
double unbiased_sample_variance1 = 0.0, unbiased_sample_variance2 = 0.0;
for (size_t x = 0; x < ARRAY1_SIZE; x++) {//1st part of added unbiased_sample_variance
unbiased_sample_variance1 += (ARRAY1[x]-fmean1)*(ARRAY1[x]-fmean1);
}
for (size_t x = 0; x < ARRAY2_SIZE; x++) {
unbiased_sample_variance2 += (ARRAY2[x]-fmean2)*(ARRAY2[x]-fmean2);
}
// printf("unbiased_sample_variance1 = %lf\tunbiased_sample_variance2 = %lf\n",unbiased_sample_variance1,unbiased_sample_variance2);//DEBUGGING
unbiased_sample_variance1 = unbiased_sample_variance1/(ARRAY1_SIZE-1);
unbiased_sample_variance2 = unbiased_sample_variance2/(ARRAY2_SIZE-1);
const double WELCH_T_STATISTIC = (fmean1-fmean2)/sqrt(unbiased_sample_variance1/ARRAY1_SIZE+unbiased_sample_variance2/ARRAY2_SIZE);
const double DEGREES_OF_FREEDOM = pow((unbiased_sample_variance1/ARRAY1_SIZE+unbiased_sample_variance2/ARRAY2_SIZE),2.0)//numerator
/
(
(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))
);
// printf("Welch = %lf DOF = %lf\n", 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);
if ((isinf(x) != 0) || (isnan(x) != 0)) {
return 1.0;
}
if ((isinf(a) != 0) || (isnan(a) != 0)) {
return 1.0;
}
 
/* Purpose:
/******************************************************************************/
/*
Purpose:
 
BETAIN computes the incomplete Beta function ratio.
Line 130 ⟶ 179:
Beta function ratio.
*/
const double beta = lgammal(a)+0.57236494292470009-lgammal(a+0.5);
{
const double acu = 0.1E-14;
double ai;
// double betain;
double cx;
int indx;
Line 147 ⟶ 195:
 
value = x;
// ifault = 0;
//Check the input arguments.
if ( (pa <= 0.0)) {// || (q0.5 <= 0.0 )){
// *ifault = 1;
// return value;
}
if ( x < 0.0 || 1.0 < x )
{
// *ifault = 2;
return value;
}
Line 164 ⟶ 212:
return value;
}
psq = pa + q0.5;
cx = 1.0 - x;
 
if ( pa < psq * x )
{
xx = cx;
cx = x;
pp = q0.5;
qq = pa;
indx = 1;
}
Line 178 ⟶ 226:
{
xx = x;
pp = pa;
qq = q0.5;
indx = 0;
}
Line 234 ⟶ 282:
 
return value;
}
 
double Pvalue (const double *restrict ARRAY1, const size_t ARRAY1_SIZE, const double *restrict ARRAY2, const size_t ARRAY2_SIZE) {//calculate a p-value based on an array
if (ARRAY1_SIZE <= 1) {
return 1.0;
} else if (ARRAY2_SIZE <= 1) {
return 1.0;
}
double fmean1 = 0.0, fmean2 = 0.0;
for (size_t x = 0; x < ARRAY1_SIZE; x++) {//get sum of values in ARRAY1
if (isfinite(ARRAY1[x]) == 0) {//check to make sure this is a real numbere
puts("Got a non-finite number in 1st array, can't calculate P-value.");
exit(EXIT_FAILURE);
}
fmean1 += ARRAY1[x];
}
fmean1 /= ARRAY1_SIZE;
for (size_t x = 0; x < ARRAY2_SIZE; x++) {//get sum of values in ARRAY2
if (isfinite(ARRAY2[x]) == 0) {//check to make sure this is a real number
puts("Got a non-finite number in 2nd array, can't calculate P-value.");
exit(EXIT_FAILURE);
}
fmean2 += ARRAY2[x];
}
fmean2 /= ARRAY2_SIZE;
// printf("mean1 = %lf mean2 = %lf\n", fmean1, fmean2);
if (fmean1 == fmean2) {
return 1.0;//if the means are equal, the p-value is 1, leave the function
}
double unbiased_sample_variance1 = 0.0, unbiased_sample_variance2 = 0.0;
for (size_t x = 0; x < ARRAY1_SIZE; x++) {//1st part of added unbiased_sample_variance
unbiased_sample_variance1 += (ARRAY1[x]-fmean1)*(ARRAY1[x]-fmean1);
}
for (size_t x = 0; x < ARRAY2_SIZE; x++) {
unbiased_sample_variance2 += (ARRAY2[x]-fmean2)*(ARRAY2[x]-fmean2);
}
// printf("unbiased_sample_variance1 = %lf\tunbiased_sample_variance2 = %lf\n",unbiased_sample_variance1,unbiased_sample_variance2);//DEBUGGING
unbiased_sample_variance1 = unbiased_sample_variance1/(ARRAY1_SIZE-1);
unbiased_sample_variance2 = unbiased_sample_variance2/(ARRAY2_SIZE-1);
const double WELCH_T_STATISTIC = (fmean1-fmean2)/sqrt(unbiased_sample_variance1/ARRAY1_SIZE+unbiased_sample_variance2/ARRAY2_SIZE);
const double DEGREES_OF_FREEDOM = pow((unbiased_sample_variance1/ARRAY1_SIZE+unbiased_sample_variance2/ARRAY2_SIZE),2.0)//numerator
/
(
(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))
);
// printf("Welch = %lf DOF = %lf\n", 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);
if ((isinf(x) != 0) || (isnan(x) != 0)) {
return 1.0;
}
if ((isinf(a) != 0) || (isnan(a) != 0)) {
return 1.0;
}
 
int *restrict z = 0;
return betain(x, a, 0.5, lgammal(a)+0.57236494292470009-lgammal(a+0.5), z);
}
//-------------------
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/258292"