/***************************************************************************
 *   Copyright (C) 2004 by Johan Maes                                      *
 *   on4qz@telenet.be                                                      *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

 /*
   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
*/
bool 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 (fabs(sxx) <= 0.001) return(FALSE);

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

   /* Compute the regression coefficient */
   if (fabs(syy) <= 0.00001)  *r = 1;
   else *r = sxy / sqrt(sxx * syy);
   return(TRUE);
}