where : ibrtses delphi
Delphi - linear approximation 2D
disclaimer
the source code of this page may not appear correctly in certain browsers
due to special characters. Have a look at the source of this HTML page
with notepad instead
also known as least square linear regression, a line is fitted
through a bunch of points.
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
In the above version I omitted a procedure showpoints.
Without this ability, the points do not have to be stored.
Klaus Huemmeler found this optimized solution :
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;
Feedback is welcome
sponsored links
Delphi
home
last updated: 31.july.01
Copyright (99,2001) Ing.Büro R.Tschaggelar