In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
plt.style.use('ggplot')
In [2]:
a1, a1_err = 2.0, 0.2
a2, a2_err = 1.0, 0.1
rho = -0.8
Die Kovarianzmatrix von $a_1$ und $a_2$ lässt sich wie folgt berechnen:
$$\mathrm{Cov}(a) = \pmatrix{\sigma^2_{a_1} & \rho\sigma_{a_1}\sigma_{a_2} \\ \rho\sigma_{a_1}\sigma_{a_2} & \sigma^2_{a_2}}$$
In [3]:
c12 = rho * a1_err * a2_err
covariance = np.matrix([[a1_err ** 2, c12], [c12, a2_err ** 2]])
covariance
Out[3]:
Dazu bestimmen wir zunaechst die Ableitungen von $y$ nach $a_1$ und $a_2$
$$\frac{\partial{}y}{\partial{}a_1} = x \,,\quad \frac{\partial{}y}{\partial{}a_2} = x^2 \,.$$Daraus können wir die Kovarianz von $y$ nach
$$\sigma^2_y = \mathrm{Cov}(y) = \sum_{ij}\frac{\partial{}y}{\partial{}a_i}\frac{\partial{}y}{\partial{}a_j}\mathrm{Cov}(a)_{ij}$$bestimmen. Einsetzen der eingangs bestimmten Ableitungen liefert
$$\sigma^2_y = c_{11}x^2 + 2c_{12}x^3 + c_{22}x^4$$wobei die $c_{ij}$ die Einträge von $\mathrm{Cov}(a)$ sind. Der Ausdruck lässt sich weiter vereinfachen zu
$$\sigma^2_y = x^2\left(\sigma^2_{a_1} + \sigma^2_{a_2}x^2 + 2\rho\sigma_{a_1}\sigma_{a_2}x\right) \,.$$Daraus ergibt sich für die Unsicherheit
$$ \sigma_y = \lvert{}x\rvert\sqrt{\sigma^2_{a_1} + \sigma^2_{a_2}x^2 + 2\rho\sigma_{a_1}\sigma_{a_2}x} \,.$$(also für $\rho = 0$) vereinfacht sich der obige Ausdruck zu
$$\sigma_y = \lvert{}x\rvert\sqrt{\sigma^2_{a_1} + \sigma^2_{a_2}x^2} \,.$$Mit den eingangs berechneten Werten aus covariance ergibt sich also
In [4]:
def err_ana_wo(x):
return 0.2 * np.abs(x) * np.sqrt(1 + 0.25 * x ** 2)
xs = np.linspace(-3, 3, 10000)
ys = 1 + a1 * xs + a2 * xs ** 2
errs = err_ana_wo(xs)
plt.plot(xs, ys)
plt.fill_between(xs, ys - errs, ys + errs, alpha=0.5)
plt.show()
In [5]:
def err_ana(x):
return 0.2 * np.abs(x) * np.sqrt(1 + 0.25 * x ** 2 - 0.4 * x)
errs = err_ana(xs)
plt.plot(xs, ys)
plt.fill_between(xs, ys - errs, ys + errs, alpha=0.5)
plt.show()
Generieren Sie Wertepaare $(a_1, a_2)$ gemäß ihrer Kovarianzmatrix und visualisieren Sie diese, z.B. mit einem Scatter-Plot.
Hinweis: Wenn $x_1$ und $x_2$ zwei gaussverteilte Zufallszahlen mit Mittelwert null und Varianz eins sind, erhält man ein Paar korrelierter gaussverteilter Zufallszahlen $(y_1, y_2)$ mit Mittelwert null und Varianz eins durch $(y_1 = x_1; y_2 = x_1ρ + x_2\sqrt{1 − \rho^2})$.
In [6]:
x1s, x2s = np.random.normal(size=(2, 10000))
plt.hist2d(x1s, x2s, bins=40)
plt.title('2-dim Normalverteilung')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.show()
y1s = x1s
y2s = x1s * rho + x2s * np.sqrt(1 - rho ** 2)
plt.hist2d(y1s, y2s, bins=40)
plt.title('2-dim Normalverteilung mit Korrelation')
plt.xlabel('$y_1$')
plt.ylabel('$y_2$')
plt.show()
a1s = a1 + y1s * a1_err
a2s = a2 + y2s * a2_err
plt.hist2d(a1s, a2s, bins=40)
plt.title('Verteilung der $a_1$ und $a_2$')
plt.xlabel('$a_1$')
plt.ylabel('$a_2$')
plt.show()
In [7]:
def y(x, a1, a2):
return 1 + a1 * x + a2 * x ** 2
def var_analytical(x):
xx = x ** 2
#return xx * (covariance[0, 0] + covariance[1, 1] * xx + 2 * covariance[0, 1] * x)
return 0.04 * xx * (1 + 0.25 * xx - 0.8 * x)
for x in (-1, 0, 1):
ys = y(x, a1s, a2s)
mean = np.mean(ys)
var = np.var(ys)
print('〈y({})〉= {:.3f}'.format(x, mean))
print(' σ² = {:.3f}'.format(var))
print(' analytical = {:.3f}'.format(var_analytical(x)))
plt.hist(ys, bins=100)
plt.xlabel('y({})'.format(x))
plt.ylabel('Absolute Häufigkeit')
plt.show()
Der Fall $x = 0$ ist hier besonders. Da alle Koeffizienten vor Potenzen von $x$ stehen, ergibt sich für den Fall immer $y=0$ unabhängig von den $a_i$. Wir können also keine Aussage über die Varianz treffen.
Wir lösen die Reparametrisierung nach Koeffizienten von Potenzen von $x$ auf. Dabei können wir den Term $1$ vernachlässigen, weil er in beiden Definitionen gleichermaßen auftritt.
\begin{align} a_1 x + a_2 x^2 &= \frac{x(1 + x)}{b_1} + \frac{x(1 - x)}{b_2} \\ &= \frac{x}{b_1} + \frac{x^2}{b_1} + \frac{x}{b_2} - \frac{x^2}{b_2} \\ &= x\left(\frac{1}{b_1} + \frac{1}{b_2}\right) + x^2\left(\frac{1}{b_1} - \frac{1}{b_2}\right) \end{align}Damit ist
$$a_1 = \left(\frac{1}{b_1} + \frac{1}{b_2}\right) \quad\text{und}\quad a_2 = \left(\frac{1}{b_1} - \frac{1}{b_2}\right)$$also
$$b_1 = \frac{2}{a_1 + a_2} \quad\text{und}\quad b_2 = \frac{2}{a_1 - a_2} \,.$$Für die Jacobimatrix der Transformation ergibt sich
\begin{equation} M = \pmatrix{ \frac{-2}{(a_1 + a_2)^2} & \frac{-2}{(a_1 + a_2)^2} \\ \frac{-2}{(a_1 - a_2)^2} & \frac{+2}{(a_1 - a_2)^2} } \quad\text{wobei}\quad m_{ij} = \frac{\partial b_i}{\partial a_j} \end{equation}\begin{equation} M^T = \pmatrix{ \frac{-2}{(a_1 + a_2)^2} & \frac{-2}{(a_1 - a_2)^2} \\ \frac{-2}{(a_1 + a_2)^2} & \frac{+2}{(a_1 - a_2)^2} } \end{equation}Die transformierte Kovarianzmatrix ist dann $\mathrm{Cov}(b) = M\mathrm{Cov}(a)M^T$.
In [8]:
b1 = 2 / (a1 + a2)
b2 = 2 / (a1 - a2)
denom1 = (a1 + a2) ** 2
denom2 = (a1 - a2) ** 2
M = np.matrix([[-2 / denom1, -2 / denom1],
[-2 / denom2, 2 / denom2]])
cov_b = M * covariance * M.T
cov_b
Out[8]:
In [9]:
b1s = 2 / (a1s + a2s)
b2s = 2 / (a1s - a2s)
print('b1 = {}'.format(np.mean(b1s)))
print('var = {}'.format(np.var(b1s)))
print('b2 = {}'.format(np.mean(b2s)))
print('var = {}'.format(np.var(b2s)))
plt.hist(b1s, bins=100)
plt.show()
plt.hist(b2s, bins=100)
plt.show()
Dabei tritt das Problem auf, dass für einige Kombinationen von Werten für $a_1$ und $a_2$ der Nenner sehr nah an 0 kommt. Dadurch ergeben sich sehr große (unrealistische) Werte für $b_2$. Wir können dem entgegenwirken, indem wir einen Bereich festlegen in dem wir die Werte für $b_2$ erwarten.
In [10]:
cut = np.logical_and(b2s < 5, b2s > 0)
b1s_ = b1s[cut]
b2s_ = b2s[cut]
print('b2 gefiltert = {}'.format(np.mean(b2s_)))
print('var = {}'.format(np.var(b2s_)))
plt.hist(b1s_, bins=100)
plt.show()
plt.hist2d(b1s_, b2s_, bins=100)
plt.show()
Wenn wir uns das Leben etwas erleichtern wollen, koennen wir auch einfach die Funktion cov aus Numpy verwenden, die uns fuer zwei Arrays direkt die Kovarianzmatrix ausrechnet
In [11]:
ncov_b = np.cov(b1s_, b2s_)
ncov_b
Out[11]:
Zunächst berechnen wir die partiellen Ableitungen von $y$ nach den Koeffizienten $b_1$ und $b_2$ und damit die Jacobimatrix $M$.
\begin{equation} M = \pmatrix{\frac{\partial y}{\partial b_1} \\ \frac{\partial y}{\partial b_2}} = \pmatrix{\frac{-x(1+x)}{b_1^2}\\ \frac{-x(1-x)}{b_2^2}} \end{equation}Damit ergibt sich für die Varianz
\begin{align} \sigma_y^2 &= M^T \mathrm{Cov}(b) M \\ &= x^2\left[\frac{c_{11}}{b_1^4}(1+x)^2 + \frac{2c_{12}}{b_1^2b_2^2}(1 - x^2) + \frac{c_{22}}{b_2^4}(1-x)^2\right] \\ &= x^2\left(\alpha + \beta x + \gamma x^2 \right) \end{align}mit den Koeffizienten
\begin{equation} \alpha = \left(\frac{c_{11}}{b_1^4} + \frac{2c_{12}}{b_1^2b_2^2} + \frac{c_{22}}{b_2^4}\right) \quad,\quad \beta = 2\left(\frac{c_{11}}{b_1^4} - \frac{c_{22}}{b_2^4}\right) \quad\text{und}\quad \gamma = \left(\frac{c_{11}}{b_1^4} - \frac{2c_{12}}{b_1^2b_2^2} + \frac{c_{22}}{b_2^4}\right) \,. \end{equation}In Zahlen ausgedrückt sind die Koeffizienten
In [12]:
c11 = cov_b[0, 0]
c12 = cov_b[0, 1]
c22 = cov_b[1, 1]
print(cov_b)
a = c11 / b1 ** 4
b = 2 * c12 / b1 ** 2 / b2 ** 2
c = c22 / b2 ** 4
alpha = a + b + c
beta = 2 * (a - c)
gamma = a - b + c
np.sqrt(alpha), alpha/alpha, beta/alpha, gamma/alpha
Out[12]:
Es ist also
\begin{equation} \sigma_y = 0.2\lvert{}x\rvert\sqrt{1 - 0.8x - 0.25x^2} \end{equation}was exakt das gleiche Ergbnis ist wie für die ursprüngliche Parametrisierung.
In [13]:
def err_ana(x):
return np.abs(x) * np.sqrt(alpha + beta * x + gamma * x ** 2)
xs = np.linspace(-10, 10, 10000)
ys = 1 + xs * (xs + 1) / b1 + xs * (xs - 1) / b2
errs = err_ana(xs)
plt.plot(xs, ys)
plt.fill_between(xs, ys - errs, ys + errs, alpha=0.5)
plt.show()
Die gleiche Rechnung mit den per Monte Carlo bestimmten Werten
In [14]:
c11_n = ncov_b[0, 0]
c12_n = ncov_b[0, 1]
c22_n = ncov_b[1, 1]
a_n = c11_n / b1 ** 4
b_n = 2 * c12_n / b1 ** 2 / b2 ** 2
c_n = c22_n / b2 ** 4
alpha_n = a_n + b_n + c_n
beta_n = 2 * (a_n - c_n)
gamma_n = a_n - b_n + c_n
np.sqrt(alpha_n), alpha_n/alpha_n, beta_n/alpha_n, gamma_n/alpha_n
Out[14]:
In [15]:
def err_num(x):
return np.abs(x) * np.sqrt(alpha_n + beta_n * x + gamma_n * x ** 2)
errs = err_num(xs)
plt.plot(xs, ys)
plt.fill_between(xs, ys - errs, ys + errs, alpha=0.5)
plt.show()
Im folgenden nennen wir das neue $y$ der Einfachheit halber $z$.
\begin{equation} z = \ln(1 + a_1x + a_2x^2) = \ln(y) \end{equation}Für die Unsicherheit von $z$ ergibt sich
\begin{align} \sigma_z &= \sqrt{\left(\frac{\partial z}{\partial y}\right)^2\sigma_y^2} \\ &= \sqrt{\left(\frac{1}{y}\right)^2\sigma_y^2} \\ &= \left\lvert\frac{\sigma_y}{y}\right\rvert \\ &= \frac{0.2\lvert x\rvert\sqrt{1 + 0.25x^2 - 0.4x}}{\lvert 1 + a_1 x + a_2 x^2 \rvert} \,. \end{align}Völlig analog ist die Rechnung für die Reparametrisierung. Hier ergibt sich
\begin{align} \sigma_z &= \left\lvert\frac{\sigma_y}{y}\right\rvert \\ &= \frac{0.2\lvert x\rvert\sqrt{1 + 0.25x^2 - 0.4x}}{\left\lvert 1 + \frac{x(1+x)}{b_1} + \frac{x(x-1)}{b_2}\right\rvert} \,. \end{align}
In [16]:
def err1(x):
return err_ana(x) / np.abs(1 + a1 * x + a2 * x ** 2)
def err2(x):
return err_ana(x) / np.abs(1 + xs * (xs + 1) / b1 + xs * (xs - 1) / b2)
xs = np.linspace(-3, 3, 100)
ys1 = 1 + a1 * xs + a2 * xs ** 2
ys2 = np.log(1 + xs * (xs + 1) / b1 + xs * (xs - 1) / b2)
errs1 = err1(xs)
errs2 = err2(xs)
plt.plot(xs, ys1)
plt.fill_between(xs, ys1 - errs1, ys1 + errs1, alpha=0.5)
plt.show()
plt.plot(xs, ys2)
plt.fill_between(xs, ys2 - errs2, ys2 + errs2, alpha=0.5)
plt.show()