© Der/die Autor(en), exklusiv lizenziert an Springer Fachmedien Wiesbaden GmbH, ein Teil von Springer Nature 2022
L. ClassenMit Jupyter durchs Physikpraktikumessentialshttps://doi.org/10.1007/978-3-658-37723-6_6

6. Modellanpassung

Lew Classen1  
(1)
Institut für Kernphysik, WWU Münster, Münster, Deutschland
 
 
Lew Classen
Email: lew.classen@uni-muenster.de

In diesem Kapitel werden Sie ...

  • die Grundzüge der Modellanpassung kennen lernen,

  • die Funktion scipy.optimize.curve_fit() für die Anpassung mit Unsicherheiten in einer Variable oder ohne Unsicherheiten verwenden,

  • Anpassungen mit Unsicherheiten in zwei Variablen mit der Funktion scipy.odr() durchführen.

In vielen Praktikumsexperimenten werden Zusammenhänge zwischen zwei (oder mehr) physikalischen Größen untersucht. Dabei wird oft eine Größe (die unabhängige Variable) variiert und eine andere Größe (abhängige Variable) in Anhängigkeit von der ersten aufgenommen.

6.1 Optimale Parameter

Eine typische Ausgangslage ist ein Datensatz und ein Modell mit freien Parametern, von dem man ausgeht, dass es die Daten sinnvoll beschreiben kann. Das Ziel ist es nun, optimale Modellparameter zu finden, welche die Daten (im Rahmen des Modells) am besten beschreiben. Dabei sind Modellparameter dann besonders gut, wenn sie den Unterschied zwischen Daten und Modell minimieren. Verschiedene Methoden der Modellanpassung unterscheiden sich darin, wie dieser Unterschied definiert wird.

Besonders nachvollziehbar ist dieses Prinzip bei der Methode der kleinsten Quadrate1 (least squares fit). Hier wird der quadrierte vertikale Abstand zwischen Datenpunkten und Modell minimiert. Bei einem Datensatz mit N Wertepaaren $$x_i, y_i$$, die wir mit der Modellfunktion f(x) beschreiben wollen, erhalten wir also die optimalen Parameter, wenn die Größe Q minimal wird:
$$ Q = \sum _{i=0}^{N} (y_{i} - f(x_i))^2 $$
Sind die Unsicherheiten $$\Delta y_i$$ der Messwerte $$y_i$$ bekannt, können diese ebenfalls in die Parametersuche einfließen:
$$ Q = \sum _{i=0}^{N} \frac{(y_{i} - f(x_i))^2 }{\left( \Delta y_i\right) ^2} $$
In diesem Fall werden die tatsächlichen Abweichungen zwischen Modell und Datenwert also zur Unsicherheit des Wertes ins Verhältnis gesetzt. Eine betragsmäßig gleiche Abweichung trägt dann bei großer Unsicherheit weniger zur Summe Q bei, als bei einer kleinen Unsicherheit.
Schauen wir uns dieses Prinzip an einem praktischen Beispiel an: In Abschn. 3.​1 haben wir, auf der Suche nach dem absoluten Temperaturnullpunkt, einen hypothetischen Behälter mit Luft erwärmt und dabei Temperatur und Druck aufgezeichnet2. Wir modellieren die Luft als ideales Gas. Der Druck p(T) als Funktion der Temperatur (in der Einheit ) lässt sich dann als
$$ p(T) = \frac{nR}{V}\cdot (T - T_0) $$
schreiben3, wobei n für die Stoffmenge, V für das Volumen und R für die ideale Gaskonstante stehen. Wir können also eine verschobene Ursprungsgerade als Modellfunktion verwenden und sind besonders am Modellparameter $$T_0$$ interessiert, dem Schnittpunkt der Geraden mit der Temperaturachse bzw. dem absoluten Temperaturnullpunkt. In den folgenden Abschnitten werden wir dieses Modell an den Datensatz anpassen, und zwar ohne und mit Berücksichtigung der Messunsicherheiten.
Das Paket scipy (der Name steht für scientific python) enthält zahlreiche Tools für das wissenschaftliche Arbeiten. Im Unterpaket scipy.optimize sind Funktionen für Optimierungsaufgaben zusammengefasst. Zu diesen gehört auch die Modellanpassung. Sie können dieses Paket folgendermaßen importieren:

