I am new to C.
This function is used to calculate a basic F-statistic.
Most of the time it works fine, once in a while it will fail with
several values
set to -1.#IND00.
double calcF(int *marker, int *groups, double *groupmean, double
*pheno, int numStrains, int Largest, double totalmean)
{
int gpctr = 0; // group iterator
/******************** Count number in each group
***************************/
int m = 0;
int tmp = 0;
for(gpctr = 1; gpctr<=Largest; gpctr++)
{ tmp = 0;
for(m=0;m<numStrains;m++)
{
if( gpctr == marker[m] )
{
tmp++;
}
}
}
/***************** Calc mean for each group ***********************/
double tmpsum = 0.0;
for(gpctr = 1; gpctr<=Largest; gpctr++)
{
tmpsum = 0.0;
for(m=0;m<numStrains;m++)
{
if( gpctr == marker[m])
{
tmpsum += pheno[m];
}
groupmean[gpctr] = (tmpsum/groups[gpctr]);
}
}
/*****************************************************************/
double sstot = 0.0; //sum of squares total
double ssbet = 0.0; //sum of squares between
double ssw = 0.0; //sum of squares within
double dfbet = 0.0; // degrees of freedom between
double dfw = 0.0; // degrees of freedom within
double f = 0.0; // F-statistic
/******** CALC sum of square total *********/
double sstmp = 0.0;
for(m=0;m<numStrains;m++)
{
sstmp += ((pheno[m] - totalmean) * (pheno[m] - totalmean) );
}
sstot = sstmp;
/******** CALC sum of square between *********/
sstmp = 0.0;
for(gpctr = 1; gpctr<=Largest; gpctr++)
{
sstmp += ((groupmean[gpctr] - totalmean) *
(groupmean[gpctr] - totalmean));
}
ssbet = sstmp; //########## THIS SOMETIMES SET TO -1.#IND00
/******** CALC sum of square within *********/
if( ssbet == 0)
{
ssw = sstot;
}else{
ssw = sstot - ssbet;
}
/******** CALC degrees of freedom *********/
dfbet = (Largest - 1);
dfw = (numStrains - Largest);
/******** Calc F-statistic ************************/
if(ssw == 0 || dfbet == 0 || dfw == 0 || ssbet == 0)
{
return f = 1000.1;
}else{
f = (ssbet/dfbet)/(ssw/dfw);
return f;
}