# Error Propagation using Random Gaussian Numbers

Example calculation of propagating uncertainties, both when adding and multiplying number, and also in the general case. The propagation can be done both analytically (using the error propagation formula) and also using simulation.

The example is based on FIRST doing the error propagation **analytically**, and then verifying it by running a so-called Monte-Carlo (MC) program, which uses random numbers for propagating errors.

## References:
- Barlow: page 48-61
- Bevington: page 36-48

## Author(s), contact(s), and dates:
- Author: Troels C. Petersen (NBI)
- Email:  petersen@nbi.dk
- Date:   25th of November 2021

***

DO THE FOLLOWING ANALYTICAL EXERCISE FIRST!!!

1. A class of students estimate by eye, that the length of the table in Auditorium A is $L = (3.5\pm 0.4)$m, and that the width is $W = (0.8\pm 0.2)$m.

   Assuming that there is no correlation between these two measurements, calculate ANALYTICALLY what the Perimeter (P), area (A), and diagonal (D) length is including (propagated) uncertainties. Repeat the calculation, given that the correlation between length and width is $\rho(L,W) = 0.5$ - not an unreasonable number, given that they are estimated by the same (uncertain) scale.
   
NOTE: This is a complete standard problem, that you will be asked to solve again and again in the course. For this reason, make sure that you understand how to do it, and become good at doing it reasonably fast.


```python
# Defining the parameters:
mu1   =  3.5
sig1  =  0.4
mu2   =  0.8
sig2  =  0.2
rho12 =  0.5           # Correlation parameter!
```


```python
if not (-1.0 <= rho12 <= 1.0): 
    raise ValueError(f"Correlation factor not in interval [-1,1], as it is {rho12:6.2f}")
```

### Note on analytic solutions:

Python includes symbolic algebra in the package *SymPy*, which is both simple and powerful (and in Python). In addition, printing with Latex can also be included (see below), which (in combination) is very nice.

Below is a SymPy and Latex example with the hope that it will wet your appetite.


```python
from IPython.core.display import Latex

def lprint(*args,**kwargs):
    """Pretty print arguments as LaTeX using IPython display system 
    
    Parameters
    ----------
    args : tuple 
        What to print (in LaTeX math mode)
    kwargs : dict 
        optional keywords to pass to `display` 
    """
    display(Latex('$$'+' '.join(args)+'$$'),**kwargs)
    
from IPython.core.display import Latex

def lprint(*args,**kwargs):
    """Pretty print arguments as LaTeX using IPython display system 
    
    Parameters
    ----------
    args : tuple 
        What to print (in LaTeX math mode)
    kwargs : dict 
        optional keywords to pass to `display` 
    """
    display(Latex('$$'+' '.join(args)+'$$'),**kwargs)
    
def myDiff(formula):
    return sqrt((formula.diff(L) * dL)**2 + (formula.diff(W) * dW)**2)

def myDiffWithCorr(formula, name = "", printNow = False):
    dd = sqrt((formula.diff(L) * dL)**2 + (formula.diff(W) * dW)**2 + 2*(formula.diff(L)*formula.diff(W)*(sigCorr**2)))
    if(printNow):
        lprint(latex(Eq(symbols('sigma_'+name), dd)))
    fd = lambdify((L,dL,W,dW,sigCorr),dd)
    return dd, fd
    
def diff_and_print(formula, name = ""):
    # Calculate uncertainty and print original relation/formula and the uncertainty
    dd = myDiff(formula)
    lprint(latex(Eq(symbols(name),formula)))
    lprint(latex(Eq(symbols('sigma_'+name), dd)))
    
def lambdifyFormula(formula, *args, name = ""):
    # Turn expression into numerical functions 
    f = lambdify((L,W),formula)
    d = myDiff(formula)
    fd = lambdify((L,dL,W,dW),d)
    return f, fd
```


