Funktionen schnell an Messwerte anfitten

Wissenschaft Physik
1

Moien

Wie kann man eine Funktion schnell an Messwerte anfitten ?

Ich hab etwa 50000 Datensätze mit jeweils 100 Messpunkten. Angefitten werden soll einmal eine Gaussglocke, einmal eine atan Funktion. Die meisten der Datensätze werden nicht sonderlich gut passen.

Mein Ansatz für Gauss wäre:

  • das Maximum suchen, seine Position als Mittelpunkt zu nehmen
  • die ersten paar Werte von links und rechts mitteln, als Basislinie nehmen
  • Simplexoptimierung von da an.

Bei atan würde ich Mittelwerte der ersten paar Werte von links und rechts bestimmen dann wieder Simplex fahren.

Geht das irgendwie schneller ?

Danke

2

Moien

Wie kann man eine Funktion schnell an Messwerte anfitten ?

Normalerweise mit einem least squarers Ansatz. Du bildest bei jedem x_i das Quadrat der Differenz zwischen Messwert und Funktionswert. Das alles summiert gibt eine Zielfunktion von der man das Minimum über die Parameter der Funktion sucht (=> Newtonverfahren).
Beispiel: Du willst die Funktion f(x_i)=arctan(p*x_i) fitten.
Deine Zielfunktion ist dann z§=SUM_i(y_i-arctan(p*x_i))²,
also die Summe über i der quadrierten Differenzen von Messwert und Funktionswert. Das lässt sich leicht nach p ableiten. Das Newtonverfahren benutzt dann die Iterationsvorschrift p_k+1=p_k-z(p_k)/z’(p_k), mit einem möglichst guten Startwert.
Sinkt p_k+1-p_k unter einen bestimmten Toleranzwert wird die Iteration abegebrochen.

Viel Erfolg !

hendrik

3

Da der arctan streng monoton steigend ist reicht es auch z§=SUM_i(y_i-p*x_i)² zu minimieren. Da das eine nach oben geöffnete Parabel ist sucht man also einfach den Scheitel.

hendrik

[Bei dieser Antwort wurde das Vollzitat nachträglich automatisiert entfernt]

4

Hallo Pumpkin!

Moien

Wie kann man eine Funktion schnell an Messwerte anfitten ?

Ich hab etwa 50000 Datensätze mit jeweils 100 Messpunkten.
Angefitten werden soll einmal eine Gaussglocke, einmal eine
atan Funktion. Die meisten der Datensätze werden nicht
sonderlich gut passen.

Das hört sich interessant an, jedenfalls sind 50000 Fits kein Pappenstiel. Was hast du denn als Programm/Tool zur Verfügung?
Wenn du’s selbst programmiert schau auf jeden Fall mal, was die „Numerical Recipies“ so bieten (Kap. 15): http://www.nrbook.com/b/bookcpdf.php

Mein Ansatz für Gauss wäre:

  • das Maximum suchen, seine Position als Mittelpunkt zu nehmen
  • die ersten paar Werte von links und rechts mitteln, als
    Basislinie nehmen
  • Simplexoptimierung von da an.

Bei atan würde ich Mittelwerte der ersten paar Werte von links
und rechts bestimmen dann wieder Simplex fahren.

Geht das irgendwie schneller ?

Danke

Dein Ansatz scheint prinzipiell ok. Ein paar Gedanken noch dazu:

  • Gute Startwerte sind speziell bei Gauß wichtig, exp(-x^2) geht sonst schnell in den Underflow und dann tut sich beim Fitten nichts mehr …

  • Wenn du eine reinen Gauß hast (ohne Untergrund) könntest Du Mittelwert und Standardabw., Summe relativ schnell berechnen unt hättest gute Startwerte.

  • Gauß mit flachem Untergrund: Untergrund näherungsweise bestimmen (wie du geschrieben hast), abziehen, den Rest (so er positiv ist) logarithmieren und mit lin. Fit Parabel bestimmen. Aus der Parabel erhältst du wieder gute Startwerte für den Gauß-Fit.

  • Für y=A*atan(B*(x-C)) den Nulldurchgang suchen (-> C), dort lokal die Steigung bestimmen (-> A*B). A hast du ja schon aus den Limits.

  • Simplex ist robust, aber evtl. nicht das schnellste Verfahren. Die meisten fertigen Fit-Programme nehmen wohl heute den Levenberg-Marquardt-Algorithmus oder eine Variante davon.
    http://de.wikipedia.org/wiki/Levenberg-Marquardt-Alg…