6.2 Unbekannte Unsicherheiten oder Unsicherheiten in einer Dimension

Sollen bei der Modellanpassung keine Unsicherheiten oder nur Unsicherheiten in einer Variable berücksichtigt werden, ist die Funktion optimize.curve_fit() das Mittel der Wahl. In der Praxis ist das der Fall, wenn die Unsicherheiten einer oder beider Variablen vernachlässigbar klein oder unbekannt sind. Der Algorithmus benötigt

  • eine Modellfunktion,

  • Daten,

  • Unsicherheiten in einer Variable (optional), sowie

  • Startwerte für die Parametersuche (ebenfalls optional).

Nachdem wir den Datensatz bereits in Abschn. 3.​1 aus einer Textdatei eingelesen und in den Arrays temperatur und druck abgelegt haben, übersetzten wir die Modellfunktion4 in Python-Schreibweise:
Ohne Berücksichtigung von Unsicherheiten lautet der Befehl für die Optimierung:
Sollen Unsicherheiten berücksichtigt werden, übergibt man diese als zusätzlichen Parameter (sigma) an die Funktion. Die Funktion geht hier automatisch davon aus, dass es sich um die Unsicherheiten der abhängigen Variable (y-Werte) handelt. Die Funktion erwartet ein Array mit je einer Unsicherheit pro Datenpunkt, also eines mit derselben Länge wie das Array der y-Werte. Diese sind im Allgemeinen nicht gleich. Hat man es, wie in unserem Beispiel, mit identischen Unsicherheitswerten zu tun, kann man sich mit der numpy-Funktion np.ones_like()5 behelfen, wie wir es hier tun6:

Der Parameter absolute_sigma=True sorgt dafür, dass curve_fit() die Werte als tatsächliche $$1\sigma $$-Unsicherheiten und nicht als eine relative Gewichtung7 interpretiert.

Die Parametersuche ist (bis auf Ausnahmen) ein iterativer Prozess. Man kann ihn sich so veranschaulichen, dass schrittweise bestimmte Kombinationen von Modellparametern getestet werden. Dabei sollte die Summe der quadratischen Abweichungen schrittweise kleiner werden. Ist das nicht der Fall, schlägt der Algorithmus in der Parameterebene einen anderen Weg ein. Ändert sich die Summe nur noch wenig, wird die Suche beendet. Das Ziel der Suche – die optimalen Modellparamter – entsprechen einem Minimum in der Parameterebene. Es kann aber passieren, dass die Suche in einem lokalen Minimum endet und so falsche Ergebnisse liefert. Eine Möglichkeit, das zu verhindern, sind Startwerte für die Suche. Diese sollten grob im Bereich der erwarteten Werte liegen8. In unserem Beispiel können wir die Steigung (durch ein Steigungsdreieck) und den Achsenabschnitt grob abschätzen. Der interaktive Modus von matplotlib, in dem die aktuellen Koordinaten des Mauszeigers immer angezeigt werden, kann hier beim Vermessen der Kurve hilfreich sein. Die Übergabe der Startwerte erfolgt als Liste an den Parameter p0. Die Reihenfolge richtet sich nach der Reihenfolge der Parameter in der von Ihnen definierten Modellfunktion:

Hier schätzen wir also die Steigung über ein Steigungsdreieck mit $$\Delta y =$$ 100 und $$\Delta x =$$ 30 und den Temperaturnullpunkt mit -300 ab. Ohne explizite Startwerte werden alle Modellparameter zunächst auf 1 gesetzt.

