# Automatisieren von Fehlerrechnung¶

$$ \sigma_{f} = \sqrt{\sum_{i=1}^m \left( \frac{\partial f}{\partial x_i}\right)^{\!2} \sigma_{x_i}^2} $$

- Gesucht: Fehler von $f(x_i)$, wenn $x_i$ Fehler haben
- Manuelle Fehlerfortpflanzung:
  1. Berechne die Ableitungen von $f$ nach allen fehlerbehafteten Größen $x_i$
  2. Ableitungen in die obere Formel einsetzen
  3. Werte und Fehler der $x_i$ einsetzen
- Probleme:
  - Kompliziert, dauert lange, man macht oft Fehler
  - Falsches Ergebnis, wenn $x_i$ korreliert sind

$$ \operatorname{cov}(\vec{f}) = \boldsymbol{J} \operatorname{cov}(\vec{x})\, \boldsymbol{J}^\top $$

$\boldsymbol{J}$ ist die Jacobimatrix von $\vec{f}$ nach $\vec{x}$.

## `uncertainties`¶

- Erlaubt es, Fehlerrechnung automatisch durchzuführen
- Datentyp: `ufloat`, repräsentiert Wert mit Fehler

In [1]:

    from uncertainties import ufloat

    x = ufloat(5, 1)
    y = ufloat(3, 1)

    x + y

Out[1]:

    8.0+/-1.4142135623730951

Korrelationen werden von uncertainties beachtet:

In [2]:

    x = ufloat(3, 1)
    y = ufloat(3, 1)

    print(x**2 - y)
    print(x**2 - x) # error is smaller!

    6.0+/-6.1
    6.0+/-5.0

`uncertainties.unumpy` ergänzt numpy:

In [3]:

    import numpy as np
    import uncertainties.unumpy as unp

    x = [1, 2, 3, 4, 5]
    err = [0.1, 0.3, 0.1, 0.8, 1.0]

    y = unp.uarray(x, err)

    unp.cos(unp.exp(y))

Out[3]:

    array([-0.9117339147869651+/-0.11166193174450133,
           0.4483562418187328+/-1.9814233218473645,
           0.3285947554325321+/-1.8970207322669204,
           -0.3706617333977958+/-40.567208903209576,
           -0.7260031145123346+/-102.06245489729305], dtype=object)

Zugriff auf Wert und Standardabweichung mit `n` und `s`:

In [4]:

    x = ufloat(5, 1)
    print(x.n)
    print(x.s)

    5.0
    1.0

Bei `unumpy` mit `nominal_values` und `std_devs`

In [5]:

    x = unp.uarray([1, 2, 3], [0.3, 0.3, 0.1])
    print(unp.nominal_values(x))
    print(unp.std_devs(x))

    [ 1.  2.  3.]
    [ 0.3  0.3  0.1]

Kann man natürlich auch abkürzen:

In [6]:

    from uncertainties.unumpy import (nominal_values as noms,
                                      std_devs as stds)

    print(noms(x))
    print(stds(x))

    [ 1.  2.  3.]
    [ 0.3  0.3  0.1]

Man muss daran denken, die Funktionen aus unumpy zu benutzen (`exp`, `cos`, etc.)

In [7]:

    np.cos(x)

    ---------------------------------------------------------------------------
    AttributeError                            Traceback (most recent call last)
    <ipython-input-7-625fa4da891e> in <module>()
    ----> 1 np.cos(x)

    AttributeError: 'Variable' object has no attribute 'cos'

# Korrelierte Werte¶

In [8]:

    from uncertainties import correlated_values

    values = [1, 2]

    cov = [[0.5, 0.25],
           [0.25, 0.2]]

    x, y = correlated_values(values, cov)

### Vorsicht bei Fits:¶

korrelierte Fit-Parameter führen zu nichts-sagenden Ergebnissen. Kontrolle: Kovarianzmatrix.

In [9]:

    %matplotlib inline

    from numpy.random import normal
    import numpy as np
    import matplotlib.pyplot as plt
    plt.rcParams['figure.figsize'] = (10, 8)
    plt.rcParams['font.size'] = 16
    from scipy.optimize import curve_fit
    from uncertainties import correlated_values, correlation_matrix

    def f1(x, a, b):
        return a * np.exp(b * x)

    def f2(x, a, b, c):
        return a * np.exp(b * x + c)

    x = np.linspace(0, 2, 30)
    y = 0.5 * np.exp(2 * x) + normal(0, 1, 30)

    #params, cov = curve_fit(f1, x, y)
    params, cov = curve_fit(f2, x, y)

    params = correlated_values(params, cov)

    for param in params:
        print(param)
    print()
    print(cov)

    plt.matshow(correlation_matrix(params), cmap='seismic')
    plt.colorbar()

    (0.0+/-2.5)e+06
    1.89+/-0.09
    (-0.0+/-1.3)e+06

    [[  6.44166811e+12  -4.66200723e+04  -3.20060186e+12]
     [ -4.66200868e+04   6.63955723e-03   2.31636064e+04]
     [ -3.20060186e+12   2.31635995e+04   1.59024838e+12]]

Out[9]:

    <matplotlib.colorbar.Colorbar at 0x7efe2c130908>

## Vorsicht¶

Man kann keine `ufloat`s plotten:

In [10]:

    x = np.linspace(0, 10)
    y = unp.uarray(np.linspace(0, 5), 1)

    #plt.plot(x, y, 'rx')
    plt.errorbar(x, unp.nominal_values(y), yerr=unp.std_devs(y), fmt='rx')

Out[10]:

    <Container object of 3 artists>

# Sympy¶

- Kann Ableitungen automatisch generieren

SymPy importieren:

In [11]:

    import sympy

Mathematische Variablen erzeugen mit `var()`:

In [12]:

    x, y, z = sympy.var('x y z')

    x + y + z

Out[12]:

    x + y + z

Differenzieren mit `diff()`:

In [13]:

    f = x + y**3 - sympy.cos(z)**2

    print(f.diff(x))
    print(f.diff(y))
    print(f.diff(z))
    print(f.diff(z, z, z))

    1
    3*y**2
    2*sin(z)*cos(z)
    -8*sin(z)*cos(z)

Eine Funktion, die automatisch die Fehlerformel generiert:

In [14]:

    import sympy

    def error(f, err_vars=None):
        from sympy import Symbol, latex
        s = 0
        latex_names = dict()

        if err_vars == None:
            err_vars = f.free_symbols

        for v in err_vars:
            err = Symbol('latex_std_' + v.name)
            s += f.diff(v)**2 * err**2
            latex_names[err] = '\\sigma_{' + latex(v) + '}'

        return latex(sympy.sqrt(s), symbol_names=latex_names)

    E, q, r = sympy.var('E_x q r')

    f = E + q**2 * r

    print(f)
    print(error(f))
    print()

    E_x + q**2*r
    \sqrt{\sigma_{E_{x}}^{2} + 4 \sigma_{q}^{2} q^{2} r^{2} + \sigma_{r}^{2} q^{4}}

$$f= E + q^2 r \quad\rightarrow\quad \sigma_f = \sqrt{\sigma_{E_{x}}^{2} + 4 \sigma_{q}^{2} q^{2} r^{2} + \sigma_{r}^{2} q^{4}}$$