<jupyter_start><jupyter_text>Musterlösung zu Übungsblatt 1  * [Aufgabe 1](Aufgabe-1) * [Aufgabe 2](Aufgabe-2) * [Aufgabe 3](Aufgabe-3)<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
plt.style.use('ggplot')<jupyter_output><empty_output><jupyter_text>--- Aufgabe 1Gegeben sei eine parametrische Funktion $y = f(x)$, $y = 1 + a_1x + a_2x^2$ mit Parametern $a_1 = 2.0 ± 0.2$, $a_2 = 1.0 ± 0.1$ und Korrelationskoeffizient $ρ = −0.8$.---<jupyter_code>a1, a1_err = 2.0, 0.2
a2, a2_err = 1.0, 0.1
rho = -0.8<jupyter_output><empty_output><jupyter_text>--- 1.1Geben Sie die Kovarianzmatrix von $a_1$ und $a_2$ an.--- 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}}$$<jupyter_code>c12 = rho * a1_err * a2_err
covariance = np.matrix([[a1_err ** 2, c12], [c12, a2_err ** 2]])
covariance<jupyter_output><empty_output><jupyter_text>--- 1.2Bestimmen Sie analytisch die Unsicherheit von $y$ als Funktion von $x$:--- 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} \,.$$ --- 1.2.1unter Vernachlässigung der Korrelation--- (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$$\sigma_y = \lvert{}x\rvert\sqrt{0.04 + 0.01x^2} = 0.2\lvert{}x\rvert{}\sqrt{1 + 0.25x^2} \,.$$<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>--- 1.2.2mit Berücksichtigung der Korrelation--- $$\sigma_y = \lvert{}x\rvert\sqrt{0.04 + 0.01x^2 - 0.016x} = 0.2\lvert{}x\rvert\sqrt{1 - 0.8x + 0.25x^2} \,.$$<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>--- 1.3Bestimmen Sie per Monte Carlo die Unsicherheit von $y$ als Funktion von $x$:  1.3.1Generieren Sie Wertepaare $(a_1, a_2)$ gemäß 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, erhä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})$.---<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>--- 1.3.2Bestimmen Sie die Verteilung von $y$ für $x = \{−1, 0, +1\}$ und vergleichen Sie Mittelwert und Varianz (Standardabweichung) mit den Resultaten der analytischen Rechnung.---<jupyter_code>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()<jupyter_output>〈y(-1)〉= 0.002
 σ² = 0.081
 analytical = 0.082<jupyter_text>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. --- Aufgabe 2Betrachten Sie folgende Reparametrisierung von $y = f(x)$:    $$y = 1 + \frac{x(1+x)}{b_1} + \frac{x(1-x)}{b_2}$$ 2.1Bestimmen Sie analytisch die transformierten Parameter $b_1$ und $b_2$ und deren Kovarianzmatrix--- 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$.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>--- 2.2Bestimmen Sie die Kovarianzmatrix der transformierten Parameter per Monte Carlo---<jupyter_code>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()<jupyter_output>b1 = 0.6681390380773017
var = 0.0009033648632514698
b2 = 2.2516701252806204
var = 1.797015387269852<jupyter_text>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.<jupyter_code>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()<jupyter_output>b2 gefiltert = 2.1462016045455883
var = 0.4459879545199094<jupyter_text>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<jupyter_code>ncov_b = np.cov(b1s_, b2s_)
ncov_b<jupyter_output><empty_output><jupyter_text>--- 2.3Bestimmen Sie analytisch die Unsicherheit von $y$ als Funktion von $x$: 2.3.1unter Verwendung der analytisch bestimmten Kovarianzmatrix von $(b_1, b_2)$--- 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<jupyter_code>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<jupyter_output>[[ 0.00088889  0.01333333]
 [ 0.01333333  0.328     ]]<jupyter_text>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.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>--- 2.3.2unter Verwendung der numerisch bestimmten Kovarianzmatrix von $(b_1, b_2)$--- Die gleiche Rechnung mit den per Monte Carlo bestimmten Werten<jupyter_code>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
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()<jupyter_output><empty_output><jupyter_text>--- Aufgabe 3Lösen Sie die obigen Teilaufgaben für $y = f(x)$ mit$$y = \ln\left(1 + a_1x + a_2x^2\right) \quad \text{bzw.} \quad y = \ln\left(1 + \frac{x(1+x)}{b_1} + \frac{x(x-1)}{b_2}\right)$$ 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}<jupyter_code>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()<jupyter_output>/opt/local/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in true_divide
  from ipykernel import kernelapp as app
{'fraction_non_alphanumeric': 0.19241022424337462, 'fraction_numerical': 0.06766504625999686, 'mean_word_length': 3.791510142749812, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '260250', 'license': 'MIT', 'n_tokens_mistral': 6154, 'n_tokens_neox': 5746, 'n_words': 1203, 'path': 'solutions/1.ipynb', 'repo_name': 'kdungs/teaching-SMD2-2016', 'size': '12754'}