/***************************************** * Module for basic statistical Functions * *****************************************/ #include #include #include double Dot_Product(int n, double *x, double *y) { int i; double temp; temp=0.0; for (i=0; i assumed.\n"); variance = (sum_sq-Sqr(sum_)/(1.0*n))/(1.0*n); } if (variance < 0.0) { // an error exists printf(" From NormalStats: negative variance %f\n", variance); *std_dev = -1.0; } else *std_dev = sqrt(variance); *avg = sum_/n; } //NormalStats /*--------------------------------------------------------- ! For data to be represented by y=ax+b, calculates linear ! regression coefficients, sample standard error of y on x, ! and sample correlation coefficients. Sets r=0 if an error ! exists. If the intercept coefficient a is set to 0 on ! input, the regression is forced through (0,0). !--------------------------------------------------------*/ void LinearReg(double *x, double *y, int n, char flag, double *a, double *b, double *s_yx, double *r) { double avg, std_dev; double sum_x,sum_y,sum_xy,sum_xx,sum_yy,temp; sum_x = Sum(n,x); sum_y = Sum(n,y); sum_xy = Dot_Product(n,x,y); sum_xx = Dot_Product(n,x,x); sum_yy = Dot_Product(n,y,y); if (*a != 0.0) { // calculate full expression temp = n*sum_xx - Sqr(sum_x); *a = (sum_y*sum_xx - sum_x*sum_xy)/temp; *b = (n*sum_xy - sum_x*sum_y)/temp; *s_yx = sqrt((sum_yy - (*a)*sum_y - (*b)*sum_xy)/n); } else { //just calculate slope *b = sum_y/sum_x; *s_yx = sqrt((sum_yy - 2.0*(*b)*sum_xy + (*b)*(*b)*sum_xx)/n); } if (flag=='s' || flag=='S') *s_yx = *s_yx * sqrt((1.0*n)/(1.0*(n-2))); // Use NormalStats to get standard deviation of y NormalStats(y,n,flag,&avg,&std_dev); if (std_dev > 0.0) { temp = 1.0 - Sqr(*s_yx/std_dev); if (temp > 0.0) *r = sqrt(temp); else { // an error exists *r = 0.0; printf(" From LinearReg: error in temp %f\n", temp); } } else // an error exists *r = 0.0; } // LinearReg /*---------------------------------------------------------- ! Generates an array of normal random numbers from pairs of ! uniform random numbers in range [0,1]. !---------------------------------------------------------*/ void NormalArray(double *a, int n) { int i; double pi, temp, u1,u2; double MAX_VALUE = 1.0; pi = 3.1415926535; // fills array with uniform random for (i=0; i