Sei [x1,y1]..[xn,yn] die Punktemenge, dann wäre (yn-y0)/(xn-x0) ein sinnvoller Startwert für die Steigung 'a' und (yn+y0)/2 ein sinnvoller Startwert für den Offset der Geraden.
Ich habe es so gelöst: In einer Schleife werden abwechselnd a und b so iteriert, das die resultierende Fläche minimiert wird. Dafür verwende ich ein stark vereinfachtes regula falsi, ein ziemlich lahmes Verfahren. Ich würde hier vermutlich mit ein wenig mehr Elan das Newtonsche Näherungsverfahren nehmen, weil es schneller ist. Newton benötigt zwar die 1.Ableitung, aber das geht schon, weil wir das mit [f(x+dx)-f(x)]/dx annähern können (x ist hier 'a' oder 'b').
Für die Fläche nehme ich einfach die Summe der einzelnen Vierecke, die durch die Punkte X_i+1/x_i und f(i+1)/f(i) aufgespannt wird. Hierbei ist f(i) = y_i - a*x_i+b die Differenz zwischen dem Kontrolpunkt und dem Punkt auf der gesuchten Gerade an dieser Stelle.
Hier mein zusammengerotzer Ansatz (der sogar vielleicht funktioniert).
Delphi-Quellcode:
Procedure Iterate();
Var
a1, a2, a, b: Double;
// calc area between control points p and a*x+b
Function _CalcArea(a, b: Double): Double;
Var
i: Integer;
dx, dy: Double;
Begin
Result := 0;
For i := 0
To NPoints - 1
Do Begin
dx := p[i + 1].x - p[i].x;
dy := p[i + 1].y - (a * p[i + 1].x + b) + p[i].y - (a * p[i].x + b);
result := result + abs(dx * dy / 2);
End
End;
Procedure _IterateA(
Var a: Double; b: double);
Var
da, area, area1: Double;
Begin
da := Max(0.1, a / 10);
area := _CalcArea(a, b);
Repeat
area1 := _CalcArea(a + da, b);
// area for new candidate
If area1 < area
Then Begin // any better?
a := a + da;
// yes
area := area1;
End
Else
da := -da / 2;
// no, reverse and lower delta
Until abs(da) < 1E-5;
End;
Procedure _IterateB(a: Double;
Var b: Double);
Var
db, area, area1: Double;
Begin
area := _CalcArea(a, b);
db := max(0.1, b / 10);
Repeat
area1 := _CalcArea(a, b +
db);
If area1 < area
Then Begin
b := b +
db;
area := area1;
End
Else
db := -
db / 2;
Until abs(
db) < 1E-5;
End;
Begin
a := (p[NPoints].Y - p[0].Y) / (p[NPoints].X - p[0].x);
b := (p[0].Y + p[NPoints].Y) / 2;
a2 := -1;
Repeat
a1 := a2;
_IterateA(a, b);
_IterateB(a, b);
a2 := _calcArea(a, b);
writeln(a: 8: 4, '
', b: 8: 4, '
', a2);
Until abs(a1 - a2) < 1E-5;
readln;
End;
So ganz sicher bin ich nicht, dass das immer funktioniert. So eine unabhängige Iteration über die beiden Parameter a und b ist normalerweise nur mit genetischer Programmierung/Iteration halbwegs sicher.
Ich würde als Startwerte ruhig die durch die lineare Regression vorgegebenen Werte verwenden und da/
db kleiner wählen, damit ist man dann auf der sicherereren Seite, denke ich.