Ich würde es zunächst damit probieren (gibt’s auch bei Num. Rec.!)

Viel Spaß!
Gruß Kurt
*der neugierig ist, was du da eigentlich machst bzw. wozu man 50000 Fits braucht …*

5

Moien

Das hört sich interessant an, jedenfalls sind 50000 Fits kein
Pappenstiel. Was hast du denn als Programm/Tool zur Verfügung?

Alles was in das Budget von 2500 Euro past … und Matlab.

Ich bin noch am planen, das eigentliche Umsetzten kommt später. Aber ich möchte die Laufzeiten von ein paar der Operationen vorher abschätzen können und mal auf den vorhandenen Rechnern testen. Denn die 2500 müssen auch für Zusatzhardware reichen.

Es ist ein Projekt an einer FH wo ich (als Informatiker ohne Ahnung von Physik) für einem Physiker (ohne Ahnung von Informatik) ein paar Ideen umsetzen soll. Mit normaler Physik komm ich eigentlich gut klar, aber das hier geht zu weit ins Mikroskopische.

Wenn du’s selbst programmiert schau auf jeden Fall mal, was
die „Numerical Recipies“ so bieten (Kap. 15):

Ja, bin dabei. Allerdings scheint keiner der Algos wirklich 100% zu passen. Die Ideen sind richtig gut, wie immer bei den NR, aber neu implementiert werden muss das Zeug wohl doch.

Dein Ansatz scheint prinzipiell ok. Ein paar Gedanken noch
dazu:

  • Gute Startwerte sind speziell bei Gauß wichtig, exp(-x^2)
    geht sonst schnell in den Underflow und dann tut sich beim
    Fitten nichts mehr …

Ich kann einen relativ engen Suchraum angeben. Wenn die Kurve da raus läuft macht es physikalisch keinen Sinn mehr (angeblich).

  • Wenn du eine reinen Gauß hast (ohne Untergrund) könntest Du
    Mittelwert und Standardabw., Summe relativ schnell berechnen
    unt hättest gute Startwerte.

Schön wärs…

  • Gauß mit flachem Untergrund: Untergrund näherungsweise
    bestimmen (wie du geschrieben hast), abziehen, den Rest (so er
    positiv ist) logarithmieren und mit lin. Fit Parabel
    bestimmen. Aus der Parabel erhältst du wieder gute Startwerte
    für den Gauß-Fit.

Sehr gute Idee, das teste ich mal.

  • Für y=A*atan(B*(x-C)) den Nulldurchgang suchen (-> C),
    dort lokal die Steigung bestimmen (-> A*B). A hast du ja
    schon aus den Limits.

Nulldurchgang suchen wird etwas komplexer.

Das sind Messwerte einer 2D Abtastung, also eine 3D Punktwolke (Position in X und Y + Intensität in Z). Da den Nulldurchgang finden wird haarig.

Die atan-Fläche soll die Richtung der Kanten genau bestimmen. Normale Kantensuche mit Roberts & Co liefert liefert bei so kleinen Punktwolken ja immer nur 0°, 45° oder 90°.

  • Simplex ist robust, aber evtl. nicht das schnellste
    Verfahren. Die meisten fertigen Fit-Programme nehmen wohl
    heute den Levenberg-Marquardt-Algorithmus oder eine Variante
    davon.

Ja, könnte ich machen. Der zusätzliche Implementierungsaufwand scheint mir minimal.

*der neugierig ist, was du da eigentlich machst bzw. wozu man
50000 Fits braucht …*

