import numpy as np
x = [2,5,7]
# standard error of the mean
from scipy.stats import sem
sem(x)
# physical constants
import scipy.constants as const
const.epsilon_0
# 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))
# more constants (including units and errors)!
list(const.physical_constants.items())[:10]
Wenn solche Konstanten genutzt werden, muss das korrekt mitgeteilt, also zitiert werden. Darauf gehen wir nächste Woche im LaTeX-Workshop ein :-)
(Quelle hier: python-scipy)
const.physical_constants["proton mass"]
# value, unit, error
Oft möchte man eine Funktion, zum Beispiel eine Erwartung aus der Theorie, an die gemessenen Werte anpassen. Dies nennt man Fit.
# fit arbitrary functions
from scipy.optimize import curve_fit
x, y = np.genfromtxt('example_data.txt', unpack=True)
def f(x, a, b):
return a * x + b
params, covariance_matrix = curve_fit(f, x, y)
errors = np.sqrt(np.diag(covariance_matrix))
print('a =', params[0], '±', errors[0])
print('b =', params[1], '±', errors[1])
# prepare plot
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10, 8)
plt.rcParams['font.size'] = 16
x_plot = np.linspace(0, 10)
plt.plot(x, y, 'k.', label="example data")
plt.plot(x_plot, f(x_plot, *params), 'r-', label='linearer Fit', linewidth=3)
plt.legend(loc="best")
Da ein Fit mit curve_fit()
intern ein Minimierungsalgorithmus ist, hängt
das Ergebnis unter Umständen von den Anfangsbedingungen ab.
Fit einer komplexeren Funktion: Sigmoidfunktion
(Ähnlich zum tanh
)
def sigmoid(x, a, b, c):
y = a / (1 + np.exp(-(x-b))) + c
return y
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('fit_data_with_init_values.txt', unpack=True)
plt.plot(x, y, 'ro', label=r'Messwerte')
plt.xlabel('Temperatur / °C')
plt.ylabel('$GP$')
Ein einfacher Fit wie oben funktioniert hier nicht so gut:
params, covariance_matrix = curve_fit(sigmoid, x, y)
errors = np.sqrt(np.diag(covariance_matrix))
print('a =', params[0], '±', errors[0])
print('b =', params[1], '±', errors[1])
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, 'ro', label='Messdaten')
plt.plot(x_plot, sigmoid(x_plot, *params), "b-", 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))
errors = np.sqrt(np.diag(covariance_matrix))
print('a =', params[0], '±', errors[0])
print('b =', params[1], '±', errors[1])
plt.xlabel('Temperatur / °C')
plt.ylabel('$GP$')
x_plot = np.linspace(0, 50, 1000)
plt.plot(x, y, 'ro', label='Messwerte')
plt.plot(x_plot, sigmoid(x_plot, *params), "b-", 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 p in parameter:
plt.plot(x_plot, sigmoid(x_plot, *p), label="f(x; {0:0.3f}, {1:0.3f}, {2:0.3f})".format(*p))
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.
B = [0, 20, 40]
x_plot = np.linspace(-50, 50, 1000)
plt.xlabel('x')
plt.ylabel('y')
for b in B:
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[0].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.