Regression / Abstand zu Punkten

[cltom]

Hallo im neuen Jahr!

gegeben sind ein paar x-y-Daten (reell, kaum mehr als 50), die annähernd in Form einer Parabel verlaufen. Gesucht ist eine Gerade, die so durch die Punkte geht, dass die Fläche zwischen der Gerade und der (durch die Punkte gebildeten) Kurve möglichst klein ist. Anbei eine symbolische Skizze, die das klarer machen sollte. Effektiv gemessen werden diskrete Punkte. Durch diese wird dann eine Kurve gelegt (die rote Linie).

Gesucht also jene Gerade, bei der die Differenz der Summen der gelben und hellblauen Flächen möglichst klein ist. Ist das einfach die lineare Regression der Punkte? Es sieht im Grunde so aus, das ist aber wohl kein gültiger Beweis. Eine lineare Regression wäre natürlich etwas einfacher als eine Optimierung der Gerade auf möglichst geringe Differenzflächen. Am schönsten wäre es wohl, beides zu programmieren und dann zu schauen, wie groß die Unterschiede sind je nach Datensatz. Aber ich vermute, diese Fragestellung ist nicht ganz neu.

danke und gruß
cltom

[jfheins]

Nein, die beiden Kriterien sind nicht identisch.

Stelle dir einfach mal vor, deine Punkte sind links auf dem aufsteigenden Ast der Parabel dicht und rechts dünn verteilt. Die Parabel lässt sich wunderbar fitten. Benutzt du nun die Parabel um die Gerade zu fitten, geht die Information flöten wo viele Punkte sind und wo wenig. Der RMS fehler der Punkte zur Gerade kann dadurch sehr groß werden.

Idealerweise solltest du direkt die Parabel weiterverwenden (wenn du weißt, dass es sich um eine Parabel handelt) oder eben eine Gerade durch die Punkte fitten.

[ASM]

Korrekte Durchführung der Regression mittels Least Squares Fitting.
Alternativ wäre noch möglich das Simplex-Verfahren.

[sx2008]

Als erstes bräuchtest du mal eine Formel für die liegende Parabel.
Ich würde mal

Code:

f(x) = a * (x-b)^c + d

versuchen wobei a bis d die gesuchten Parameter sind. (startwerte: a=1, b=0, c=0.5, d=0)
Evtl. könnte man das Koordinatensystem so drehen dass man mit

Code:

f(x) = a*x*x + b*x + c

arbeiten kann.

Zu Beginn setzt man die Parameter auf einigermassen sinnvolle Werte und beginnt die Berechnungen.
Man berechnet für jeden Punkt die y-Abweichung zur Funktion und quadriert diese Abweichung.
Alle quadrierten Abweichungen werden aufsummiert.
Dann verändert man der Reihe nach die verschiedenen Parameter zufallsgesteuert zwischen 0 bis 10%
und rechnet erneut.
Hat sich die Summe der quadrierten Abweichungen verkleinert dann werden diese Parameter zum neuen Ausgangspunkt.
Nach vielleicht 10000 Iterationen sollte die Kurve nahezu ideal zu den Punkten passen.
Mit zunehmender Zahl der Iterationen muss dann auch die zufällige Abweichung verringert werden (z.B. 10% , 5%, 2.5%,...).

