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() - input[i].x ());
  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();
  26         Real y = input[i].y();
  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 }