```python
# Import SymPy: 
from sympy import * 
    
# Define variables:
L,W,P,A,D = symbols("L, W, P, A, D")
dL,dW,dP,dA,dD = symbols("sigma_L, sigma_W, sigma_P, sigma_A, sigma_D")

# Define relations:
# Perimeter:
P = 2*L + 2*W
# Area:
A = L*W
# Diagonal
D = sqrt(L**2 + W**2)

# Try writing a simple function to not repeat yourself! (See cell above)
diff_and_print(P,"P")
diff_and_print(A,"A")
diff_and_print(D,"D")

dP = myDiff(P)
dA = myDiff(A)
dD = myDiff(D)

# Turn expressions into numerical functions 
fP, fdP = lambdifyFormula(P,"P")
fA, fdA = lambdifyFormula(A,"A")
fD, fdD = lambdifyFormula(D,"D")

# Define values and their errors
vL, vdL = mu1,sig1
vW, vdW = mu2,sig2

# Numerically evaluate expressions and print 
vP = fP(vL,vW)
vdP = fdP(vL,vdL,vW,vdW)
vA = fA(vL,vW)
vdA = fdA(vL,vdL,vW,vdW)
vD = fD(vL,vW)
vdD = fdD(vL,vdL,vW,vdW)

lprint(fr'P = ({vP:.1f} \pm {vdP:.1f})\,\mathrm{{m}}')
lprint(fr'A = ({vA:.1f} \pm {vdA:.1f})\,\mathrm{{m}}')
lprint(fr'D = ({vD:.1f} \pm {vdD:.1f})\,\mathrm{{m}}')

#Adding correlations (and also derivation, printing and lambdifying)
sigCorr = symbols("sigma_LW")
rho = symbols("rho_LW")

dP, fdP = myDiffWithCorr(P, "P", True)
dA, fdA = myDiffWithCorr(A, "A", True)
dD, fdD = myDiffWithCorr(D, "D", True)

sCorr = sqrt(rho*dL*dW)
fSC = lambdify((rho,dL,dW),sCorr)

vSigmaCorr = fSC(rho12,vdL,vdW)

# Numerically evaluate expressions and print 
vdP = fdP(vL,vdL,vW,vdW,vSigmaCorr)
vdA = fdA(vL,vdL,vW,vdW,vSigmaCorr)
vdD = fdD(vL,vdL,vW,vdW,vSigmaCorr)

lprint(fr'P = ({vP:.1f} \pm {vdP:.1f})\,\mathrm{{m}}')
lprint(fr'A = ({vA:.1f} \pm {vdA:.1f})\,\mathrm{{m}}')
lprint(fr'D = ({vD:.1f} \pm {vdD:.1f})\,\mathrm{{m}}')
```


$$P = 2 L + 2 W$$



$$\sigma_{P} = \sqrt{4 \sigma_{L}^{2} + 4 \sigma_{W}^{2}}$$



$$A = L W$$



$$\sigma_{A} = \sqrt{L^{2} \sigma_{W}^{2} + W^{2} \sigma_{L}^{2}}$$



$$D = \sqrt{L^{2} + W^{2}}$$



$$\sigma_{D} = \sqrt{\frac{L^{2} \sigma_{L}^{2}}{L^{2} + W^{2}} + \frac{W^{2} \sigma_{W}^{2}}{L^{2} + W^{2}}}$$



$$P = (8.6 \pm 0.9)\,\mathrm{m}$$



$$A = (2.8 \pm 0.8)\,\mathrm{m}$$



$$D = (3.6 \pm 0.4)\,\mathrm{m}$$



$$\sigma_{P} = \sqrt{4 \sigma_{L}^{2} + 8 \sigma_{LW}^{2} + 4 \sigma_{W}^{2}}$$



$$\sigma_{A} = \sqrt{L^{2} \sigma_{W}^{2} + 2 L W \sigma_{LW}^{2} + W^{2} \sigma_{L}^{2}}$$



