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

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 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;

var p:Point2D;
begin
p:=Point2D.create(x,y);
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.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;
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;

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