计算t-test 的C程序

/* gdb output   程序还未调试成功：http://ubuntuforums.org/archive/index.php/t-412096.html   */

/*(gdb) run

Starting program: /home/nrh1/code/testt


Program received signal SIGSEGV, Segmentation fault.

0x0804967f in var ()   */

/* function to calculate ttest ttest.c
*/

#include <stdio.h>

#include <math.h>

#define square(x) (x*x)
void ttest(double *data1, unsigned int n1, double *data2, unsigned int n2, double *tstatistic, double *tprob);
int main()


{

/* declare arrays of x and y data from statistical methods in cell biology 2nd edition pg 45*/
double a[8]= {3.2,1.6,5.7,2.8,5.5,1.2,6.1,2.9};
double b[8]= {3.8,1.0,8.4,3.6,5.0,3.5,7.3,4.8};

unsigned int na=8;

unsigned int nb=8;
double* t;
double *prob;


ttest( a, na, b, nb, &t, &prob);

printf("tvalue %f\n", &t);


printf("tvalue %f\n", &prob);

} /* program closing brace    */

/*uses part of stats lib stats.c the relevant functions are shown below ftest code is nt mine but works very well I've used it in other contexts   */

/*

ttest assumes homogenity of variance, uses ftest probability code to calculate probability of the significance of the t statistic

*/


double var(double *data, unsigned int n)

{

unsigned int i;
double average=0.0;
double s, temp, var;

for (i=0; i<n; n++)

average=average+data[n]; /* calculate average   */


average=average/n;


var=temp=s=0.0; /* now variance  */
for (i=0; i<n; n++)

{

s=data[n]-average;

temp=temp+s;

var=var+(s*s);


}
return (var-temp*temp/n)/(n-1);

} /* function closing brace   */
/* end of function */

void ttest(double *data1, unsigned int n1, double *data2, unsigned int n2, double *tstatistic, double *tprob)

{

double nmean(unsigned int count, double *data); /* function declarations within function    */
double pof(double F, unsigned int df1, unsigned int df2);
double var(double *data, unsigned int n);
double svar=0.0;
double var1=0.0;
double var2=0.0;
double mean1=0.0;
double mean2=0.0;
double df=0.0;
/*double *tstatistic=0.0;  */
/*double *tprob=0.0     */


mean1=nmean(n1, data1); /*calculate means for individual data sets    */

mean2=nmean(n2, data2);

var1=var(data1, n1); /*calculate variances for individual data sets     */

var2=var(data2, n2);

df=n1+n2-2; /* calculate overall df                 */

svar=((n1-1)*var1+(n2-1)*var2)/df; /* pooled variance  */
*tstatistic=(mean1-mean2)/sqrt(svar*(1.0/n1+1.0/n2)); /* calculate t value    */
*tprob=pof(((*tstatistic)*(*tstatistic)), 1, df); /*probability        */

} /* function closing brace      */

/* end of function */


/*FUNCTION pof: probability of F */
/*ALGORITHM Compute probability of F ratio.

*/

/*

example if F=1.432923 and df1=3 and df2=5 p=0.337564
*/
double pof(double F, unsigned int df1, unsigned int df2)

{

int i, j;
int a, b;
double w, y, z, d, p;

if (F < F_EPSILON || df1 < 1 || df2 < 1)

p=1.0;
else

{

a = df1%2 ? 1 : 2;

b = df2%2 ? 1 : 2;

w = (F * df1) / df2;


z = 1.0 / (1.0 + w);
if (a == 1)
if (b == 1)

{

p = sqrt (w);

y = I_PI; /* 1 / 3.14159 */

d = y * z / p;

p = 2.0 * y * atan (p);

}
else

{

p = sqrt (w * z);

d = 0.5 * p * z / w;

}
else if (b == 1)

{

p = sqrt (z);

d = 0.5 * z * p;

p = 1.0 - p;

}
else

{

d = z * z;

p = w * z;

}

y = 2.0 * w / z;

for (j = b + 2; j <= df2; j += 2)

{

d *= (1.0 + a / (j - 2.0)) * z;

p = (a == 1 ? p + d * y / (j - 1.0) : (p + w) * z);

}


y = w * z;

z = 2.0 / z;

b = df2 - 2;
for (i = a + 2; i <= df1; i += 2)

{

j = i + b;

d *= y * j / (i - 2.0);

p -= z * d / j;

}
/* correction for approximation errors suggested in certification */
if (p < 0.0)

p = 0.0;
else if (p > 1.0)

p = 1.0;

p= (1.0-p);
return p;

}
return p;

}