[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!

[jfheins]

Zitat:

Zitat von Furtbichler (Beitrag 1242027)

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

Zum Glück als Frage, von daher: Nein! :stupid:

Die Fläche des Integral ist ja vorzeichenbehaftet. Wenn ich das hier richtig interpretiere:

Zitat:

Gesucht also jene Gerade, bei der die Differenz der Summen der gelben und hellblauen Flächen möglichst klein ist.

Dann meinst er den Absolutbetrag der Differenz. Und der wird minimal gleich Null.

Ansonsten dürfte die Differenz gegen minus unendlich gehen - das ist aber keine zufriedenstellenden Lösung.

Oben habe ich jedoch noch einen klitzekleinen Denkfehler gemacht: Nicht die beiden Integrale sollen gleich werden, sondern das Integral der Differenz soll 0 werden!

für y=m*(x-d) + e und y=a*x^2+b*x+c ergibt sich somit: integrate(a*x^2 + (b-m) * x + c + m*d - e) (x, 0, 2) = 0
das ergibt: a*(8/3 - 2*d^2) - 2*(d-1)*(b-m) = 0 für d <> 1.

Alles natürlich ohne Gewähr :stupid:

[Furtbichler]

Aber die Fläche soll doch minimiert werden und sich nicht bezüglich 'oberhalb' und 'unterhalb' der Geraden aufheben, oder ist das das gleiche? :oops: Ach so...

[BUG]

Wenn die Flächen nur durch die Gerade und die Parabel begrenzt sind, sollten die gelben Flächen unendlich groß sein.
Damit macht die Frage nicht viel Sinn ... oder habe ich da irgendwas übersehen :gruebel:

Ich vermute mal, dass der gezeigte Ausschnitt (auf dem Bild) auch eine Begrenzung darstellt.
Dann stellt sich die Frage, ob es vielleicht viele solche Geraden gibt und ob da eine Bestimmte gesucht ist.
jfheins legt zusätzlich noch einen Punkt fest ... vielleicht hast du ja eine andere Beschränkung? Es fehlen also Details.

[Furtbichler]

Wir haben Stützwerte und ich würde annehmen, das es darum geht, die Fläche im Interval x1..xn zu minimieren. Das macht zumindest meine Frickeliteration.

[cltom]

vielen dank für die Antworten. Die analytischen Ansätze muss ich mir mal durchsehen. Als Zwischenbericht von meiner Seite: nach einigen Stunden mit realen Daten bin ich mittlerweile der Überzeugung, dass eine lineare Regression durch die Daten im betreffenden Parabelabschnitt ausreichen müsste. Die Daten streuen so sehr, dass die theoretische Forderung nach den gleich großen Integralen in der Praxis kaum Gewicht hat. Da muss ich mich mehr um die Frage kümmern, welche von den Messadaten ich überhaupt nehme ...

Zitat:

Zitat von BUG (Beitrag 1242326)

Wenn die Flächen nur durch die Gerade und die Parabel begrenzt sind, sollten die gelben Flächen unendlich groß sein.
Damit macht die Frage nicht viel Sinn ... oder habe ich da irgendwas übersehen :gruebel:

Ich vermute mal, dass der gezeigte Ausschnitt (auf dem Bild) auch eine Begrenzung darstellt.
Dann stellt sich die Frage, ob es vielleicht viele solche Geraden gibt und ob da eine Bestimmte gesucht ist.
jfheins legt zusätzlich noch einen Punkt fest ... vielleicht hast du ja eine andere Beschränkung? Es fehlen also Details.

hmm, hoffe eigentlich nicht, dass etwas fehlt. Warum sollten die gelben Flächen unendlich groß sein? Die Gerade liegt über oder unter der Parabel (je nach betrachtetem Intervall). Und in jedem Fall kann ich eine endliche Differenz zwischen einer Geradenfunktion und einer Parabelfunktion ausrechnen. Und auch das Integral ist in einem beschränkten Intervall endlich. Vielleicht war das offen. Es geht um einen diskreten Abschnitt der Parabel, nicht um den gesamten Verlauf von minus unendlich bis plus unendlich.

Zu ein paar anderen Fragen: die Bedingung, dass die gelben und blauen Flächen gleich sind, würde ich auch so auslegen, dass die Differenz der Integrale nicht Null ist, sondern nur möglichst klein.

[bernhard_LA]

unter
http://www.tu-ilmenau.de/num/team/we...publikationen/

gibt Algorithmen zu diesem Problem ( Kapitel 8 )

[jfheins]

Zitat:

Zitat von BUG (Beitrag 1242326)

jfheins legt zusätzlich noch einen Punkt fest ... vielleicht hast du ja eine andere Beschränkung? Es fehlen also Details.

Ja, das stand hier:

Zitat:

Zitat von cltom (Beitrag 1241898)

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.

Zitat:

Zu ein paar anderen Fragen: die Bedingung, dass die gelben und blauen Flächen gleich sind, würde ich auch so auslegen, dass die Differenz der Integrale nicht Null ist, sondern nur möglichst klein.

Du bist da einfach zu zurückhaltend :stupid:
Die Forderung "Differenz möglichst klein" führt zunächst einmal zu einer Differenz gegen unendlich. Die Forderung "Betrag der Differenz möglichst klein" führt dann direkt zu dieser Gleichung:

minimiere abs(a*(8/3 - 2*d^2) - 2*(d-1)*(b-m))

Der Inhalt des Betrags für (beispielsweise d=0) ist dann (a*(8/3) + 2*b - 2*m)) und nimmt damit für m in ℝ ebenfalls Werte aus ganz ℝ an.

Aus dem Zwischenwertsatz folgt damit unmittelbar, dass es einen Wert m gibt, der den Betrag zu 0 werden lässt.

Die Forderung "Betrag minimal" ist also eine Formulierung, die durch "Betrag gleich 0" präzisiert werden kann ohne die Lösungsmenge einzuschränken.

Anschaulich gesprochen: Für eine Gerade mit einer Steigung gegen unendlich wird die Fläche (Parabel-gerade) sehr negativ. Für eine Steigung gegen minus unendlich wird sie sehr positiv. Dazwischen muss eine Steigung existieren, für die diese Differenz 0 wird.

Das ganze gilt allerdings nur falls der Stützpunkt nicht in der Mitte des Intervalls ligt. Denn dann hängt die Fläche der Gerade nicht mehr von der Steigung ab. Das führt auch zu einem neuen, interessanten Kriterium: Das Integral über das Quadrat der Differenz soll minimal werden. (Da alle Untersummen positiv sind, kann es nicht 0 werden.)
In diesem Fall kann man Integration und Differenz aber nicht mehr vertauschen.