\documentclass[12pt]{article}
\usepackage{amsmath,amsfonts,amsthm,html,fullpage}

\newcommand{\R}{\mathbb{R}}
\def\e{{\mathrm e}}

\begin{document}
\title{Propagation of Error or Uncertainty}
\author{Marcel Oliver}
\date{December 4, 2015}
\maketitle

\theoremstyle{definition}
\newtheorem{example}{Example}
\newtheorem*{warning}{Warning}
\newtheorem*{remark}{Remark}

\section{Introduction}

\subsection{Motivation}

All measurements are subject to error or uncertainty.  Common causes
are noise or external disturbances, imperfections in the experimental
setup and the measuring devices, coarseness or discreteness of
instrument scales, unknown parameters, and model errors due to
simplifying assumptions in the mathematical description of an
experiment.  An essential aspect of scientific work, therefore, is
quantifying and tracking uncertainties from setup, measurements, all
the way to derived quantities and resulting conclusions.

In the following, we can only address some of the most common and
simple methods of error analysis.  We shall do this mainly from a
calculus perspective, with some comments on statistical aspects later
on.


\subsection{Absolute and relative error}

When measuring a quantity with true value $x_{\operatorname{true}}$,
the measured value $x$ may differ by a small amount $\Delta x$.  We
speak of $\Delta x$ as the \emph{absolute error} or \emph{absolute
uncertainty} of $x$.  Often, the magnitude of error or uncertainty is
most naturally expressed as a fraction of the true or the measured
value $x$.  Since the true value $x_{\operatorname{true}}$ is
typically not known, we shall define the \emph{relative error} or
\emph{relative uncertainty} as $\Delta x/x$.

In scientific writing, you will frequently encounter measurements
reported in the form $x=3.3 \pm 0.05$.  We read this as $x=3.3$ and
$\Delta x=0.05$.

\subsection{Interpretation}

Depending on context, $\Delta x$ can have any of the following three
interpretations which are different, but in simple settings lead to
the same conclusions.

\begin{enumerate}
\item An exact specification of a deviation in one particular trial
(instance of performing an experiment), so that
\[
  x_{\operatorname{true}} = x + \Delta x \,.
\]
Of course, typically there is no way to know what this true value
actually is, but for the analysis of error propagation, this
interpretation is useful.

\item An error bound (which may or may not be strict).  In this
interpretation, we understand that $x_{\operatorname{true}}$ may lie
anywhere in the interval $[x-\Delta x, x+\Delta x]$. 

\item Uncertainty in the statistical sense.  If repeated measurements
of $x$ are expected to follow a probability distribution, then
$\Delta x$ is usually taken as the \emph{sample standard deviation} of
a set of measurements.

\end{enumerate}
Our focus is on analyzing how the error or uncertainty in the input of
a calculation propagates through the calculation, leading to error or
uncertainty in the output.

\section{Error propagation in one variable}
\label{s.1d}

Suppose that $x$ is the result of a measurement and we are calculating
a dependent quantity
\begin{equation}
  y = f(x) \,.
\end{equation}
Knowing $\Delta x$, we must derive $\Delta y$, the associated error or
uncertainty of $y$.

Let us recall the equation for the tangent line to $f$ at point $x$,
\begin{equation}
  \ell(x+\Delta x) = f(x) + f'(x) \, \Delta x \,.
\end{equation}
We use the tangent line equation as a (linear) approximation to $f$;
when $\Delta x$ is not too large, we expect this approximation to be
good.  Thus,
\begin{equation}
  \Delta y 
  = f(x+\Delta x) - f(x)
  \approx \ell(x+\Delta x) - f(x)
  = f'(x) \, \Delta x \,. 
\end{equation}
If we think of $\Delta x$ as a positive number which specifies the
(expected) magnitude of error, we shall write the error propagation
formula in the form
\begin{equation}
  \Delta y \approx \lvert f'(x) \rvert \, \Delta x \,.
  \label{e.1d}