Wenn ich die Physik hinter der Messung verstanden hab kann ich’s dir erzählen. Bis jetzt ich mir ziemlich unklar was die Werte eigentlich aussagen… nur der Physiker weis genau was für Fits von welchen Funktionen er braucht.

cu

6

Hallo Pumpkin!

Es ist ein Projekt an einer FH wo ich (als Informatiker ohne
Ahnung von Physik) für einem Physiker (ohne Ahnung von
Informatik) ein paar Ideen umsetzen soll.

Na, mit einem Informatiker, NR, C und notfalls Mathlab muss das zu machen sein. Ein selbst gestricktes C-Progr. für die Startwerte und der Lev-Marq.-Alg. von NR sollten es eigentlich tun (*). Aber was machst Du dann am Schluß mit den 2500 EURO :smile:)

(* wenn die 50000 Fits nicht in 1 Min. fertig sein müssen - wieviel Rechenzeit darf’s denn brauchen ?)

Das sind Messwerte einer 2D Abtastung, also eine 3D Punktwolke
(Position in X und Y + Intensität in Z). Da den Nulldurchgang
finden wird haarig.

Die atan-Fläche soll die Richtung der Kanten genau bestimmen.
Normale Kantensuche mit Roberts & Co liefert liefert bei so
kleinen Punktwolken ja immer nur 0°, 45° oder 90°.

Wenn du nur die Richtung willst: Formulier den Problem um, fitte den tan(alpha) und rechne später daraus alpha aus. Das wäre einfacher, weil linear (oder „linearer“).
Aber Kantendetektion und atan-Fit, das hört sich so an, als hätte da einer die Kante in der Intensitätsverteilung vorher differenziert.
Tip, falls dem so ist:
Überrede deinen Physiker, vorher nicht zu differenzieren. Fitte direkt an die Stufe eine passende Funktion. Differenziern erhöht statistische Fluktuationen, macht korrelierte Fehler , bringt nur Ärger. Direkter Fit ist besser. Hab das mal (mit einem FH-Diplomand) gemacht. Wir haben da das Integral der Gauß-Fkt. drangefittet. Wir hatten auch 2-d „Bilder“. Die Kante kann man ja erst zeilenweise bestimmen, anschließend per Geradenfit die Richtung der Kanten.

Das Gauß-Integral kann man zwar nicht analytisch angeben, aber für den Fit gibt es sehr gute Näherungsfunktionen, die in jedem einschlägigen Programmpaket enthalten sind (wir hatten IDL verwendet, in NR gibt es das natürlich auch). Zum Probieren tut es gnuplot (hat auch die ganzen Funktionen, einschl. Fit mit Lev.-Marq.). Es geht sogar mit Excel und dem Solver. Das habe ich kürzlich einem Doktoranden empfohlen. Der hat das Gauss-Integral dann an seine Testerg. (auch mit einer Art verschmiereter Stufe in den Daten) mit dem Excel-Solver gefittet und das hat prima geklappt. Für 50000 Fits ist das aber nichts!!!
Gruß Kurt

7

Hallo,

Das hört sich interessant an, jedenfalls sind 50000 Fits kein
Pappenstiel. Was hast du denn als Programm/Tool zur Verfügung?

Alles was in das Budget von 2500 Euro past … und Matlab.

Meine Erfahrungen mit der Geschwindigkeit von Matlab sind mehr als ernüchternd, insbesondere (aber nicht nur) bei I/O-Operationen. Deswegen würde ich im Zweifelsfall C nehmen und das Geld in Hardware stecken :wink: - kommt aber natürlich auch auf die anderen Randbedingungen an.

Grüße,
Moritz

8

Hi,

Wie kann man eine Funktion schnell an Messwerte anfitten ?

Ich verwende für meine fits inzwischen fast ausschliesslich python (das scipy.optimize.leastsq modul) und bin sehr zufrieden damit. Ich muss ne ganze Menge Lorentz Kurven fitten (ungefähr 200 Punkte pro Datensatz), und schaff mehr als 10 fits pro Sekunde (ich habs nie gemessen, dürfte aber so ungefähr stimmen).
Falls du python noch nicht kennst kann ich dir empfehlen gleich die enthought python edition zu installieren, da hast du gleich ne ganze Menge mathematischer module dabei (ich benutz vor allem numpy, scipy und pylab/matplotlib).

