SciPy¶
- Baut auf NumPy auf
- Kann numerisch integrieren, DGLs lösen, optimieren, minimieren, …
- Enthält auch physikalische Konstanten und wichtige mathematische Funktionen
Überblick über (fast) alles Weitere, was SciPy kann, gibt es hier:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
Inhalt¶
Abrufen von Naturkonstanten¶
Wenn solche Konstanten genutzt werden, muss das korrekt mitgeteilt, also zitiert werden. Darauf gehen wir nächste Woche im LaTeX-Workshop ein :-)
(Quelle hier: scipy + version)
import numpy as np
# physical constants
import scipy.constants as const
planck_value = const.Planck # returns the value
planck_full = const.physical_constants[
"Planck constant"
] # returns value, unit and uncertainty
print(f"{planck_value = }", f"{planck_full = }", sep="\n")
# a list of all included constants
# search the web for "scipy physical constants reference guide" to find the same list
# We don't print the full list in the html file, to keep the notebook more managable
# const.physical_constants
# convert temperatures:
print(const.convert_temperature(100, "c", "K"))
print(const.convert_temperature(100, "kelvin", "Celsius"))
print(const.convert_temperature(100, "f", "C"))
# convert angles:
print(np.rad2deg(np.pi))
print(np.deg2rad(90))
Fitten¶
Oft möchte man eine Funktion mit freien Parametern, zum Beispiel eine Erwartung aus der Theorie, an die gemessenen Werte anpassen. Dies nennt man Fit.
Die Funktion scipy.optimize.curve_fit nutzt die numerische Methode der kleinsten Quadrate, die arbiträre Funktionen fitten kann.
Für Funktionen, die eine Linearkombination von Einzelfunktionen sind, also
$$ f(x) = \sum_i^N a_i \cdot f_i(x), $$
existiert eine analytische Lösung. Deswegen sollten in solchen Fällen (z.B. alle Polynome) entsprechende Funktionen genutzt werden (z.B. np.polyfit).
Lineare Regression bzw. Regression von Polynomen¶
Im Folgenden wird eine lineare Regression mit np.polyfit durchgeführt:
$$ f(x) = a + b \cdot x $$
# prepare plot
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (10, 8)
plt.rcParams["font.size"] = 16
# load data
x, y = np.genfromtxt("data/example_data_linear.txt", unpack=True)
fig, ax = plt.subplots(1, 1, layout="constrained")
ax.plot(x, y, "k.", label="example data")
ax.set(xlabel=r"$t \,/\, \mathrm{s}$", ylabel=r"$s \,/\, \mathrm{m}$");
# Fit a polynomial of degree 1, return covariance matrix
params, covariance_matrix = np.polyfit(x, y, deg=1, cov=True)
uncertainties = np.sqrt(np.diag(covariance_matrix))
for name, value, uncertainty in zip("ab", params, uncertainties):
print(f"{name} = {value:.3f} ± {uncertainty:.3f}")
x_plot = np.linspace(0, 10)
fig, ax = plt.subplots(1, 1, layout="constrained")
ax.plot(x, y, "k.", label="Messwerte")
ax.plot(
x_plot,
params[0] * x_plot + params[1],
label="Lineare Regression",
linewidth=3,
color="tab:orange",
)
ax.legend()
ax.set(xlabel=r"$t \,/\, \mathrm{s}$", ylabel=r"$s \,/\, \mathrm{m}$");
Nichtlineare Funktionen der Parameter¶
Wenn eine Funktion nicht linear bezüglich der freien Parameter ist, muss die Lösung numerisch gefunden werden.
Als erstes Beispiel dafür schauen wir uns folgende Daten an:
x, y = np.genfromtxt("data/sin_fit.txt", unpack=True)
ax.cla()
ax.plot(x, y, "o", label="Messwerte")
ax.set(xlabel=r"$t$ / s", ylabel=r"$U$ / V")
ax.legend()
fig
Diese Daten sollen an in der Theorie eine Sinus-Kurve beschreiben. Wir schreiben den Sinus allgemein als $$U(t) = a \cdot \sin(b \cdot t) + c~.$$ und definieren uns die entsprechende Funktion im Code.
def sinus(x, a, b, c):
return a * np.sin(b * x) + c
Für das fitten nutzen wir jetzt curve_fit und bekommen wieder die Parameter und Kovarianz-Matrix als Rückgabewerte.
from scipy.optimize import curve_fit
params, covariance_matrix = curve_fit(sinus, x, y)
uncertainties = np.sqrt(np.diag(covariance_matrix))
for name, value, uncertainty in zip("abc", params, uncertainties):
print(f"{name} = {value:.3f} ± {uncertainty:.3f}")
Jetzt können wir die gefittete Funktion in den Plot der Messwerte hinzufügen.
ax.cla()
x_sin_fit = np.linspace(0, 2 * np.pi)
ax.plot(x, y, "o", label="Messwerte")
ax.plot(x_sin_fit, sinus(x_sin_fit, *params), label="Sinus Fit")
ax.set(xlabel=r"$t$ / s", ylabel=r"$U$ / V")
ax.legend()
fig
Bei dieser Variante des Fits kann es passieren, dass der Minimierungsalgorithmus in lokale Maxima hineinläuft und unsinnige Ergebnisse liefert, dies kann mit guten Startwerten meistens vermieden werden.
Als Beispiel einer komplexeren Funktion wird im Folgenden die Sigmoidfunktion verwendet (ähnlich zum tanh):
$$ f(x; a, b, c) = \frac{a}{1 + \exp(-(x-b))} + c$$
def sigmoid(x, a, b, c):
return a / (1 + np.exp(-(x - b))) + c
x_plot = np.linspace(-50, 50, 1000)
ax.cla()
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.plot(x_plot, sigmoid(x_plot, 1, 0, 0), label="Sigmoid")
ax.plot(x_plot, np.tanh(x_plot), label="tanh")
ax.legend()
ax.set(xlabel=r"$x$", ylabel=r"$f(x)$")
fig
Die Messwerte aus einem Praktikumsversuch:
x, y = np.loadtxt("data/fit_data_with_init_values.txt", unpack=True)
ax.cla()
ax.plot(x, y, "o", label=r"Messwerte")
ax.set(xlabel=r"Temperatur / °C", ylabel=r"$GP$")
fig
Ein einfacher Fit wie oben funktioniert hier nicht so gut:
params, covariance_matrix = curve_fit(sigmoid, x, y)
uncertainties = np.sqrt(np.diag(covariance_matrix))
for name, value, uncertainty in zip("abc", params, uncertainties):
print(f"{name} = {value:8.3f} ± {uncertainty:.3f}")
Schaut man sich die berechnete Ausgleichskurve an sieht man auch, dass das nicht stimmen kann:
ax.cla()
ax.plot(x, y, "o", label="Messdaten")
ax.plot(x_plot, sigmoid(x_plot, *params), "-", label=r"Sigmoid Fit")
ax.legend()
ax.set(xlabel=r"Temperatur / °C", ylabel=r"$GP$")
axins = ax.inset_axes(
[0.60, 0.4, 0.3, 0.3],
xlim=(24, 46), ylim=(0.3, 0.9)
)
axins.plot(x, y, "o")
axins.plot(x_plot, sigmoid(x_plot, *params), "-")
ax.indicate_inset_zoom(axins, edgecolor="black", alpha=0.4)
fig
Was macht man jetzt?
Bei solchen Fragen hilft die Dokumentation der Pythonmodule (hier: scipy) oder Stack Overflow weiter.
Folgendes Google-Muster ist ein guter Anfang (beachte englische Sprache):
python <module-name> <function-name> <What went wrong?>Also in diesem Fall: python scipy curve_fit fails
Damit dieser Fit funktioniert müssen die Startwerte für den internen Minimierungsalgorithmus angepasst werden.
Aus der Dokumentation/Stack Overflow wissen wir jetzt, dass man mit dem
keyword argument p0 (Standardwert is p0=(1,1,1)) die Startwerte einstellt:
params, covariance_matrix = curve_fit(sigmoid, x, y, p0=(-1, 40, 1))
uncertainties = np.sqrt(np.diag(covariance_matrix))
for name, value, uncertainty in zip("abc", params, uncertainties):
print(f"{name} = {value:8.3f} ± {uncertainty:.3f}")
ax.cla()
x_plot = np.linspace(0, 50, 1000)
ax.plot(x, y, "o", label="Messwerte")
ax.plot(x_plot, sigmoid(x_plot, *params), "-", label="Sigmoid Fit")
ax.legend()
ax.set(xlabel=r"Temperatur / °C", ylabel=r"$GP$")
fig
Zum Vergleich der beiden Anfangswerte (seeds) kann man sich die einmal ansehen und mit den angepassten Parametern vergleichen:
default_seed = (1, 1, 1)
good_seed = (-1, 40, 1)
parameter = [default_seed, good_seed, params]
x_plot = np.linspace(-80, 80, 1000)
ax.cla()
for a, b, c in parameter:
ax.plot(
x_plot, sigmoid(x_plot, a, b, c), label=f"f(x; {a:0.3f}, {b:0.3f}, {c:0.3f})"
)
ax.legend()
ax.set(xlabel=r"$x$", ylabel=r"$f(x)$")
fig
Die richtigen Startwerte findet man entweder durch
trial and error => einfach ausprobieren bis es klappt
nachdenken => siehe unten
Im obigen Beispiel musste nur der Parameter b angepasst werden,
weil der für die Form der Kurve sehr wichtig ist.
$$ f(x; a, b, c) = \frac{a}{1 + \exp(-(x-b))} + c$$
bs = [0, 20, 40]
x_plot = np.linspace(-50, 50, 1000)
ax.cla()
for b in bs:
(line,) = ax.plot(x_plot, sigmoid(x_plot, 1, b, 0), label=f"f(x; 1, {b}, 0)")
ax.plot(
b,
sigmoid(b, 1, b, 0),
"o",
color=line.get_color(),
ms=20,
label=f"f(x={b}) = {sigmoid(b, 1, b, 0)}",
)
ax.legend()
ax.set(xlabel=r"$x$", ylabel=r"$f(x)$")
fig
Der Parameter $b$ gibt den $x$-Wert an, bei dem die Funktion auf die Hälfte des Maximums abgefallen ist.
Bei den Messwerten oben ist die Stelle ungefähr bei $x=40$, also ist b = 40 ein guter Startwert.
Das lässt sich auch automatisieren:
ax.cla()
ax.plot(x, y, "o")
idx = np.argmin(np.abs(y - 0.5))
ax.plot(x[idx], y[idx], "o")
print(x[idx], y[idx])
ax.set(xlabel=r"$x$", ylabel=r"$f(x)$")
fig
Weitere nützliche Funktionen für das Praktikum¶
Folgende Eigenschaften von Scipy sollen nur die Idee vermitteln, was noch alles möglich ist.
Statistische Verteilungen und Funktionen¶
Das scipy.stats Modul enthält viele Wahrscheinlichkeitsverteilungen und -funktionen. Als Beispiel wird hier jedoch nur die Standardabweichung des Mittelwertes berechnet.
import numpy as np
x = np.array([2, 5, 7]) # create test array
from scipy.stats import sem
print(sem(x)) # standard error of the mean
Finden von Peaks¶
Das scipy.signal Modul enthält Funktionen zur Signalverarbeitung. Das kann interessant sein zum automatischen Finden von Peaks, oder Bestimmung der Periodizität.
# example data
from scipy.datasets import electrocardiogram
x = electrocardiogram()[2000:4000]
# find peaks
from scipy.signal import find_peaks
peaks, _ = find_peaks(x, distance=150)
fig, ax = plt.subplots(1, 1, layout="constrained")
ax.plot(x, label="EKG")
ax.plot(peaks, x[peaks], "x", ms=10, label="peaks")
ax.set(xlabel=r"$t$", ylabel=r"$I \,/\, a.u.$")
ax.legend();
