import marimo __generated_with = "0.24.2" app = marimo.App() @app.cell def _(): import marimo as mo return (mo,) @app.cell(hide_code=True) def _(mo): # Cell tags: task mo.md(r""" # Sheet 1: Error Propagation """) return @app.cell def _(): # Cell tags: task import matplotlib.pyplot as plt import numpy as np # create instance of pseudo-random number generator with fixed seed = 1 rng = np.random.default_rng(seed=1) return np, plt, rng @app.cell(hide_code=True) def _(mo): # Cell tags: task mo.md(r""" --- ## Task 1: Error Propagation Consider the function $y = f(x)$, $y = 1 + a_1x + a_2x^2$ with parameters $a_1 = 2.0 ± 0.2$, $a_2 = 1.0 ± 0.1$ and correlation coefficient $ρ = −0.8$. ### 1.1 Write down the covariance matrix for $a_1$ and $a_2$. --- """) return @app.cell(hide_code=True) def _(mo): mo.md(r""" The covariance matrix is $$ \mathrm{Cov}(a) = \begin{pmatrix}\sigma^2_{a_1} & \rho\sigma_{a_1}\sigma_{a_2} \\ \rho\sigma_{a_1}\sigma_{a_2} & \sigma^2_{a_2} \end{pmatrix} $$ """) return @app.cell def _(np): a1, a1_err = 2.0, 0.2 a2, a2_err = 1.0, 0.1 rho = -0.8 c12 = rho * a1_err * a2_err cov_a = np.array([[a1_err**2, c12], [c12, a2_err**2]]) cov_a return a1, a2, cov_a @app.cell(hide_code=True) def _(mo): # Cell tags: task mo.md(r""" --- ### 1.2 Compute the uncertainty of $y$ analytically using error propagation. --- """) return @app.cell(hide_code=True) def _(mo): mo.md(r""" We first compute derivatives of $y$ by $a_1$ and $a_2$ $$\frac{\partial{}y}{\partial{}a_1} = x \qquad \frac{\partial{}y}{\partial{}a_2} = x^2 \,.$$ We then apply the generic propagation formula using the Jacobi matrix and covariance matrix of inputs. Since the output of $f(x)$ is one-dimensional, the Jacobi matrix is reduced to a vector and the covariance matrix of $y$ is a scalar equal to the standard deviation squared, $$\sigma^2_y = \mathrm{Cov}(y) = \sum_{ij}\frac{\partial{}y}{\partial{}a_i}\frac{\partial{}y}{\partial{}a_j}\mathrm{Cov}(a)_{ij} \,.$$ Inserting the covariance matrix of $a_1$ and $a_2$ yields $$\sigma^2_y = c_{11}x^2 + 2c_{12}x^3 + c_{22}x^4 \,,$$ where the $c_{ij}$ are the entries of $\mathrm{Cov}(a)$. The expression can be simplified to: $$\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) \,.$$ Taking the square root: $$ \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} \,.$$ """) return @app.cell(hide_code=True) def _(mo): # Cell tags: task mo.md(r""" --- ### 1.3 Compare the numerical result for $\sigma_y$ for the two cases where you a) ignore the correlation $\rho$ by setting it to zero, and b) correctly include the correlation. Compute $\sigma_y$ for $x \in [-3, 3]$ and plot the results. In practice, correlations are often unknown. Is it conservative to set the correlation to zero? A conservative error estimate is always larger than the exact error. --- """) return @app.cell def _(cov_a, np): def compute_y(x, a1, a2): return 1 + a1 * x + a2 * x**2 def compute_err_y_no_correlation(x): return np.abs(x) * np.sqrt(cov_a[0, 0] + cov_a[1, 1] * x**2) def compute_err_y(x): return np.abs(x) * np.sqrt(cov_a[0, 0] + cov_a[1, 1] * x**2 + 2 * cov_a[0, 1] * x) return compute_err_y, compute_err_y_no_correlation, compute_y @app.cell def _(a1, a2, compute_err_y, compute_err_y_no_correlation, compute_y, np, plt): xs = np.linspace(-3, 3, 1000) _ys = compute_y(xs, a1, a2) # element-wise computation here, since we use numpy arrays errs_nc = compute_err_y_no_correlation(xs) errs = compute_err_y(xs) _fig, _ax = plt.subplots(1, 2, figsize=(10, 4), sharex=True) plt.sca(_ax[0]) plt.plot(xs, _ys) plt.fill_between(xs, _ys - errs_nc, _ys + errs_nc, alpha=0.5) plt.fill_between(xs, _ys - errs, _ys + errs, alpha=0.5) plt.sca(_ax[1]) plt.plot(xs, errs_nc / errs) plt.axhline(1, ls='--', color='0.5') plt.ylabel('error w/o correlation / full error') plt.show() return @app.cell(hide_code=True) def _(mo): mo.md(r""" Setting an unknown correlation to zero does not produce a conservative estimate of the error in general. We see in this particular example that setting the correlation to zero is only conservative for $x > 0$, but the true error on $y$ is underestimated for $x < 0$. """) return @app.cell(hide_code=True) def _(mo): # Cell tags: task mo.md(r""" --- ### 1.4 Compute the uncertainty of $y$ numerically using Monte-Carlo simulation under the assumption that $a_1$ and $a_2$ are normally distributed. Generate pairs $(a_{1,i}, a_{2,i})$ from the given values of $a_1$ and $a_2$ and their covariance matrix and visualize these pairs with a scatter plot or 2D histogram. Use `help(rng)` and make yourself familiar with the available methods to draw random numbers from a variety of statistical distributions. Locate a method that generates correlated normally distributed numbers directly. --- """) return @app.cell def _(a1, a2, cov_a, plt, rng): # help(rng) a1s, a2s = rng.multivariate_normal((a1, a2), cov_a, size=(100_000)).T plt.figure() plt.hist2d(a1s, a2s, range=((1.3, 2.7), (0.7, 1.3)), bins=40) plt.title("distribution of $a_1$ and $a_2$") plt.xlabel("$a_1$") plt.ylabel("$a_2$"); plt.show() return a1s, a2s @app.cell(hide_code=True) def _(mo): # Cell tags: task mo.md(r""" --- ### 1.5 Determine the distribution of $y_i$-samples for $x = \{-1, 0, +1\}$ and compare their mean and variance with the results of the analytical calculation. --- """) return @app.cell def _(a1s, a2s, compute_err_y, compute_y, np, plt): # we previously computed y and err_y in 1.2.2 for _x in (-1, 0, 1): _ys = compute_y(_x, a1s, a2s) mean = np.mean(_ys) # a1s, a2s are the Gaussian distributions from before var = np.var(_ys) plt.figure() plt.hist(_ys, bins=100) plt.xlabel(f'y({_x})') plt.ylabel('counts') plt.title(f'y({_x})= {mean:.3f}, variance = {var:.3f}, analytical = {compute_err_y(_x) ** 2:.3f}') plt.show() return @app.cell(hide_code=True) def _(mo): mo.md(r""" The case $x = 0$ is special here. Since all coefficients precede powers of $x$, the case always yields $y=0$ regardless of the $a_i$. """) return @app.cell(hide_code=True) def _(mo): # Cell tags: task mo.md(r""" --- ## Task 2: Error Propagation with Transformation Now consider the following reparametrisation of $y = f(x)$: $$y = 1 + \frac{x(1+x)}{b_1} + \frac{x(1-x)}{b_2}$$ ### 2.1 Determine analytically the transformed parameters $b_1$ and $b_2$ and their covariance matrices. --- """) return @app.cell(hide_code=True) def _(mo): mo.md(r""" The transformed covariance matrix is $\mathrm{Cov}(b) = J \,\mathrm{Cov}(a) \, J^T$, where $J$ is the Jacobi matrix containing the derivatives of the new parameters as a function of the old parameters. We first express the $b_i$ as function of the $a_i$. Here we can neglect the term $+1$ because it occurs equally in both definitions. \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} Thus $$a_1 = \left(\frac{1}{b_1} + \frac{1}{b_2}\right) \quad\text{and}\quad a_2 = \left(\frac{1}{b_1} - \frac{1}{b_2}\right)$$ and $$b_1 = \frac{2}{a_1 + a_2} \quad \text{and} \quad b_2 = \frac{2}{a_1 - a_2} \,.$$ For the Jacobian matrix of the transformation we get \begin{equation} J = \begin{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{pmatrix} \quad \text{where} \quad J_{ij} = \frac{\partial b_i}{\partial a_j} \, . \end{equation} """) return @app.cell def _(a1, a2, cov_a, np): def trafo(a): a1, a2 = a b1 = 2 / (a1 + a2) b2 = 2 / (a1 - a2) return np.array([b1, b2]) denom1 = (a1 + a2) ** 2 denom2 = (a1 - a2) ** 2 J = np.array([[-2 / denom1, -2 / denom1], [-2 / denom2, 2 / denom2]]) cov_b = J @ cov_a @ J.T trafo((a1, a2)), cov_b return cov_b, trafo @app.cell(hide_code=True) def _(mo): mo.md(r""" --- ### 2.2 Determine the covariance matrix of the transformed parameters by Monte-Carlo simulation and compare with the analytical calculation. _Hint_: If the covariance matrix is not what you expect, plot histograms of the transformed parameters and look for suspicious properties. --- """) return @app.cell def _(a1s, a2s, cov_b, np, trafo): # generate b-distributions from previously generated a-distributions b1s, b2s = trafo((a1s, a2s)) _ncov_b = np.cov(b1s, b2s) print(_ncov_b) print(_ncov_b / cov_b) return b1s, b2s @app.cell(hide_code=True) def _(mo): mo.md(r""" We observe large deviations of the simulation from the analytical solution, which seem to originate from the parameter $b_2$. Further investigation of the distribution reveals outliers in the $b_2$ sample. """) return @app.cell def _(b1s, b2s, np, plt): _fig, _ax = plt.subplots(1, 2, figsize=(10, 4)) _ax[0].hist(b1s, bins=60) _ax[0].set(title=f'mean = {np.mean(b1s):.3f}, var = {np.var(b1s):.3f}', xlabel='b1s') _ax[1].hist(b2s, bins=100) _ax[1].set(title=f'mean = {np.mean(b2s):.3f}, var = {np.var(b2s):.3f}', xlabel='b2s') plt.show() plt.figure() # plot y-axis logarithmically to reveal outliers plt.hist(b2s, bins=100) plt.xlabel('b2s') plt.semilogy() plt.show() return @app.cell(hide_code=True) def _(mo): mo.md(r""" The problem arises that for some combinations of values for $a_1$ and $a_2$ the denominator comes very close to `0`. This results in very large (unrealistic) values for $b_2$. We can counteract this by removing the smallest and largest 1% of the distribution. """) return @app.cell def _(b1s, b2s, np, plt): b2_low, b2_high = np.percentile(b2s, (1, 99)) print(f'b2_low: {b2_low:.3f},', f'b2_high: {b2_high:.3f}') mask = np.logical_and(b2_low < b2s, b2s < b2_high) b1s_cut = b1s[mask] # all values that meet this cut are written into a new array b2s_cut = b2s[mask] _fig, _ax = plt.subplots(1, 2, figsize=(10, 4)) _ax[0].hist(b1s_cut, bins=60) _ax[0].set(title=f'mean = {np.mean(b1s):.3f} var = {np.var(b1s):.3f}', xlabel='b1s(cut)') _ax[1].hist(b2s_cut, bins=100) _ax[1].set(title=f'mean = {np.mean(b2s):.3f} var = {np.var(b2s):.3f}', xlabel='b2s(cut)') print(f'Cut efficiency: {len(b1s_cut) / len(b1s)}') plt.show() return b1s_cut, b2s_cut @app.cell(hide_code=True) def _(mo): mo.md(r""" The covariance matrix of the cutted distribution is closer to the analytical result. They are not expected to be equal because the transformation is non-linear. """) return @app.cell def _(b1s_cut, b2s_cut, cov_b, np): _ncov_b = np.cov(b1s_cut, b2s_cut) (_ncov_b, _ncov_b / cov_b) return @app.cell(hide_code=True) def _(mo): mo.md(r""" --- ## Task 3: Bias Consider the transformation $y = f(x) = \sin(x)$ for a normally distributed $x$ with $\bar x = 0.8 \pm 0.3$. In this case, $f(\bar x)$ is not equal to $\bar y$, the mean of the $y$-distribution, since $f(x)$ is not linear. ### 3.1 Calculate the bias $f(\bar x) - \bar y$ using the approximate formula from the lecture and subtract the bias from the naive result $f(\bar x)$. --- """) return @app.cell(hide_code=True) def _(mo): mo.md(r""" The general formula to compute the bias for a function $y = f(x)$ to using a second order Taylor expansion is: $$f(\bar x) - \bar y \approx -\frac{1}{2}\frac{\partial^2\!f}{\partial x^2}\, \sigma^2_{x}$$ """) return @app.cell def _(np): _x = 0.8 err_x = 0.3 var_x = err_x ** 2 _y = np.sin(_x) d2f_dx2 = -np.sin(_x) bias = -0.5 * d2f_dx2 * var_x print(f'Naive function result: {_y:.3f}') print(f'Bias: {bias:.3f}') print(f'corrected function result: {_y - bias:.3f}') return @app.cell(hide_code=True) def _(mo): mo.md(r""" --- ## 3.2 Check your calculation by Monte-Carlo simulation, by comparing the corrected function result with the mean of the transformed sample. --- """) return @app.cell def _(np, plt, rng): _x = rng.normal(0.8, 0.3, size=100000) _y = np.sin(_x) plt.figure() plt.hist(_y, bins=100, label='Sample') plt.axvline(np.mean(_y), c='C1', label=f'$\\bar y = {np.mean(_y):.2f}$') plt.legend() plt.xlabel('y') plt.show() return @app.cell def _(): return if __name__ == "__main__": app.run()