\end{equation}

\begin{example}
If you measure $x=49 \pm 4$, what should you report for $y=\sqrt x$
together with its uncertainty?  To answer this question, we use
\eqref{e.1d} to compute
\begin{equation}
  \Delta y \approx \frac1{2 \sqrt x} \, \Delta x
           = \frac1{2 \sqrt 49} \cdot 4 
           = \frac27 
           \approx 0.3 \,.
\end{equation}
Thus, $y = 7 \pm 0.3$.  (Note that it does not make sense to compute
an estimate for $\Delta y$ to a large number of significant digits as
there is no basis for reporting the resulting error with great
accuracy.)
\end{example}

\begin{example}[Area of circle]
Suppose you can measure the radius of a circle with relative
uncertainty $\Delta r/r = 5\%$.  We ask for the associated relative
uncertainty of the computed area
\begin{equation}
  A = \pi \, r^2 \,.
\end{equation}
By \eqref{e.1d}, $\Delta A \approx \pi \, 2r \, \Delta r$.  Thus, the
relative uncertainty of the area is given by
\begin{equation}
  \frac{\Delta A}{A}
  \approx \frac{\pi \, 2r \, \Delta r}{\pi \, r^2}
  = 2 \, \frac{\Delta r}r = 10\% \,.
\end{equation}
\end{example}

\begin{example}[Allometric relation] \label{ex.allo}
More generally, when $y$ is allometrically related to $x$, i.e., 
\begin{equation}
  y = k \, x^\alpha \,,
\end{equation}
the same type of computation (fill in the details!)\ gives
\begin{equation}
  \frac{\Delta y}{y} \approx \alpha \, \frac{\Delta x}x \,.
\end{equation}
\end{example}

\begin{example}[Exponential relation]
When $y$ is exponentially related to $x$, i.e.,
\begin{equation}
  y = c \, \e^{\lambda x} \,,
\end{equation}
we find
$\Delta y \approx \lvert c \lambda \rvert \, \e^{\lambda x} \, \Delta
x$, so that
\begin{equation}
  \frac{\Delta y}{y}
  \approx \frac{\lvert c \lambda \rvert \, \e^{\lambda x} \, \Delta x}%
               {c \, \e^{\lambda x}}
  = \lvert \lambda \rvert \, \Delta x \,.
\end{equation}
(We take $c>0$ for simplicity.)  Thus, in an exponential relation it
is most natural to express the relative error of the output in terms
of the absolute error of the input data!
\end{example}


\section{Error propagation in several variables}

\subsection{Linear approximation and simple error propagation}

Suppose, for simplicity, that $z=f(x,y)$ is a function of two
variables.  Once again, we may use linear approximation, where
\begin{equation}
  \ell(x+\Delta x, y+\Delta y)
  = f(x,y) + \frac{\partial f}{\partial x}(x,y) \, \Delta x
           + \frac{\partial f}{\partial y}(x,y) \, \Delta y
\end{equation}
so that
\begin{equation}
  \Delta z 
  \approx \ell(x+\Delta x, y+\Delta y) - \ell(x,y)
  = \frac{\partial f}{\partial x}(x,y) \, \Delta x
    + \frac{\partial f}{\partial y}(x,y) \, \Delta y \,.
  \label{e.linearapprox}
\end{equation}
As in Section~\ref{s.1d}, this leads to the error propagation formula
\begin{equation}
  \Delta z 
  \approx \biggl\vert \frac{\partial f}{\partial x}(x,y) \biggr\vert \, 
          \Delta x
        + \biggl\vert \frac{\partial f}{\partial y}(x,y) \biggr\vert \, 
          \Delta y \,.   
  \label{e.2d1}
