/* Program for mean, standard deviation */ /* Program for least-squares fit to a straight line */ /* Program for computing correlation coefficient "r" */ #include #include #include int main() { /* Initialize all the parameters and type variables */ FILE *fp; /* file pointer */ static double x[] = { 1981., 1982., 1983., 1984., 1985., 1986., 1987., 1988., 1989., 1990., 1991. }; /* x values (years) */ static double y[] = { 90.08, 90.57, 90.76, 91.30, 91.57, 92.44, 92.87, 93.68, 93.85, 94.49, 94.88 }; /* y values (reliabilities) */ /* The variables sum_{something} are used as accumulators. */ /* a, b are the y-intercept and slope, respectively of the least-squares line. */ /* r is the correlation coefficient (the so-called Pearson's r). */ double Delta, a, b, r, sum_x, sum_x2, sum_y, sum_y2, sum_xy; /* mu_{something} is a mean, sigma_{something} is a standard deviation (unbiased for least-squares = n-2 degrees of freedom). */ double mu_x, mu_y, sigma_x, sigma_y, sigma_xy; /* The integer i is an index. (unsigned, short integer) */ int i; /* The long integer n is the count of the number of data points. */ long int n; /* Begin computations */ /* Initialize n, x[n], y[n], and give the output file a name. Then you are ready to begin execution. */ n = 11; /* Make sure the accumulators are all initialized to zero. */ sum_x = (double) 0.0; sum_y = (double) 0.0; sum_y2 = (double) 0.0; sum_x2 = (double) 0.0; sum_xy = (double) 0.0; sigma_x = (double) 0.0; sigma_y = (double) 0.0; sigma_xy = (double) 0.0; /* The so-called "for loop" computes the sum_{something}s */ for (i=0; i < n; ++i){ /* i marches from 0 to n-1. */ sum_y += (double) y[i]; /* All the arrays in "C" */ sum_y2 += (double) y[i]*y[i]; /* programs begin with index */ sum_x += (double) x[i]; /* zero. */ sum_x2 += (double) x[i]*x[i]; sum_xy += (double) x[i]*y[i]; } /* Display the results of the computation of all the sums, Delta, the slope, and the y-intercept for y = a + b * x, the least-squares straight line */ printf("\n sum_x = %.16lf", sum_x); printf("\n sum_y = %.16lf", sum_y); printf("\n sum_x2 = %.16lf", sum_x2); printf("\n sum_y2 = %.16lf", sum_y2); printf("\n sum_xy = %.16lf", sum_xy); Delta = (double) n*sum_x2 - sum_x*sum_x; printf("\n Delta = %.16lf", Delta); a = (double) (sum_x2*sum_y - sum_x*sum_xy)/Delta; b = (double) (n*sum_xy - sum_x*sum_y)/Delta; printf("\n A (y-intercept) = %.16lf", a); /* Use Capital "A" here */ printf("\n B (slope) = %.16lf", b); /* to be consistent with */ mu_y = (double) sum_y/n; /* the text book (Taylor). */ mu_x = (double) sum_x/n; /* Display (print to screen) the means only if needed. */ /* printf("\n mu_y = %.16lf", mu_y); */ /* printf("\n mu_x = %.16lf", mu_x); */ for (i=0; i < n; ++i){ /* Biased estimates for sigma_x and sigma_y */ sigma_y += (mu_y-y[i])*(mu_y-y[i]); /* We are just using the */ sigma_x += (mu_x-x[i])*(mu_x-x[i]); /* sigma_{something}s here */ sigma_xy += (mu_x-x[i])*(mu_y-y[i]); /* as accumulators. */ } /* compute the correlation coefficient from the values in the accumulators. This is not the actual formula. */ r = sigma_xy/sqrt(sigma_x*sigma_y); printf("\n correlation coefficient = %.16lf", r); sigma_y = 0.0; for (i=0; i < n; ++i){ sigma_y += (y[i] - a - b*x[i])*(y[i] - a - b*x[i]); /* This is the least-squares data to be plotted. */ printf("\n %d %.4lf %.4lf %.4lf", i, x[i], y[i],a + b*x[i]); } /* Unbiased estimates for sigma_x and sigma_y */ /* Note that there are n-2 degrees of freedom (not n-1). */ sigma_y = sqrt((double) sigma_y/(n-2.0)); printf("\n sigma_y (unbiased) = %.16lf", sigma_y); sigma_x = sigma_y * sqrt((double) sum_x2/Delta); printf("\n sigma_x (unbiased) = %.16lf", sigma_x); /* Output to file. Make sure you have a unique filename. */ fp = fopen("pgm3.txt","w"); /* Open the file */ fprintf(fp,"\n A (y-intercept) = %.16lf", a); fprintf(fp,"\n B (slope) = %.16lf", b); fprintf(fp,"\n correlation coefficient = %.16lf", r); fprintf(fp,"\n sigma_y (unbiased) = %.16lf", sigma_y); fprintf(fp,"\n x[i] y[i] a+b*x[i]"); for (i=0; i < n; ++i){ fprintf(fp,"\n %.4lf %.4lf %.4lf",x[i],y[i], a + b*x[i]); } fclose(fp); /* Close the file */ return(0); } /* End Of File */