I've got this one little function that calculates standard deviation. All the separate peices work, but together it blows up. I simply cannot figure out why.

I'm trying to get this formula (its hard to show in ascii): Sx = sqrt( (#SUM OF FOR I TO N# Xi^2 - n(average of X)^2) / N)

Here's the function. FYI I put in a zillion extra ()'s to make sure it did things right...

See, the output with DEBUG set is as follows:Code:double standev( void ) { //standard deviation int i; double top=0; for (i=0;i<N;i++) { top += ((x[i]*x[i]) - (((double)(N))*(avgx*avgx))); if (DEBUG) printf("[+] top += %f %f*%f=%f - %f * %f * %f = %f\n",(x[i]*x[i]) - (N*(avgx*avgx)),x[i],x[i],(x[i]*x[i]),(double)(N),avgx,avgx,(avgx*avgx)); } if (DEBUG) printf("[+] standard dev: %f\n[+] top: %f, divide by N : %f\n",sqrt(top/(double)(N)),top,(double)(N)); return sqrt(top/(double)(N)); }

See, all the various peices hold the right values, but the final is crazy way off.Code:[+] top += -3252.785233 3.000000*3.000000=9.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3245.785233 4.000000*4.000000=16.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3236.785233 5.000000*5.000000=25.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3225.785233 6.000000*6.000000=36.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3212.785233 7.000000*7.000000=49.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3197.785233 8.000000*8.000000=64.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3180.785233 9.000000*9.000000=81.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3161.785233 10.000000*10.000000=100.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3140.785233 11.000000*11.000000=121.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3117.785233 12.000000*12.000000=144.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3092.785233 13.000000*13.000000=169.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3065.785233 14.000000*14.000000=196.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3036.785233 15.000000*15.000000=225.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -3005.785233 16.000000*16.000000=256.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -2972.785233 17.000000*17.000000=289.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -2937.785233 18.000000*18.000000=324.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -2900.785233 19.000000*19.000000=361.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -2861.785233 20.000000*20.000000=400.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -2820.785233 21.000000*21.000000=441.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262 [+] top += -2777.785233 22.000000*22.000000=484.000000 - 20.000000 * 12.770641 * 12.770641 = 163.089262

Here's the full source:

Thanks a lot, I'm pulling my hair out.Code:#include <stdio.h> #include <stdlib.h> #include <math.h> #define N 20 #define DEBUG 1 double x[N],y[N],avgx,avgy,avgxy,r; void getvals( void ) { //read from file... int i; for (i=0;i<20;i++) { x[i]=(double)(i) + 3.0; y[i]=(double)(i); if (DEBUG) printf("[+] (%f,%f)\n",x[i],y[i]); } } void getavg( void ) { int i; double sumx, sumy, sumxy; //before here, validate so that sizeof(x)=sizeof(y) for (i=0;i<N;i++) { sumx += x[i]; sumy += y[i]; sumxy += x[i]*y[i]; } avgx = sumx/(double)(N); avgy = sumy/(double)(N); avgxy = sumxy/(double)(N); if (DEBUG) printf("[+] avgx %f avgy %f avgxy %f\n",avgx,avgy,avgxy); } double standev( void ) { //standard deviation int i; double top=0; for (i=0;i<N;i++) { top += ((x[i]*x[i]) - (((double)(N))*(avgx*avgx))); if (DEBUG) printf("[+] top += %f %f*%f=%f - %f * %f * %f = %f\n",(x[i]*x[i]) - (N*(avgx*avgx)),x[i],x[i],(x[i]*x[i]),(double)(N),avgx,avgx,(avgx*avgx)); } if (DEBUG) printf("[+] standard dev: %f\n[+] top: %f, divide by N : %f\n",sqrt(top/(double)(N)),top,(double)(N)); return sqrt(top/(double)(N)); } void findlinear( void ) { double slope,yint,sx; //here we need to find the slope, then change the //y-int so that the line runs through (avgx,avgy) getvals(); getavg(); sx = standev(); slope = (avgxy - (avgx*avgy))/(sx*sx); yint = avgy - (slope*avgx); printf("y = %fx + %f\n"); } int main(int argc, char *argv[]) { findlinear(); return 0; }