\end{equation}
This formula is in fact the best we can do when $\Delta x$ and
$\Delta y$ represent true interval constraints on the error.  It is,
however, unnecessarily pessimistic in the case when the true error can
be though of as lying in a ``ball of uncertainty'' with a known radius
in the $x$-$y$ plane.  This is the case, for example, when measuring
distances in the plane where the estimated uncertainty should not
depend on the particular choice of coordinates.  Formula
\eqref{e.2d1}, however, is not invariant under rotation of the
coordinate axes.

\subsection[Not used]{A geometric view on error
propagation\footnote{Supplementary material, skip on
first reading.}}

Suppose you know that the true values $x_{\operatorname{true}}$ and
$y_{\operatorname{true}}$ lie in an ellipse with center $(x,y)$ and
semi-axes $\Delta x$ and $\Delta y$.\footnote{The analogous geometric
picture behind \eqref{e.2d1} is that $x_{\operatorname{true}}$ and
$y_{\operatorname{true}}$ lie in a box with center $(x,y)$ and
semi-widths $\Delta x$ and $\Delta y$.}  Then the maximum error in
\eqref{e.linearapprox} is the maximum of
\begin{equation}
  F(\xi,\eta) = a \, \xi + b \, \eta
  \label{e.F}
\end{equation}
with
\begin{equation}
  a = \frac{\partial f}{\partial x}(x,y) \, \Delta x
  \qquad \text{and} \qquad
  b = \frac{\partial f}{\partial y}(x,y) \, \Delta y \,.
\end{equation}
under the constraint that $(\xi,\eta)$ is a point on the unit circle,
i.e., 
\begin{equation}
  g(\xi,\eta) \equiv \xi^2 + \eta^2 = 1 \,.
  \label{e.constraint}
\end{equation}
This problem is easily solved using Lagrange multipliers.  Setting up
the Lagrange equations,
\begin{equation}
  \frac{\partial F}{\partial \xi} 
  = \lambda \, \frac{\partial g}{\partial \xi}
  \qquad \text{and} \qquad 
  \frac{\partial F}{\partial \eta} 
  = \lambda \, \frac{\partial g}{\partial \eta} \,,
\end{equation}
we find $a = \lambda \, 2 \xi$ and $b = \lambda \, 2 \eta$, so that
\begin{equation}
  \frac{a}b = \frac\xi\eta \,.
\end{equation}
Using this relation to eliminate $\eta$ from the constraint
$\xi^2+\eta^2=1$, we find
\begin{equation}
  \xi^2 + \xi^2 \, \frac{b^2}{a^2} = 1
  \qquad \text{so that} \qquad
  \xi^2 = \frac{a^2}{a^2+b^2} \,.
\end{equation}
A similar expression is obtained by eliminating $\xi$ and solving for
$\eta$ so that, taking square roots, we have
\begin{equation}
  \xi = \frac{\pm a}{\sqrt{a^2+b^2}}
  \qquad \text{and} \qquad
  \eta = \frac{\pm b}{\sqrt{a^2+b^2}} \,.
\end{equation}
Inserting this expression into \eqref{e.F} and maximizing the value of
$F$ by choosing the positive sign on each root, we obtain
\begin{equation}
  F_{\max}
  = a \, \frac{a}{\sqrt{a^2+b^2}}
  + b \, \frac{b}{\sqrt{a^2+b^2}}
  = \frac{a^2 + b^2}{\sqrt{a^2+b^2}}
  = {\sqrt{a^2+b^2}} \,.
\end{equation}
Substituting the definitions of $a$ and $b$ back into this expression
yields a formula for $\Delta z$ which we write out in the next
section.\footnote{Using vector calculus, there is an alternative, much
shorter argument.  The estimate for $\Delta z$ is the maximal rate of
ascend of $\phi(\xi,\eta) = f(x + \xi\Delta x, y + \eta \Delta y)$, so
it is given by the length of the gradient of $\phi$.  This directly
yields \eqref{e.2d2}.  From this point of view, the argument in this
section can be seen as a proof that the gradient is the vector
pointing in the direction of steepest ascend, its magnitude giving the
rate of change in this direction.  This goes beyond what is covered in
this class, but can be found in any standard text on vector calculus,
e.g.\ \cite[Chapter~13]{K}.}