Wenn Sie die obigen Code-Beispiele ohne Fehlermeldung ausgeführt haben, ist das ein Zeichen, dass die Modellanpassung prinzipiell funktioniert hat. Allerdings haben wir bisher noch keine Ergebnisse der Parametersuche gesehen. Das liegt daran, dass diese in den obigen Code-Beispielen nicht auf den Schrim ausgegeben, sondern in den Variablen opt und cov abgelegt werden. Schauen wir uns diese nun genauer an. Die Variable opt ist ein Array. Es enthält die optimalen Parameter und hat in unserem Fall9 die Länge 2. Die Variable cov, ebenfalls ein Array, enthält die Kovarianzmatrix. Auf der Hauptdiagonalen stehen die Varianzen der optimalen Parameter. Wir interessieren uns, als Maß für die Unsicherheit auf die Modellparamter, für die Standardabweichungen der Parameter, also die Wurzeln der Diagonaleinträge. Man kann diese beispielsweise folgendermaßen erhalten:
Hier sehen Sie eine Verschachtelung von numpy-Funktionen: Die Wurzelfunktion np.sqrt() wird direkt auf das Ergebnis von np.diag() angewendet, welches die Diagonaleinträge von cov liefert.
Abb. 6.1

Komplette Modellanpassung mit optimize.curve_fit() mit Unsicherheiten in einer Dimension (sigma) und Startwerten für die Modellparameter (p0)

In Abb. 6.1 finden Sie den kompletten Code für einen erfolgreichen Fit mit curve_fit() auf einen Blick. In diesem Fall werden sowohl die Unsicherheiten in der Druckmessung berücksichtigt, als auch Startwerte für die Optimierung der Modellparameter angegeben. Für den gesuchten Temperaturnullpunkt erhalten wir hier den Wert $$T_0 = (-277 \pm 4)\,^{\circ }{\text {C}}$$.

Mit den matplotlib-Funktionen aus Abschn. 4.​1 können wir das Ergebnis auch graphisch darstellen (siehe Abb. 6.2). Die Plausibilität der Unsicherheitsangabe und die Passgenauigkeit eines Modells können Sie übrigens schnell grob überprüfen, indem Sie nachzählen, welcher Anteil der Messwerte im Rahmen der Unsicherheiten mit der Modellfunktion verträglich ist. Handelt es sich um $$1\sigma $$-Unsicherheiten, sollte unsere Gerade nur knapp zwei Drittel der Unsicherheitsbereiche schneiden.
Abb. 6.2

Diagramm mit dem Ergebnis der Modellanpassung mit optimize.curve_fit() aus Abb. 6.1

6.3 Unsicherheiten in zwei Dimensionen

Im allgemeinsten Fall einer Modellanpassung müssen Unsicherheiten in beiden Variablen beachtet werden. Da das im Rahmen der Methode der kleinsten Quadrate nicht möglich ist, kommt hier z. B. die Orthogonale Regression zum Einsatz10. Die entsprechenden Tools bietet das scipy-Unterpaket scipy.odr. Die Vorgehensweise bei der Modellanpassung ähnelt im Prinzip der von optimize.curve_fit(), hat aber ein paar Besonderheiten. Ein Fit mit scipy.odr11 gliedert sich in folgende Schritte:

  • Definition einer Modellfunktion

  • Erstellen des Modells

  • Formatieren der Messdaten, inklusive Unsicherheiten

  • Übergeben der Daten, des Modells und der Startwerte für die Parametersuche an den Algorithmus

  • Ausführen des Fits

  • Abholen der Ergebnisse

In Abb. 6.3 sehen Sie die Anwendung von scipy.odr auf unser obiges Beispiel, die Suche nach dem Temperaturnullpunkt. Die Messreihe ist hier bereits importiert und die Unsicherheiten in den Variablen u_temperatur und u_druck abgelegt. Im Folgenden werden wir uns den Code für die einzelnen Schritte getrennt anschauen. Sie können entsprechend jeden Schritt in einer separaten Code-Zelle ausführen. Das ist aber nicht zwingend: Haben Sie den Algorithmus einmal aufgesetzt, können Sie alle Befehle auch in eine gemeinsame Zelle platzieren, um beim wiederholten Ausführen Zeit zu sparen.
Abb. 6.3