$$\sigma_{D} = \sqrt{\frac{L^{2} \sigma_{L}^{2}}{L^{2} + W^{2}} + \frac{2 L W \sigma_{LW}^{2}}{L^{2} + W^{2}} + \frac{W^{2} \sigma_{W}^{2}}{L^{2} + W^{2}}}$$



$$P = (8.6 \pm 1.1)\,\mathrm{m}$$



$$A = (2.8 \pm 0.9)\,\mathrm{m}$$



$$D = (3.6 \pm 0.4)\,\mathrm{m}$$



```python
import numpy as np                                     # Matlab like syntax for linear algebra and functions
import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab
import seaborn as sns                                  # Make the plots nicer to look at
from iminuit import Minuit                             # The actual fitting tool, better than scipy's
import sys                                             # Modules to see files and folders in directories
```


```python
sys.path.append('../../../External_Functions')
from ExternalFunctions import Chi2Regression
from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure
```

## Error propagation - Simulation

Now we want to try to see, if we can solve the above error propagation problem using simulation. The method is relatively straight forward: You simply take "realistic" values of the input parameters x (here Length (x1) and Width (x2)), calculate the resulting value y (here Perimeter, Area, and Diagonal), and do this many times. The resulting distribution of y should be centered around the value y(x1,x2), and the standard deviation should reflect the uncertainty in y from the uncertainties in the input variables.

This is a much more clumsy way of calculating the uncertainty, but comes with the advantage, that if the resulting uncertainty is not Gaussian, then one can actually see this (i.e. it avoids the assumptions used in the usual error propagation formula).


```python
# First we set the parameters of the program:
N_exp = 10000           # Number of "experiments" (i.e. drawing from random distributions)
save_plots = False
r = np.random
r.seed(42)
```

Define parameters for two random numbers (Gaussianly distributed):


```python
mu1   =  3.5
sig1  =  0.4
mu2   =  0.8
sig2  =  0.2
rho12 =  0.5           # Correlation parameter!
```


```python
if not (-1.0 <= rho12 <= 1.0): 
    raise ValueError(f"Correlation factor not in interval [-1,1], as it is {rho12:6.2f}")
```

Now we calculate numbers that allows the transform from uncorrelated variables `u` and `v` to correlated random numbers `x1` and `x2` below (see Barlow page 42-44 for method):


```python
# Define the parameters needed for the transformation:
theta = 0.5 * np.arctan( 2.0 * rho12 * sig1 * sig2 / ( np.square(sig1) - np.square(sig2) ) )
sigu = np.sqrt( np.abs( (((sig1*np.cos(theta))**2) - (sig2*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )
sigv = np.sqrt( np.abs( (((sig2*np.cos(theta))**2) - (sig1*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )

# Produce random numbers with the (possible) correlation:
u = r.normal(0.0, sigu, N_exp)
v = r.normal(0.0, sigv, N_exp)
x1_all = mu1 + np.cos(theta)*u - np.sin(theta)*v
x2_all = mu2 + np.sin(theta)*u + np.cos(theta)*v
x12_all = np.array([x1_all, x2_all])

y_all = x1_all - 2*x2_all         # Silly formula - you have to put this in yourself!
```

Note that the above is nothing more than a matrix multiplication written out! Also note that the absolute value is taken before the square root to avoid `np.sqrt(x)` with `x<0`.

***
Plot both input distribution and resulting 2D-histogram on screen:


```python
fig, ax = plt.subplots(figsize=(14, 8))
counts, xedges, yedges, im = ax.hist2d(x1_all, x2_all, bins=[120, 80], range=[[0.0, 6.0], [-1.0, 3.0]], cmin=1)
ax.plot([0.0, 6.0], [0.0, 0.0], "--k")   # NOTE: This draws a line from [x1, x2], [y1, y2] with dashed line ("--") and in black ("k")
fig.colorbar(im) # ticks=[-1, 0, 1]

ax.set(title='Histogram of lengths (x) and widths (y)',
       xlabel='x', 
       ylabel='y',
       aspect='equal', # NOTE: This forces the x- and y-axis to have the SAME scale!!!
      )

d = {'Entries': len(x12_all),
     'Mean x': x1_all.mean(),
     'Mean y': x2_all.mean(),
     'Std x': x1_all.std(ddof=1),
     'Std y': x2_all.std(ddof=1),
    }

text = nice_string_output(d, extra_spacing=2, decimals=3)
add_text_to_ax(0.02, 0.97, text, ax, fontsize=15);

fig.tight_layout()
fig

if save_plots :
    fig.savefig("Dist_2Dgauss.pdf")
```

Now we bin `y_all` and fit it with a Gaussian distribution:


```python
def gaussian(x, N, mu, sigma):
    return N * binwidth / (sigma*np.sqrt(2*np.pi)) * np.exp(-0.5* (x-mu)**2/sigma**2)
```


```python
nbins = 100
xmin, xmax = 0.0, 4.0
binwidth = (xmax-xmin)/nbins

fig2, ax2 = plt.subplots(figsize=(16, 6));
counts, bin_edges, _ = ax2.hist(y_all, nbins, range=(xmin, xmax), histtype='step', linewidth=2)
bin_centers = (bin_edges[1:] + bin_edges[:-1])/2
s_counts = np.sqrt(counts)
```

This is the distribution of "whatever you put into it" (initially x1-2*x2), which shows what output you get and what uncertainty to expect (given by the width - think about this!). We can thus get the result by simply recording the mean and width (SD):


```python
mean = y_all.mean()
std = y_all.std(ddof=1)
print(f"  Mean = {mean:5.3f},    Std = {std:5.3f}")
```

      Mean = 1.895,    Std = 0.402


However, we are in principle not even sure, if this distribution is Gaussian, so in order to check this, we draw a Gaussian on top using the above mean and standard deviation values.


```python
xaxis = np.linspace(xmin, xmax, 1000)
yaxis = gaussian(xaxis, N_exp, mean, std)
ax2.plot(xaxis, yaxis, linewidth=2)
fig2.tight_layout()
fig2
```


```python
if save_plots:
    fig2.savefig("Dist_ErrorProp.pdf")
```

---

# Questions:

0. First solve the problem of obtaining the Perimeter, Area & Diagonal with uncertainty ANALYTICALLY.

**Example answer 0**: 

See cell 4. No need to differentiate by hand, SymPy can do it easily for you.

---

1. Now look at the program, and assure yourself that you understand what is going on. Put in the correct expression for y in terms of x1=L and x2=W in order to calculate the perimeter, area, and diagonal length, and run the program. Does the output correspond well with the results you expected from your analytical calculations to begin with?

**Example answer 1**:
Yes, both the mean value and the uncertainty are in a very good agreement to the results from analytical calculations. In the area calculation, there is a slight discrepancy, simply because multiplying two Gaussians does not give (exactly) a Gaussian... only when adding! 
Also, one little thing to note is, that the width can actually come out negatively in rare cases (i.e. beyond 4 sigma low). When producing a simulation of cases, one should of course be aware of such things!

---

2. Imagine that you wanted to know the central value and uncertainty of y1 and y2, given the
   same above PDFs for `x1`=$L$ and `x2`=$W$:
   
     `y1 = log(square(x1*tan(x2))+sqrt((x1-x2)/(cos(x2)+1.0+x1)))`
     
     `y2 = 1.1+sin(20*x1)`

   Get the central value of y, and see if you can quickly differentiate this with
   respect to `x1` and `x2`, and thus predict what uncertainty to expect for y using
   the error propagation formula. It is (for once) OK to give up on the first expression :-)
   Next, try to estimate the central value and uncertainty using random numbers
   like above - do you trust this result more? And are the distributions Gaussian?