\subsection{Propagation of independent uncertainties}

The geometrically motivated computation in the previous section leads
to the following error propagation formula for $z=f(x,y)$:
\begin{equation}
  \Delta z 
  \approx \sqrt{\biggl(\frac{\partial f}{\partial x} (x,y) \biggr)^2 \,
                \Delta x^2 +
                \biggl(\frac{\partial f}{\partial y} (x,y) \biggr)^2 \,
                \Delta y^2} \,.
  \label{e.2d2}
\end{equation}
The formula extends to more than two variables in the obvious way.
Moreover, probability theory shows that this formula is also correct
in the case where $x$ and $y$ are \emph{independent random variables}
with standard deviations $\Delta x$ and $\Delta y$.  A
first-principles discussion is beyond the scope of this class, but the
concept is very frequently used in practice.  We will discuss some
examples and implications.

\begin{example} \label{ex.5} If you measure $x=6.0\pm 0.5$ and
$y=2.0\pm 0.1$ independently, what should you report for $z=x/y$
together with its uncertainty?  We compute
\begin{equation}
  \frac{\partial z}{\partial x} = \frac1y
  \qquad \text{and} \qquad
  \frac{\partial z}{\partial y} = -\frac{x}{y^2}
\end{equation}
and substitute these into \eqref{e.2d2} to find
\begin{equation}
  \Delta z 
  \approx \sqrt{\biggl( \frac1{2.0} \, 0.5 \biggr)^2
              + \biggl( \frac{6.0}{2.0^2} \, 0.1 \biggr)^2}
  \approx 0.3 \,.
\end{equation}
Thus, $z=3.0\pm 0.3$.
\end{example}

\begin{example}[Sum of independently measured quantities]
\label{ex.sum}
Suppose you measure $x$ and $y$ independently with uncertainties
$\Delta x$ and $\Delta y$.  What is the uncertainty of their sum
$z=x+y$?  Here, the partial derivatives $\partial z/\partial x
= \partial z/\partial y = 1$, so that
\begin{equation}
  \Delta z \approx \sqrt{\Delta x^2 + \Delta y^2} \,.
\end{equation}
This formula suggests that, geometrically, we may think of the
absolute uncertainty $\Delta z$ as the Euclidean length of the
absolute uncertainty vector $(\Delta x, \Delta y)$.
\end{example}

\begin{example}[Product of independently measured quantities]
\label{ex.prod}
In the setting of Example~\ref{ex.sum}, let us look at the product
$z=xy$.  Here, $\partial z/\partial x=y$ and $\partial z/\partial
y=x$, so that
\begin{equation}
  \Delta z \approx \sqrt{y^2 \, \Delta x^2 + x^2 \, \Delta y^2} \,.
\end{equation}
Dividing through by $z$ on both sides, we find that 
\begin{equation}
  \frac{\Delta z}z 
  \approx \sqrt{\frac{\Delta x^2}{x^2} + \frac{\Delta y^2}{y^2}} \,.
\end{equation}
In other words, when taking products, it is the \emph{relative}
uncertainty $\Delta z/z$ which is the Euclidean length of the
\emph{relative} uncertainty vector $(\Delta x/x,\Delta y/y)$.
\end{example}

\begin{example}[Atwood machine]
The acceleration of two masses $m_1$ and $m_2$ in an Atwood
machine\cite{a} is given by the formula
\begin{equation}
  a = g \, \frac{m_1-m_2}{m_1+m_2} \,,
