Regression / Abstand zu Punkten

[Furtbichler]

Polynomiale Regression ist doch bekannt, womit das Polynom f(x)=a*x^2+b*x + c definiert ist, ergo die Kurve. Das simple Näherungsverfahren mit Regula Falsi kann man nun auf die Gerade beschränken. Ist weniger Arbeit. Man kann das auch analytisch lösen, aber das ist mir ne halbe Nummer zu hoch.

[cltom]

vielen Dank für die Antworten. Um die Fragestellung etwas zu präzisieren: es geht nicht um die Findung der Parabel, polynomiale Regression ist bekannt, da gibt es ja genug Material dazu. Gesucht ist aber eben jene Gerade, die die Parabel (eigentlich die ursprünlichen Datenpunkte) in der beschriebenen Weise schneidet (also wo Schnittflächen über und unter der Parabel gleich sind). Da ist schon die Frage, ob man sich mit dem Polynom-Fit einen Gefallen tut. Weil ich dort ja eine gewisse Abweichung erzeuge und dann später, beim Finden der Gerade mit gefitteten Werten arbeite.

Da wäre wohl der Weg von sx2008 denkbar, die Parameter per Zufall oft genug zu variieren. Denkbar wäre wohl auch, die Gerade in kleinen Winkelschritten zu rotieren und ganz einfach die Differenzen zu den ursprünglichen Datenpunkten zu variieren.

Eine analytische Lösung wär natürlich interessant, mir aber wohl eher zwei Nummern zu hoch ...

dank und gruß
cltom

[Furtbichler]

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.

[cltom]

hoi, super, danke für den Ansatz, das schaut gut aus.

Was ich nicht beachtet/erwähnt hatte, was das Problem etwas vereinfacht: es ist ein Punkt der Geraden bekannt, nämlich ein bestimmter Punkt auf der Parabel (der sich aus einer anderen Bedingung ergibt), dh. man braucht im Grunde nur die Steigung variieren.

Aber Dein Verfahren ist mal allgemein gültig, was nützlich ist, weil die Bedingung mit dem gegebenen Punkt mitunter nicht immer existiert.

[Furtbichler]

Wenn nur die Steigung anzupassen ist, dann reicht ein Durchgang der Methode ':IterateA()'.

[jfheins]

Ich würde hier mal den analytischen Ansatz vorschlagen.

Gegeben sei eine Parabel y=a*x^2+bx+c und ein Punkt (px, py) der auf der Parabel liegt. Die gerade soll nun durch den Punkt gehen, und ein Integral soll zu Null werden. Dazu müssen noch die Intervallgrenzen definiert werden - ich nehme hier mal 0-2 an.

Dann sieht das einfach so aus: int(a*x^2+bx+c)(dx, 0, 2) = int(m*(x-xp)+yp)(dx, 0, 2)
Also flux die Integrale gelöst:

(8a)/3+2 (b+c) = 2*(m*(1-px)+py)

und schon hat man eine Gleichung die man einfach nach m auflösen kann. m ist die gesuchte Steigung

[Furtbichler]

Das Integral soll nicht null werden, sondern so klein wie möglich, oder? Ansonsten: Weiter so!