/* Program for mean, standard deviation */
/* Program for least-squares fit to a straight line */
/* Program for computing correlation coefficient "r" */

#include <stdio.h>
#include <math.h>
#include <float.h>

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 */