SciPy

SciPy

  • Baut auf NumPy auf
  • Kann numerisch integrieren, DGLs lösen, optimieren, minimieren, …
  • Enthält auch physikalische Konstanten und wichtige mathematische Funktionen
In [1]:
import numpy as np
x = [2,5,7]
In [2]:
# standard error of the mean
from scipy.stats import sem
sem(x)
Out[2]:
1.4529663145135578
In [3]:
# physical constants
import scipy.constants as const
const.epsilon_0
Out[3]:
8.854187817620389e-12
In [4]:
# convert temperatures:
print(const.convert_temperature(100,'c','K'))
print(const.convert_temperature(100,'kelvin','Celsius'))

373.15
-173.15
In [5]:
# convert angles:
print(np.rad2deg(np.pi))
print(np.deg2rad(90))

180.0
1.57079632679
In [6]:
# more constants (including units and errors)!
list(const.physical_constants.items())[:10]
Out[6]:
[('Wien displacement law constant', (0.0028977685, 'm K', 5.1e-09)),
 ('atomic unit of 1st hyperpolarizablity',
  (3.20636151e-53, 'C^3 m^3 J^-2', 2.8e-60)),
 ('atomic unit of 2nd hyperpolarizablity',
  (6.2353808e-65, 'C^4 m^4 J^-3', 1.1e-71)),
 ('atomic unit of electric dipole moment', (8.47835309e-30, 'C m', 7.3e-37)),
 ('atomic unit of electric polarizablity',
  (1.648777274e-41, 'C^2 m^2 J^-1', 1.6e-49)),
 ('atomic unit of electric quadrupole moment',
  (4.48655124e-40, 'C m^2', 3.9e-47)),
 ('atomic unit of magn. dipole moment', (1.8548019e-23, 'J T^-1', 1.6e-30)),
 ('atomic unit of magn. flux density', (235051.755, 'T', 0.0014)),
 ('deuteron magn. moment', (4.33073482e-27, 'J T^-1', 3.8e-34)),
 ('deuteron magn. moment to Bohr magneton ratio',
  (0.0004669754567, '', 5e-12))]

Achtung

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)

In [7]:
const.physical_constants["proton mass"]
# value, unit, error
Out[7]:
(1.672621898e-27, 'kg', 2.1e-35)

Fitten

Oft möchte man eine Funktion, zum Beispiel eine Erwartung aus der Theorie, an die gemessenen Werte anpassen. Dies nennt man Fit.

In [8]:
# 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])

a = 4.97742635551 ± 0.0214390034102
b = 0.070084504726 ± 0.123809150457
In [9]:
# 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")
Out[9]:
<matplotlib.legend.Legend at 0x7f8a4964ceb8>

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)

$$ f(x; a, b, c) = \frac{a}{1 + \exp(-(x-b))} + c$$
In [10]:
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()
Out[10]:
<matplotlib.legend.Legend at 0x7f8a4957a3c8>

Die Messwerte aus einem Praktikumsversuch:

In [11]:
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$')
Out[11]:
<matplotlib.text.Text at 0x7f8a4ce279b0>

Ein einfacher Fit wie oben funktioniert hier nicht so gut:

In [12]:
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])

a = -307.055916178 ± inf
b = 1.0 ± inf

/home/maxnoe/.local/anaconda3/lib/python3.6/site-packages/scipy/optimize/minpack.py:779: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)

Schaut man sich die berechnete Ausgleichskurve an sieht man auch,
dass das nicht stimmen kann:

In [13]:
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')
Out[13]:
<matplotlib.legend.Legend at 0x7f8a491720b8>

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:

In [14]:
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])

a = -0.493935581363 ± 0.0138538297848
b = 40.6675535294 ± 0.136960842835
In [15]:
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')
Out[15]:
<matplotlib.legend.Legend at 0x7f8a490bb358>

Zum Vergleich der beiden Anfangswerte (seeds) kann man sich die einmal ansehen
und mit den angepassten Parametern vergleichen:

In [16]:
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()
Out[16]:
<matplotlib.legend.Legend at 0x7f8a49072c50>

Die richtigen Startwerte findet man entweder durch

  1. trial and error => einfach ausprobieren bis es klappt

  2. nachdenken => siehe unten

Im obigen Beispiel musste nur 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$$
In [17]:
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()
Out[17]:
<matplotlib.legend.Legend at 0x7f8a49072160>

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.