Contents

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)\)

(20)\[\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

#%pylab inline --no-import-all

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

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

Setup repeated observations of two variables x, y

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

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)')

Text(0, 0.5, 'p(x)')
../_images/investigating-propagation-of-errors_7_1.png

same thing with an interactive widget

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))

<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\)

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)')

Text(0, 0.5, 'q(x)')
../_images/investigating-propagation-of-errors_11_1.png 
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)')

Text(0, 0.5, 'p(x)')
../_images/investigating-propagation-of-errors_12_1.png

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

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)')

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

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

Check propagation of errors under division

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))

<function __main__.plot_x_over_y(mean_x, std_x, mean_y, std_y, N)>
Transformation properties of the likelihood and posterior Revisiting error propagation with automatic differentiation