lily/leastsquares.cc

   1 #include "leastsquares.hh"
   2 #include "warn.hh"
   3
   4 void
   5 Least_squares::OK() const
   6 {
   7   assert (input.size() > 1);
   8   Real dx = 0.0;
   9   for (int i=1; i < input.size(); i++)
  10     dx += abs (input[i-1][X_AXIS] - input[i][X_AXIS]);
  11
  12   assert (dx);
  13 }
  14
  15 void
  16 Least_squares::minimise (Real &coef, Real &offset)
  17 {
  18   Real sx = 0.0;
  19   Real sy = 0.0;
  20   Real sqx =0.0;
  21   Real sxy = 0.0;
  22
  23   for (int i=0; i < input.size();i++)
  24     {
  25       Real x=input[i][X_AXIS];
  26       Real y = input[i][Y_AXIS];
  27       sx += x;
  28       sy += y;
  29       sqx += sqr (x);
  30       sxy += x*y;
  31     }
  32   int N = input.size();
  33
  34   coef =0.0;
  35   offset =0.;
  36
  37   Real den = (N*sqx - sqr (sx));
  38   if (!N || !den)
  39     programming_error ("Least_squares::minimise():  Nothing to minimise");
  40   coef = (N * sxy - sx*sy)/den;
  41   offset = (sy - coef * sx)/N;
  42 }