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+/-6
6+/-5

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()

1.34+/-0.10
2.00+/-0.10
-1.01+/-0.13

[[ -3.72741161e+12   1.88047118e+04   2.77899213e+12]
 [  1.88047014e+04   8.02145082e-03  -1.40199625e+04]
 [  2.77899213e+12  -1.40199704e+04  -2.07189279e+12]]
Out[9]:
<matplotlib.colorbar.Colorbar at 0x7fe283bc7ac8>

Vorsicht

Man kann keine ufloats 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}}$$