Komplette Modellanpassung mit scipy.odr. Die Datengrundlage ist identisch zu Abb. 6.1

Importieren wir zunächst das Paket:
Anschließend wird die Modellfunktion definiert.
Es handelt sich um die gleiche verschobene Ursprungsgerade wie in dem vorangegangenen Beispiel. Allerdings müssen wir hier der Konvention des Algorithmus folgen und alle Modellparameter in einer Liste (hier B) zusammenfassen. Die Steigung ist also B[0] und der gesuchte Achsenabschnitt B[1]. Für die weitere Verwendung muss die Modellfunktion zu einem odr.Model() konvertiert werden:
Auch die Daten und ihre Unsicherheiten werden speziell formatiert:
Mit den so erzeugten Variablen bzw. Objekten12 mymodel und mydata übergeben wir nun Daten, das Modell und die Startwerte für die Parametersuche an den Algorithmus. Das geschieht mit der folgenden Codezeile:
Anders als bei curve_fit() sind hier Startwerte für die Parameteroptimierung (beta0) verpflichtend. Sie werden als Liste (entsprechend der Reihenfolge der Modellparameter in der weiter oben definierten Modellfunktion) durch den Parameter beta0 übergeben. Damit haben wir nun einen Algorithmus für unser spezielles Optimierungsproblem – den Geradenfit – aufgesetzt. Sie steuern ihn über das Objekt myodr. Nun ist es an der Zeit, dass myodr seine Arbeit tut und eine Modellanpassung durchführt. Mit
geben Sie den Befehl dazu und legen gleichzeitig die Fitergebnisse in der Variable myoutput ab. Diese sind etwas umfangreicher als bei der Modellanpassung mit optimize.curve_fit(). Eine Übersicht der wichtigsten Ergebnisse können Sie mit
auf den Bildschirm ausgeben (siehe Abb. 6.4).
Abb. 6.4

Übersicht der Fitergebnisse von scipy.odr

Der Output enthält – neben den optimalen Parametern Beta und den zugehörigen Unsicherheiten Beta Std Error – auch weitere Informationen. Speziell der Punkt Reason(s) for Halting ist hier relevant. Hat der Algorithmus tatsächlich ein Minimum im Parameterraum gefunden, sind die Ergebnisse also verlässlich, sollte hier Sum of squares convergence stehen. In unserem Beispiel ist das der Fall, so dass wir den gesuchten Temperaturnullpunkt $$T_0 = (- 277 \pm 5)\,^{\circ }{\text {C}}$$ notieren können. Anderenfalls sollte der Fit z. B. mit veränderten Startwerten wiederholt werden. Eine weitere Möglichkeit besteht darin, die Anzahl der Iterationsschritte zu erhöhen. Dafür enthält odr.ODR() den Paramter maxit. Standardmäßig werden maximal 50 Iterationen durchgeführt. Durch
erhöhen Sie das Maximum auf 60.
Abb. 6.5

Ergebnis der Modellanpassung mit scipy.odr

Natürlich können Sie die Modellparameter und ihre Unsicherheiten auch einzeln ausgeben. Der Befehl
liefert ein Array mit den optimalen Modellparametern. Die Reihenfolge der Parameter richtet sich nach Ihrer Definition der Modellfunktion: Die optimale Steigung ist myoutput.beta[0] und der optimale Achsenabschnitt myoutput.beta[1]. Für einen Plot des Modells (siehe Abb. 6.5) können Sie direkt das komplette Array verwenden. Die entsprechenden Unsicherheiten erhalten Sie durch:

Diese Werte können Sie nun für weitere Berechnungen in Variablen ablegen.

Mit dem Paket uncertainties, zu dem Sie in Kap. 5 mehr erfahren, können Sie sogar Wert und Unsicherheit in einer gemeinsamen Variable ablegen und für weitere Berechnungen nutzen:

Und damit haben Sie es geschafft! Nun sind Sie bereit für Ihre erste Versuchsauswertung mit Python. Viel Erfolg dabei!