Überblick über (fast) alles weitere was SciPy kann gibt es hier:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
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
const.epsilon_0
# a list of all included constants
# search the web for "scipy physical constants reference guide" to find the same list
const.physical_constants
# convert temperatures:
print(const.convert_temperature(100, 'c', 'K'))
print(const.convert_temperature(100, 'kelvin', 'Celsius'))
# convert angles:
print(np.rad2deg(np.pi))
print(np.deg2rad(90))
const.physical_constants["proton mass"]
# value, unit, error
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
existiert eine analytische Lösung. Deswegen sollten in solchen Fällen (z.B. alle Polynome) entsprechende Funktionen genutzt werden (z.B. np.polyfit
)
Im folgenden wird eine lineare Regression mit np.polyfit
durchgeführt:
# 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)
plt.plot(x, y, 'k.', label="example data")
# Fit a polynomial of degree 1, return covariance matrix
params, covariance_matrix = np.polyfit(x, y, deg=1, cov=True)
errors = np.sqrt(np.diag(covariance_matrix))
for name, value, error in zip('ab', params, errors):
print(f'{name} = {value:.3f} ± {error:.3f}')
x_plot = np.linspace(0, 10)
plt.plot(x, y, '.', label="Messwerte")
plt.plot(
x_plot,
params[0] * x_plot + params[1],
label='Lineare Regression',
linewidth=3,
)
plt.legend(loc="best")
Wenn eine Funktion nicht linear bezüglich der freien Parameter ist, muss die Lösung numerisch gefunden werden.
Hierbei kann es sein, 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
):
def sigmoid(x, a, b, c):
return a / (1 + np.exp(-(x - b))) + c
x_plot = np.linspace(-50, 50, 1000)
plt.xlabel('x')
plt.ylabel('y')
plt.plot(x_plot, sigmoid(x_plot, 1, 0, 0), label="Sigmoid")
plt.plot(x_plot, np.tanh(x_plot), label="tanh")
plt.legend()
Die Messwerte aus einem Praktikumsversuch:
x, y = np.loadtxt('data/fit_data_with_init_values.txt', unpack=True)
plt.plot(x, y, 'o', label=r'Messwerte')
plt.xlabel('Temperatur / °C')
plt.ylabel('$GP$')
Ein einfacher Fit wie oben funktioniert hier nicht so gut:
from scipy.optimize import curve_fit
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:
plt.xlabel('Temperatur / °C')
plt.ylabel('$GP$')
plt.plot(x, y, 'o', label='Messdaten')
plt.plot(x_plot, sigmoid(x_plot, *params), "-", label=r'Sigmoid Fit')
plt.legend(loc='best')
Was macht man jetzt?
Bei solchen Fragen hilft die Dokumentation der Pythonmodule (hier: scipy) oder Stackoverflow 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/Stackoverflow 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}')
plt.xlabel('Temperatur / °C')
plt.ylabel('$GP$')
x_plot = np.linspace(0, 50, 1000)
plt.plot(x, y, 'o', label='Messwerte')
plt.plot(x_plot, sigmoid(x_plot, *params), "-", label='Sigmoid Fit')
plt.legend(loc='best')
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)
for a, b, c in parameter:
plt.plot(x_plot, sigmoid(x_plot, a, b, c), label=f"f(x; {a:0.3f}, {b:0.3f}, {c:0.3f})")
plt.legend()
Die richtigen Startwerte findet man entweder durch
trial and error => einfach ausprobieren bis es klappt
nachdenken => siehe unten
Im obigen Beispiel musste nur Parameter b
angepasst werden,
weil der für die Form der Kurve sehr wichtig ist.
bs = [0, 20, 40]
x_plot = np.linspace(-50, 50, 1000)
plt.xlabel('x')
plt.ylabel('y')
for b in bs:
line, = plt.plot(x_plot, sigmoid(x_plot, 1, b, 0), label=f"f(x; 1, {b}, 0)")
plt.plot(
b,
sigmoid(b, 1, b, 0),
"o",
color=line.get_color(),
ms=20,
label=f"f(x={b}) = {sigmoid(b, 1, b, 0)}"
)
plt.legend()
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:
plt.plot(x, y, 'o')
idx = np.argmin(np.abs(y - 0.5))
plt.plot(x[idx], y[idx], 'o')
x[idx], y[idx]
Folgende Eigenschaften von Scipy sollen nur die Idee vermitteln, was noch alles möglich ist.
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
Das scipy.signal
Modul enthält Funktionen zur Signalverarbeitung
# example data
from scipy.misc import electrocardiogram
x = electrocardiogram()[2000:4000]
# find peaks
from scipy.signal import find_peaks
peaks, _ = find_peaks(x, distance=150)
plt.plot(x, label='EKG')
plt.plot(peaks, x[peaks], "x", ms=10, label='peaks')
plt.legend(loc='best')
plt.show()