Schätzen von Parabel-Parametern anhand von Messwerten
-
Hallo,
ich möchte gerne Messwerte, die einer auf dem Kopf stehenden Parabel ähneln, fitten.
Ich nutze dafür die C-Bibliothek lmfit. Hierfür muss man allerdings Startwerte für die Parabel-Parameter übergeben und wenn die Startwerte schlecht sind, dann kann es sein, dass Blödsinn herauskommt.Nun würde ich gerne herausfinden, wie ich an gute Startwerte herankomme.
Wenn die Parabelfunktion folgendermaßen aussieht:
f(x) = a*x^2 + b*x + c
Dann kann ich sagen, dass wenn die Parabel auf dem Kopf steht und alle Funktionswerte positiv sind, dass dann
a
negativ undb
wahrscheinlich positiv ist. Überc
bin ich mir im unklaren.
Das ist schon mal kein schlechter Anfang, aber das reicht nicht immer aus. Kann man das noch besser eingrenzen, auch was die Größenordnung angeht?Danke im Voraus.
-
@Steffo
Das ist doch eine klassiche Extremwertaufgabe.Du sucht dir das Minimum deiner Messwerte. Sagen wir mal dieses läge bei (xm, ym)
Das Minimum einer Funktion ist definiert durch: f`(x) == 0 und f``(x) > 0 .
f'(x) = 2*a*x+b f``(x) = 2*a f'(x) = 0 -> 2*a*x+b = 0 x ist in diesem Fall xm -> 2*a*xm + b = 0 -> b = -2*a*xm f''(x) > 0 -> 2*a > 0
Ich würde a > 0 wählen, b = -2 * a * xm ausrechnen und c auf ym setzen.
-
@Quiche-Lorraine
Danke für deine Antwort. Ich habe bemerkt, dass lmfit Overkill ist, da das dafür gemacht ist alle möglichen von Funktionen zu fitten (du kannst z. B. eine Gauß-Funktion übergeben). Aber für beliebige Polynome, wie eben eine Parabel, kann ich eine multiple lineare Regression anwenden, was ich nun in C++ nachimplementiert habe und ich sehr gute Ergebnisse bekomme!!!
-
Ich benutze für sowas aus Nostaliegründen (und weil ich damit klarkomme) immer ROOT:
root [0] double x[] = {1,2,3,4,5}; root [1] double y[] = {-5,1,2,0.5,-6}; root [2] auto gr = new TGraph(std::size(x), x, y); root [3] gr->Fit("pol2") **************************************** Minimizer is Linear Chi2 = 0.357143 NDf = 2 p0 = -14.5 +/- 0.906327 p1 = 11.5357 +/- 0.690681 p2 = -1.96429 +/- 0.112938 (TFitResultPtr) <nullptr TFitResult>
Braucht hier auch keine Startwerte.
-
Hat jemand einen Tipp, wie ich die numerische Stabilität erhöhen kann? Ich sehe, dass das für numerisch problemlose Datensätze wunderbar funktioniert, aber sobald der x-Wert sehr klein wird und der y-Wert sehr hoch, dann bekomme ich Quatsch heraus.
Was gibt es da für Möglichkeiten? Vielleicht Daten zentrieren, in dem ich den Mittelwert abziehe? Evtl. normieren? Wie kann ich sicherstellen, dass dann immer noch die richtige Lösung herauskommt?@wob Du scheinst da Erfahrungen zu haben.
-
@Steffo sagte in Schätzen von Parabel-Parametern anhand von Messwerten:
Was gibt es da für Möglichkeiten? Vielleicht Daten zentrieren, in dem ich den Mittelwert abziehe? Evtl. normieren?
Ja, das kann sehr gut helfen. Gerade wenn du mein Beispiel z.B. um 100.000 nach links und 100.000 nach oben schiebst - dann funktioniert der Fit so nicht mehr richtig. Daher: machen! Also verschieben ja, normieren: weiß nicht so recht. Kommt drauf an, eher nicht. Aber ich kenne deine Werteranges nicht. Und ich habe immer gerne ausprobiert.
In manchen Fällen kann es auch sinnvoll sein, Constraints anzugeben - also dass eine Variable nur zwische a und b liegen darf. Das ist allerdings oft gefährlich, wenn das intern z.B. durch eine Transformationsfunktion gelöst ist, die dann dazu führt, dass der fit ungenau wird.
Wie kann ich sicherstellen, dass dann immer noch die richtige Lösung herauskommt?
Grafisch angucken Oder auf das chi2/ndf bzw. die Unsicherheiten des Fits gucken - dann siehst du zumindest mal, wenn der nicht funktioniert hat.
-
@wob Ich bin gerade Zuhause und spiele hier mit python und numpy herum.
Ich scheine Quatsch herauszubekommen, wenn ich x vorher zentriere (so wie das in der Doku von numpy.polyfit empfohlen wird) und anschließend anhand der gefitteten Parameter über die Parabel-Ableitung den Extrempunkt berechne und dann wieder den Mittelwert addiere.polyfit issues a RankWarning when the least-squares fit is badly conditioned. This implies that the best fit is not well-defined due to numerical error. The results may be improved by lowering the polynomial degree or by replacing x by x - x.mean().
Hier meine Berechnung für das Extrema:
f(x) = a*x^2 + b*x + c f'(x) = 2*a*x + b = 0 --> x = -b / (2*a)
Wenn ich den x-Vektor vorher zentriert habe, dachte ich, dass beim Extremum einfach x_mean drauf addieren muss:
x = (-b / (2*a)) + x_mean
Aber leider kommt nicht dasselbe Ergebnis heraus als wenn ich vorher nicht den x-Vektor zentriert habe (979654289,4994783 statt 1032758932).
-
@wob
EDIT: Ich hatte einen Klammerfehler in meinem Code. Nun funktioniert es. Habe den Klammerfehler weiter oben im Beitrag nun korrigiert.