Delphi - linear approximation 2D

theory

be the line equation : f(x):= {y:=a*x+b}
the square error for the ith point is : Sqr(f(xi)-yi) = Sqr(a*xi+b-yi)
for all points : sum(i:=1 to N, Sqr(f(xi)-yi) = sum(i:=1 to N, Sqr(a*xi+b-yi)) =
                 sum(i:=1 to N, a^2*xi^2+b^2+yi^2+2a*xi*b-2a*xi*yi-2b*yi)
as the variables are a and b, derive by a and b :

d/da sum(i:=1 to N, Sqr(a*xi+b-yi))= 
d/da sum(i:=1 to N, a^2*xi^2+b^2+yi^2+2a*xi*b-2a*xi*yi-2b*yi) =
 sum(i:=1 to N, 2a*xi^2 + 2xi*b - 2xi*yi)
 

d/db sum(i:=1 to N, Sqr(a*xi+b-yi))= 
d/db sum(i:=1 to N, a^2*xi^2+b^2+yi^2+2a*xi*b-2a*xi*yi-2b*yi) =
 sum(i:=1 to N, 2b +2a*xi - 2yi)

For the error to be minmal, set the derivatives to zero :

 sum(i:=1 to N, 2a*xi^2 + 2xi*b - 2xi*yi) = 0
 sum(i:=1 to N, 2b +2a*xi - 2yi) =0

 2a*sum(xi^2) + 2b*sum(xi) -2sum(xi*yi) = 0
 2Nb + 2a*sum(xi) - 2sum(yi) =0

this gives the equation system for a and b:

 | N      sum(xi)  |*|b| =| sum(yi)    |
 |sum(xi) sum(xi^2)| |a|  | sum(xi*yi) |

And some care has to be taken when the line is vertical :
when sum(xi)=0, the x and y are swapped.

Tom Berry inspired me to add this theory

the code

type
  Point2D= class
  private
   fx,fy:extended;
  public
   constructor create(x,y:extended);
   destructor destroy;   override;
  published
   property x:extended read fx write fx;
   property y:extended read fy write fy;
  end;


  TLinApprox2D = class(TObject)
  private
   pl:TList;
  public
    constructor create;
    destructor destroy;    override;
    procedure add(x,y:extended);
    procedure approx(var slope, offset:extended; var swap:boolean);
  end;


constructor Point2D.create(x,y:extended);
begin 
 inherited create; 
 fx:=x; fy:=y; 
end;

destructor Point2D.destroy;
begin 
 inherited destroy; 
end;

constructor TLinApprox2D.create;
begin
 inherited create;
 pl:=TList.create;
end;

destructor TLinApprox2D.destroy;
var i:integer;
begin
 for i:=0 to pl.count-1 do Point2D(pl.items[i]).destroy;
 pl.destroy;
 inherited destroy;
end;

procedure TLinApprox2D.add(x,y:extended);
var p:Point2D;
begin
 p:=Point2D.create(x,y);
 pl.add(p);
end;

procedure TLinApprox2D.approx(var slope, offset:extended; var swap:boolean);
var sx,sxx,sy,sxy,syy,a11,a12,a21,a22,c1,c2 :extended;
 i,n:integer;
begin
 n:=pl.count;
 sx:=0;
 for i:=0 to n-1 do sx:=sx+Point2D(pl.items[i]).x;
 sxx:=0;
 for i:=0 to n-1 do sxx:=sxx+sqr(Point2D(pl.items[i]).x);
 sy:=0;
 for i:=0 to n-1 do sy:=sy+Point2D(pl.items[i]).y;
 sxy:=0;
 for i:=0 to n-1 do sxy:=sxy+Point2D(pl.items[i]).y*Point2D(pl.items[i]).x;
 if sx=0 then begin
   syy:=0;
   for i:=0 to n-1 do syy:=syy+sqr(Point2D(pl.items[i]).y);
   a11:=n; a12:=sy; a21:=sy; a22:=syy; c1:=sx; c2:=sxy; swap:=true;
  end
 else begin
   a11:=n; a12:=sx; a21:=sx; a22:=sxx; c1:=sy; c2:=sxy; swap:=false;
 end;
 // solve equation system
 a12:=a12/a11; c1:=c1/a11; //a11:=1;
 a22:=a22/a21; c2:=c2/a21; //a21:=1;
 a22:=a22-a12; c2:=c2-c1;  //a21:=0;
 c2:=c2/a22; a22:=1;
 slope:=a22;
 c1:=c1-a12*c2; //a12:=0;
 offset:=c1;
end;

notes

 the swap in the proc approx is set true when the line is vertical, 
 meaning the coordinates are swapped.
 Use :
 MyApprox:=TLinApprox2D.create;
 MyApprox.add(0,0);
 MyApprox.add(1,0.5);
 ..
 MyApprox.approx(a,b,swap);
 
 The line will be y:=ax+b  unless swap, then x:=ay+b

storage- and runtimeoptimized version

TLinApprox = class(TObject)
private
   sx, sxx, sy, sxy : extended;
   n : integer;
public
   constructor create;
   procedure add(x,y:extended);
   function approx(var slope, offset:extended) : boolean;
end;

constructor TLinApprox.create;
begin
   inherited create;
   sx := 0; sxx := 0; sy := 0; sxy := 0; n := 0;
end;

procedure TLinApprox.add(x,y:extended);
begin
  sx  := sx  + x;
  sxx := sxx + sqr(x);
  sy  := sy  + y;
  sxy := sxy + (x*y);
  inc(n);
end;

function TLinApprox.approx(var slope, offset:extended) : boolean;
begin
  // solve equation system
  Result := false;
  if (n = 0) or (sx = 0) then EXIT;
  try
    Slope := (n*sxy-sx*sy)/(n*sxx-sqr(sx));
    offset:=(sxx*sy-sx*sxy)/(n*sxx-sqr(sx));
    Result := true;
  except end;
end;