**Example answer 2**:
The first expression "y1" looks horrible, and it is! This is an obvious case where simulation is probably a better way forward, and at least a good cross check.
The second expression "y2" is simple, but made to make the error propagation formula break down! The function varies wildly, and as a result, the error propagation formula assumption is broken (badly). The resulting error distribution is also very far from Gaussian, as can be seen below.

For code solution example, see cell 19 below.

### Advanced questions:

3. Try to generate `x1` and `x2` with non-linear correlation, which yields zero linear correlation,
   and see that despite not having any linear correlation, the result on perimeter, area, and diagonal
   length is still affected.



```python
y1,y2,x1,x2 = symbols('y_1,y_2,x_1,x_2')
dy1,dy2,dx1,dx2 = symbols('sigma_y_1,sigma_y_2,sigma_x_1,sigma_x_2')

y1 = log(sqrt(x1+tan(x2)))+sqrt((x1-x2))/(cos(x2)+1+x1)
lprint(latex(Eq(symbols('y_1'),y1)))

y2 = 1.1 + sin(20*x1)
lprint(latex(Eq(symbols('y_2'),y2)))

dy1 = sqrt((y1.diff(x1) * dx1)**2 + (y1.diff(x2) * dx2)**2)
lprint(latex(Eq(symbols('sigma_y_1'),dy1)))

dy2 = sqrt((y2.diff(x1) * dx1)**2 + (y2.diff(x2) * dx2)**2)
lprint(latex(Eq(symbols('sigma_y_2'),dy2)))

fy1 = lambdify((x1,x2),y1)
fy2 = lambdify((x1,x2),y2)

fdy1 = lambdify((x1,dx1,x2,dx2),dy1)
fdy2 = lambdify((x1,dx1,x2,dx2),dy2)

vx1,vdx1 = mu1,sig1
vx2,vdx2 = mu2,sig2

vy1 = fy1(vx1,vx2)
vdy1 = fdy1(vx1,vdx1,vx2,vdx2)
vy2 = fy2(vx1,vx2)
vdy2 = fdy2(vx1,vdx1,vx2,vdx2)

lprint(fr'y1 = ({vy1:.1f} \pm {vdy1:.1f})\,\mathrm{{m}}')
lprint(fr'y2 = ({vy2:.1f} \pm {vdy2:.1f})\,\mathrm{{m}}') 

u = r.normal(0.0, sigu, N_exp)
v = r.normal(0.0, sigv, N_exp)
x1_all = mu1 + np.cos(theta)*u - np.sin(theta)*v
x2_all = mu2 + np.sin(theta)*u + np.cos(theta)*v
x12_all = np.array([x1_all, x2_all])

y_all = 1.1 + np.sin(20*x1_all)

xmin, xmax = 0.0, 3.0

fig3, ax3 = plt.subplots(figsize=(10, 6));
counts, bin_edges, _ = ax3.hist(y_all, 100, range=(xmin, xmax), histtype='step');
bin_centers = (bin_edges[1:] + bin_edges[:-1])/2
s_counts = np.sqrt(counts)

print(f"  Mean = {y_all.mean():5.3f},    RMS = {y_all.std(ddof=1):5.3f}")
```

830590f1605bf2ec965aeaa4350985f7b11ae06b
ipynb
Jupyter Notebook
AppStat2022/Week1/ExampleSolutions/ErrorPropagation/ErrorPropagation_ExampleSolution.ipynb
EiriniNBI/eirinipy
84fcf539c2a35399858e407db89ba6a1f3c73f32
MIT
AppStat2022/Week1/ExampleSolutions/ErrorPropagation/ErrorPropagation_ExampleSolution.ipynb
EiriniNBI/eirinipy
84fcf539c2a35399858e407db89ba6a1f3c73f32
MIT
AppStat2022/Week1/ExampleSolutions/ErrorPropagation/ErrorPropagation_ExampleSolution.ipynb
EiriniNBI/eirinipy
84fcf539c2a35399858e407db89ba6a1f3c73f32
MIT
Qwen/Qwen-72B
1. YES
2. YES