# Investigating propagation of errors

This notebook was made to investigate the propagation of errors formula.
We imagine that we have a function $q(x,y)$ and we want to propagate the
uncertainty on $x$ and $y$ (denoted $\sigma_x$ and $\sigma_y$, respectively) through to the quantity $q$.

The most straight forward way to do this is just randomly sample $x$ and $y$, evaluate $q$ and look at it's distribution. This is really the definition of what we mean by propagation of uncertianty. It's very easy to do with some simply python code.

The calculus formula for the propagation of errors is really an approximation. This is the formula for a general $q(x,y)$
\begin{equation}
\sigma_q^2 = \left( \frac{\partial q}{\partial x} \sigma_x \right)^2 + \left( \frac{\partial q}{\partial y}\sigma_y \right)^2
\end{equation}

In the special case of addition  $q(x,y) = x\pm y$ we have $\sigma_q^2 = \sigma_x^2 + \sigma_y^2$.

In the special case of multiplication $q(x,y) = x y$ and division $q(x,y) = x / y$ we have $(\sigma_q/q)^2 = (\sigma_x/x)^2 + (\sigma_y/y)^2$, which we can rewrite as $\sigma_q = (x/y) \sqrt{(\sigma_x/x)^2 + (\sigma_y/y)^2}$

Let's try out these formulas and compare the direct approach of making the distribution to the prediction from these formulas


```python
#%pylab inline --no-import-all
```


```python
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as scs
```


```python
from ipywidgets import widgets  
from ipywidgets import interact, interactive, fixed
```

## Setup repeated observations of two variables x, y


```python
mean_x = .8
std_x = .15
mean_y = 3.
std_y = .9
N = 1000
x = np.random.normal(mean_x, std_x, N)
y = np.random.normal(mean_y, std_y, N)
```

## Check propagation of errors under addition 


```python
x = np.random.normal(mean_x, std_x, N)
y = np.random.normal(mean_y, std_y, N)

q_of_x_y = x+y

pred_mean_q = mean_x+mean_y
pred_std_q = np.sqrt(std_x**2+std_y**2)

counts, bins, patches = plt.hist(q_of_x_y, 
                                 bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), 
                                 density=True, alpha=0.3)

plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)
plt.legend(('pred','hist'))
plt.xlabel('q(x)')
plt.ylabel('p(x)')
```

### same thing with an interactive widget


```python
def plot_x_plus_y(mean_x, std_x, mean_y, std_y, N):
    x = np.random.normal(mean_x, std_x, N)
    y = np.random.normal(mean_y, std_y, N)

    q_of_x_y = x+y

    pred_mean_q = mean_x+mean_y
    pred_std_q = np.sqrt(std_x**2+std_y**2)

    counts, bins, patches = plt.hist(q_of_x_y, 
                                     bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), 
                                     density=True, alpha=0.3)

    plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)
    plt.legend(('pred','hist'))
    plt.xlabel('q(x)')
    plt.ylabel('p(x)')


# now make the interactive widget
interact(plot_x_plus_y,
         mean_x=(0.,3.,.1), std_x=(.0, 2., .1), 
         mean_y=(0.,3.,.1), std_y=(.0, 2., .1),
         N=(0,10000,1000))
```


    interactive(children=(FloatSlider(value=1.5, description='mean_x', max=3.0), FloatSlider(value=1.0, descriptio…





    <function __main__.plot_x_plus_y(mean_x, std_x, mean_y, std_y, N)>



## Single variable example: division

As a warm up, let's consider $q(x) = 1/x$


```python
counts, bins, patches = plt.hist(x, bins=50, density=True, alpha=0.3)
gaus_x = scs.norm.pdf(bins, mean_x,std_x)

q_for_plot = 1./bins

plt.plot(bins, gaus_x, lw=2)
plt.plot(bins, q_for_plot, lw=2)

plt.xlabel('x')
plt.ylabel('q(x)')
```


```python
plt.xlabel('x')

q_of_x = 1./x

pred_mean_q = 1./mean_x
pred_std_q = np.sqrt((std_x/mean_x)**2)/mean_x

counts, bins, patches = plt.hist(q_of_x, bins=50, density=True, alpha=0.3)

plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)
plt.legend(('pred','hist'))
plt.xlabel('x')
plt.ylabel('p(x)')
```

### Now let's do the same thing with an interactive widget!


```python
def plot_1_over_x(mean_x, std_x, N):
    x = np.random.normal(mean_x, std_x, N)

    q_of_x = 1./x

    pred_mean_q = 1./mean_x
    pred_std_q = np.sqrt((std_x/mean_x)**2)/mean_x

    counts, bins, patches = plt.hist(q_of_x, 
                                     bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), 
                                     density=True, alpha=0.3)

    plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)
    plt.legend(('pred','hist'))
    plt.xlabel('q(x)')
    plt.ylabel('p(x)')
```


```python
# now make the interactive widget
interact(plot_1_over_x,mean_x=(0.,3.,.1), std_x=(.0, 2., .1), N=(0,10000,1000))
```


    interactive(children=(FloatSlider(value=1.5, description='mean_x', max=3.0), FloatSlider(value=1.0, descriptio…





    <function __main__.plot_1_over_x(mean_x, std_x, N)>



### Check propagation of errors under division 


```python
def plot_x_over_y(mean_x, std_x, mean_y, std_y, N):
    x = np.random.normal(mean_x, std_x, N)
    y = np.random.normal(mean_y, std_y, N)

    q_of_x_y = x/y

    pred_mean_q = mean_x/mean_y
    pred_std_q = np.sqrt((std_x/mean_x)**2+(std_y/mean_y)**2)*mean_x/mean_y

    counts, bins, patches = plt.hist(q_of_x_y, 
                                     bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), 
                                     density=True, alpha=0.3)


    plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)
    plt.legend(('pred','hist'))
    plt.xlabel('q(x)')
    plt.ylabel('p(x)')



interact(plot_x_over_y,mean_x=(0.,3.,.1), std_x=(.0, 2., .1), mean_y=(0.,3.,.1), std_y=(.0, 2., .1),N=(0,100000,1000))
```


    interactive(children=(FloatSlider(value=1.5, description='mean_x', max=3.0), FloatSlider(value=1.0, descriptio…





    <function __main__.plot_x_over_y(mean_x, std_x, mean_y, std_y, N)>




```python

```


```python

```


```python

```


```python

```

d5e65bc56c9558005b79e11ebe6357a9f81ad720
ipynb
Jupyter Notebook
book/error-propagation/investigating-propagation-of-errors.ipynb
willettk/stats-ds-book
06bc751a7e82f73f9d7419f32fe5882ec5742f2f
MIT
2020-08-18T12:14:43.000Z
2022-03-31T16:37:17.000Z
book/error-propagation/investigating-propagation-of-errors.ipynb
willettk/stats-ds-book
06bc751a7e82f73f9d7419f32fe5882ec5742f2f
MIT
2020-08-19T04:22:24.000Z
2020-12-22T15:18:24.000Z
book/error-propagation/investigating-propagation-of-errors.ipynb
willettk/stats-ds-book
06bc751a7e82f73f9d7419f32fe5882ec5742f2f
MIT
2020-08-19T02:57:47.000Z
2022-03-03T15:24:07.000Z
Qwen/Qwen-72B
1. YES
2. YES