Mein Ansatz für Gauss wäre:

  • das Maximum suchen, seine Position als Mittelpunkt zu nehmen
  • die ersten paar Werte von links und rechts mitteln, als
    Basislinie nehmen

Ich würde ausserdem noch die Halbwertsbreite bestimmen.

Falls ihr übrigens auch einen Imaginärteil messt, von dem ihr die Funktion kennt, ist es sehr ratsam diesen gleichzeitig mitzufitten. Im Fall meiner Lorentzkurven verbessert das das Ergbenis erheblich.

Gruss, Johannes

9

Moien

Alles was in das Budget von 2500 Euro past … und Matlab.

Meine Erfahrungen mit der Geschwindigkeit von Matlab sind mehr
als ernüchternd, insbesondere (aber nicht nur) bei
I/O-Operationen.

Meine auch. Aber zum durch- und quertesten von Implementierugen ist das Teil richtig gut. Bei Verfahren die verdammt komplex zu implementieren sind hat mir das Ding schon mehrmals den Arsch gerettet. Es ist halt einfacher eine m-Datei zu debuggen als ein Stück C-Code.

cu

10

Moien

(* wenn die 50000 Fits nicht in 1 Min. fertig sein müssen -
wieviel Rechenzeit darf’s denn brauchen ?)

Wie immer: So wenig wie möglich. Aber mit deinem linear-voroptimieren Ansatz für den Gauss komm ich unter 10 min für die Gaussfit. Das ist ein Quantensprung gegenüber dem Code meines Vorgängers (~2 Tage, k.A. was das Teil tut). Um das atan-Teil kümmer ich mich wenn der Physiker hier ist.

Die atan-Fläche soll die Richtung der Kanten genau bestimmen.
Normale Kantensuche mit Roberts & Co liefert liefert bei so
kleinen Punktwolken ja immer nur 0°, 45° oder 90°.

Wenn du nur die Richtung willst: Formulier den Problem um,
fitte den tan(alpha) und rechne später daraus alpha aus. Das
wäre einfacher, weil linear (oder „linearer“).
Aber Kantendetektion und atan-Fit, das hört sich so an, als
hätte da einer die Kante in der Intensitätsverteilung vorher
differenziert.

Ich komm im Moment nicht an dem Man und weiss nicht was er vor hat. Übernächste Woche fängt das Projekt offizel an, dann seh ich klarer.

Überrede deinen Physiker, vorher nicht zu differenzieren.
Fitte direkt an die Stufe eine passende Funktion.
Differenziern erhöht statistische Fluktuationen, macht
korrelierte Fehler , bringt nur Ärger. Direkter Fit ist
besser. Hab das mal (mit einem FH-Diplomand) gemacht. Wir
haben da das Integral der Gauß-Fkt. drangefittet. Wir hatten
auch 2-d „Bilder“. Die Kante kann man ja erst zeilenweise
bestimmen, anschließend per Geradenfit die Richtung der
Kanten.

Danke, ich behalts im Hinterkopf.

11

Moien

Falls du python noch nicht kennst kann ich dir empfehlen
gleich die enthought python edition zu installieren, da hast
du gleich ne ganze Menge mathematischer module dabei (ich
benutz vor allem numpy, scipy und pylab/matplotlib).

Ich kenne python schon, aber unser Cluster nicht. Zum lokalen testen wärs OK.

Falls ihr übrigens auch einen Imaginärteil messt, von dem ihr
die Funktion kennt, ist es sehr ratsam diesen gleichzeitig
mitzufitten. Im Fall meiner Lorentzkurven verbessert das das
Ergbenis erheblich.

Da muss ich nachfragen. Bis jetzt war davon nicht die Rede.

Danke.

12

Danke an alle, Problem soweit gelöst (owt)
.