计算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;
}