/*
   Linear Regression
   y(x) = a + b x, for n samples
   The following assumes the standard deviations are unknown for x and y
   Return a, b and r the regression coefficient
*/
int LinRegress(double *x,double *y,int n,double *a,double *b,double *r)
{
   int i;
   double sumx=0,sumy=0,sumx2=0,sumy2=0,sumxy=0;
   double sxx,syy,sxy;

   *a = 0;
   *b = 0;
   *r = 0;
   if (n < 2)
      return(FALSE);

   /* Conpute some things we need */
   for (i=0;i<n;i++) {
      sumx += x[i];
      sumy += y[i];
      sumx2 += (x[i] * x[i]);
      sumy2 += (y[i] * y[i]);
      sumxy += (x[i] * y[i]);
   }
   sxx = sumx2 - sumx * sumx / n;
   syy = sumy2 - sumy * sumy / n;
   sxy = sumxy - sumx * sumy / n;

   /* Infinite slope (b), non existant intercept (a) */
   if (ABS(sxx) == 0)
      return(FALSE);

   /* Calculate the slope (b) and intercept (a) */
   *b = sxy / sxx;
   *a = sumy / n - (*b) * sumx / n;

   /* Compute the regression coefficient */
   if (ABS(syy) == 0)
      *r = 1;
   else
      *r = sxy / sqrt(sxx * syy);

   return(TRUE);
}