---

lsutil.c

Go to the documentation of this file.

/*
  lsutil.c

  Least squares linear fit code

  MODIFICATIONS

  28-Jun-2000 WPS added __attribute__((unused)) to ident.

  */

/* ident tag */
#ifndef __GNUC__
#ifndef __attribute__
#define __attribute__(x)
#endif
#endif

static char __attribute__((unused)) ident[] = "$Id: lsutil.c,v 1.3 2000/10/27 20:34:42 terrylr Exp $";

static double sum(double *x, int n)
{
  double retval = 0.0;

  n--;
  while (n >= 0)
    {
      retval += x[n];
      n--;
    }

  return retval;
}

static double sumSq(double *x, int n)
{
  double retval = 0.0;

  n--;
  while (n >= 0)
    {
      retval += x[n]*x[n];
      n--;
    }

  return retval;
}

static double sumProd(double *x, double *y, int n)
{
  double retval = 0.0;

  n--;
  while (n >= 0)
    {
      retval += x[n]*y[n];
      n--;
    }

  return retval;
}

/* returns the sum of (x[i] - d) y[i]*/
static double sumProdOD(double *x, double *y, int n, double d)
{
  double retval = 0.0;

  n--;
  while (n >= 0)
    {
      retval += (x[n]-d)*y[n];
      n--;
    }

  return retval;
}

static double sumSqOD(double *x, int n, double d)
{

  double retval = 0.0;

  n--;
  while (n >= 0)
    {
      retval += (x[n]-d)*(x[n]-d);
      n--;
    }

  return retval;
}

/*
  given arrays of points x, y, and number of points n,
  computes the least squares best fit line y = a x + b,
  putting the coefficients in the passed pointer args.

  Returns 0 if OK, -1 if only one point, or a = inf
  */
int lsCoeff(double *x, double *y, int n, double *a, double *b)
{
  double sumX;
  double sumY;
  double sumXY;
  double sumXSq;
  double denom;

  if (n <= 1)
    {
      return -1;
    }

  sumX = sum(x, n);
  sumY = sum(y, n);
  sumXY = sumProd(x, y, n);
  sumXSq = sumSq(x, n);

  denom = n*sumXSq - sumX*sumX;
  if (denom == 0.0)
    {
      return -1;
    }
  *a = (n*sumXY - sumX*sumY)/denom;
  *b = (sumY - (*a)*sumX)/n;

  return 0;
}

/*
  given arrays of points x, y, number of points n, and the
  point d where y must be 0, computes the least squares best
  fit line y = a (x - d), putting the coefficient in the passed
  pointer arg a.

  Returns 0 if OK, -1 if only one point, or a = inf
  */
int lsCoeffOD(double *x, double *y, int n, double d, double *a)
{
  double sumXODY;
  double sumXODSq;

  if (n <= 1)
    {
      return -1;
    }

  sumXODY = sumProdOD(x, y, n, d);
  sumXODSq = sumSqOD(x, n, d);

  if (sumXODSq == 0.0)
    {
      return -1;
    }

  *a = sumXODY/sumXODSq;

  return 0;
}

---