\end{equation}
where $g = 9.81 \, \text{m}/\text{s}^2$ is the constant of gravity.
When $m_1=100 \pm 1 \, \text{g}$ and $m_2=50 \pm 1 \, \text{g}$,
estimate the uncertainty $\Delta a$.  After a short computation (fill
in the details!), we find
\begin{equation}
  \frac{\partial a}{\partial m_1} 
  = g \, \frac{2\, m_2}{(m_1+m_2)^2}
  \qquad \text{and} \qquad
  \frac{\partial a}{\partial m_2} 
  = g \, \frac{-2\, m_1}{(m_1+m_2)^2}
\end{equation}
so that
\begin{equation}
  \Delta z 
  \approx \frac{2 \, g}{(m_1+m_2)^2}
          \sqrt{(m_2 \, \Delta m_1)^2 + (m_1 \, \Delta m_2)^2} \,.
\end{equation}
Inserting the numbers gives $a = 3.27\, \text{m}/\text{s}^2$ and
$\Delta a = 0.097\, \text{m}/\text{s}^2$.  We might want to report
this result as
$a = 3.27 \pm 0.1 \text{m}/\text{s}^2$.\footnote{Example adapted from
\url{https://courses.washington.edu/phys431/propagation_errors_UCh.pdf}.
There, it is suggested to report the final result as
$a = 3.3 \pm 0.1 \text{m}/\text{s}^2$, but then the reported
confidence interval is clearly outside of the computed confidence
interval, so I do not recommend rounding too aggressively.}
\end{example}

\begin{example}[Volume of cylinder with independent errors]
Suppose that you measure the radius $r$ and the height $h$ of a
cylinder each with an independent relative error of $1\%$.  What is
the relative error of the computed volume $V=\pi h r^2$?  You could
apply \eqref{e.2d2} directly and go from there (try it!), but at this
point it is quicker to refer to Example~\ref{ex.allo} which shows that
$r^2$ has a relative error of $2\%$.  Further, the relative error of
$h$ and the relative error of $r^2$ add according to
Example~\ref{ex.prod}, so the relative error of $V$ is $\sqrt{5}\%$.
Thus, the relative error increases by a factor $\sqrt 5 \approx 2.24$.
\end{example}

\begin{example}[Volume of cylinder with dependent errors]
Now suppose that the errors in the measurements are actually
dependent.  Such a situation might arise, for example, due to
temperature sensitivity of the equipment which affects each
measurement in exactly the same way.  In this case, we need to refer
to \eqref{e.2d1}.  Leaving out the details of the computation (try
it!), we arrive at an increase of the relative error by a factor of
$3$.
\end{example}

\begin{warning}
The last two examples show that assuming independence of uncertainty
when this is not actually guaranteed leads to an underestimate of the
uncertainty of the result.  This is dangerous as it may lead to
conclusions which are not supported by the measured data.
\end{warning}

\begin{remark}
When the quantities involved can be described in statistical terms,
then \eqref{e.2d2} may be corrected by covariance terms which
represent correlations between the different errors, see \cite{c1}.
This, however, is beyond the scope of this class and also rarely used
in practice as the required covariances are typically difficult to
estimate.
\end{remark}

\subsection{Repeated experiments}

Sometimes, an experiment is repeated many times and the values
recorded are averaged over all trials.  Suppose we have collected $N$
readings $x_1, \dots, x_N$ and we compute their mean
\begin{equation}
  \overline x = \frac{x_1 + \dots + x_N}N \,.
\end{equation}
When an individual reading has uncertainty
$\Delta x_1= \dots =\Delta x_N \equiv \Delta x$, what is the
uncertainty of $\overline x$?  Extending \eqref{e.2d1} to $N$
variables, we find
\begin{equation}
  \Delta \overline x
  \approx \frac{\Delta x_1}N + \dots + \frac{\Delta x_N}N 
  = \Delta x \,,
  \label{e.dlx1}
\end{equation}
so there is no advantage in repeating the experiment.  However, when
the uncertainties are known to be independent, we may use
\eqref{e.2d2}, extended to $N$ variables, instead, so that
\begin{equation}
  \Delta \overline x
  \approx \sqrt{\biggl(\frac{\Delta x_1}N\biggr)^2 + \dots + 
                \biggl(\frac{\Delta x_N}N\biggr)^2}
  = \frac{\Delta x}{\sqrt N} \,.
  \label{e.dlx2}
\end{equation}
Here, by taking the average over many measurements, errors partially
cancel with high probability and the uncertainty of $\overline x$
converges to zero at a rate of $1/\sqrt N$ as the number of
measurements goes to infinity.  

Note that the benefit of repeated trials may be destroyed by unnoticed
dependencies between successive measurements---it is not helpful to
make the same error twice!



\section{Monte-Carlo testing}

Monte-Carlo testing is a computational technique for estimating the
propagation of uncertainties.  The idea is to perform the computation
many times with randomly perturbed input and use the statistics of the
output to determine its uncertainty.  With Scientific Python/Pylab or
similar mathematical software, Monte-Carlo testing is extremely easy
and quick so that it should be part of the tool set of every working
scientist.  We describe the procedure in the setting of
Example~\ref{ex.5}.

Rather than assigning a single value for each input variable, we
create a vector of values that follow a probability distribution whose
mean equals the measured value and whose standard deviation equals the
corresponding uncertainty.  Typically, we choose a normal distribution
via the \verb+normal+ function; its first argument is the mean, the
second argument is the standard deviation, and the third argument the
number of samples which we set to a moderately large integer value:
\begin{verbatim}
  In [1]: x = normal(6,0.5,10000)
  In [2]: y = normal(2,0.1,10000)
\end{verbatim}
Now we can perform the computation \emph{with all samples at once}:
\begin{verbatim}
  In [3]: z = x/y
\end{verbatim}
The uncertainty of the output is then given by the standard deviation
of $z$,
\begin{verbatim}
  In [4]: std(z)
  Out[4]: 0.29531173914816444
\end{verbatim}
This corresponds well with the result suggested by our error
propagation formula.

Monte-Carlo testing has a number of advantages: it can be done with
ease in a routine manner, it does not use linearization, thus will
also work as expected in situations when linear approximation is not
appropriate, and it is possible to include dependencies between the
different errors in the input statistics if necessary.

The drawbacks are that it requires many repeated computations.  This
is not an issue when testing simple function evaluations on
contemporary computers as in the example above.  However, when dealing
with much more complex models, computational cost may become an issue.
Another drawback is that the procedure is purely numerical.  Thus, it
is not possible to see how the output uncertainty depends on
parameters of the setup, which is sometimes of interest.  For this,
the calculus-based error propagation formulas are better suited.





\begin{thebibliography}{9}

\bibitem{K}
  H.J. Keisler, ``Elementary Calculus,'' 2nd edition, Prindle, Weber
  \& Schmidt, 1986. \\
  {\tiny\url{http://www.math.wisc.edu/~keisler/calc.html}}

\bibitem{T}
  J.R. Taylor, "An Introduction to Error Analysis,'' 2nd edition, University
  Science Books, 1997. 
  {\tiny\url{http://people.westminstercollege.edu/faculty/kphilippi/An_Introduction_To_Error_Analysis_The_Study_Of_Uncertainties_In_Physical_Measurements_Taylor_John.pdf}}

\bibitem{a}
  "Atwood machine." Wikipedia, The Free Encyclopedia.  December 6,
  2015, 16:52 UTC. \\
  {\tiny\url{http://en.wikipedia.org/w/index.php?title=Atwood_machine&oldid=670698954}}

\bibitem{c1}
  ``Propagation of uncertainty.'' Wikipedia, The Free Encyclopedia.
  December 5, 2015, 14:35 UTC. \\
  {\tiny\url{http://en.wikipedia.org/w/index.php?title=Propagation_of_uncertainty&oldid=693873085}}




\end{thebibliography}

\end{document}
