source: tags/release-0.9/Utils/linear_regression.cc @ 813

#include <iostream>
#include <math.h>
void linear_regression(int num, float *x, float *y, int ilow, int ihigh, float &slope, float &errSlope, float &intercept, float &errIntercept, float &r)
{
  /**
   *   linear_regression()
   *
   * Computes the linear best fit to data - y = a*x + b, where x and y are 
   * arrays of size num, a is the slope and b the y-intercept.
   * The values used in the arrays are those from ilow to ihigh. 
   * (ie. if the full arrays are being used, then ilow=0 and ihigh=num-1.
   * Returns the values of slope & intercept (with errors) 
   * as well as r, the regression coefficient. 
   */

  if (ilow>ihigh) {
    std::cerr << "Error! linear_regression.cc :: ilow (" << ilow 
              << ") > ihigh (" << ihigh << ")!!\n";
    std::exit(1);
  }
  if (ihigh>num-1) {
    std::cerr << "Error! linear_regression.cc :: ihigh (" <<ihigh
              << ") out of bounds of array (>" << num-1 << ")!!\n";
    std::exit(1);
  }
  if(ilow<0){
    std::cerr << "Error! linear_regression.cc :: ilow (" << ilow
              << ") < 0. !!\n";
    std::exit(1);
  }

  double sumx,sumy,sumxx,sumxy,sumyy;
  sumx=0.;
  sumy=0.;
  sumxx=0.;
  sumxy=0.;
  sumyy=0.;
  int count=0;
  for (int i=ilow;i<=ihigh;i++){
    count++;
    sumx = sumx + x[i];
    sumy = sumy + y[i];
    sumxx = sumxx + x[i]*x[i];
    sumxy = sumxy + x[i]*y[i];
    sumyy = sumyy + y[i]*y[i];
  }
  if(count!=(ihigh-ilow+1)){
    std::cerr << "Whoops!! count (" << count <<") doesn't equal what it should ("
              << (ihigh-ilow+1) << ")!!\n";
    std::exit(1);
  }

  slope = (count*sumxy - sumx*sumy)/(count*sumxx - sumx*sumx);
  errSlope = count / (count*sumxx - sumx*sumx);
  intercept = (sumy*sumxx - sumxy*sumx)/(count*sumxx - sumx*sumx);
  errIntercept = sumxx / (count*sumxx - sumx*sumx);

  r = (count*sumxy - sumx*sumy)/(sqrt(count*sumxx-sumx*sumx)*sqrt(count*sumyy-sumy*sumy));
}