/* Program for mean, standard deviation */
/* Program for least-squares fit to a straight line */
/* Program for computing correlation coefficient "r" */
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);
}
