# 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.