How does uncertainty/error propagate with differentiation?

Question

I have a noisy temperature (T) vs. time (t) measurement and I want to calculate dT/dt. If I approximate $dT/dt = \Delta T/\Delta t$ then the noise in the derivative gets too high and the derivative becomes useless. So I fit a smoothing spline (smoothing parameter say 'p') to the measured data and get $dT/dt$ by piecewise differentiation of the spline. Is there a way to obtain uncertainty in this $dT/dt$ based on uncertainty in T?

Answer 1

Model two consecutive measurements as the real values plus some noise. Call the first measured temperature $T_1$ and the second $T_2$. Call the measured noises $\gamma_1$ and $\gamma_2$, and suppose that they are drawn from a distribution $\Gamma(\gamma)$ and are uncorrelated.

The (approximation to the) derivative is

$$\text{Derivative} \approx \frac{(T_2 + \gamma_2) - (T_1 + \gamma_1)}{\Delta t} \, .$$

Note that the derivative is itself a random variable because the $\gamma$'s are random variables. What is the probability distribution of this new random variable?

Focus first on the numerator. Here we have a deterministic part $T_2 - T_1$ and a stochastic part $\gamma_2 - \gamma_1$. The trick thing you may not know is how to figure out the probability distribution of the sum or difference of two random variables; in fact the answer is not at all trivial. Given two random variables $x$ and $y$ with distributions $X(x)$ and $Y(y)$, the random variable $z$ defined by $z = x + y$ has distribution

$$ Z(z) = (X \otimes Y)(z) \equiv \int_{-\infty}^\infty X(w) Y(z - w) \, dw \, .$$

This integral is called a convolution. Anyway, the point is that the probability distribution $P_{\gamma_2 - \gamma_1}$ of $\gamma_2 - \gamma_1$ is the convolution of the distributions of $\gamma_2$ and $-\gamma_1$, which is

$$P_{\gamma_2 - \gamma_1}(\gamma) = \int_{-\infty}^\infty \underbrace{\Gamma(-\gamma')}_{\text{from }-\gamma_1} \underbrace{\Gamma(\gamma - \gamma')}_{\text{from }\gamma_2} \, d \gamma' \, .$$

As an example, suppose the noise is Gaussian distributed with standard deviation $\sigma$,

$$\Gamma(\gamma) = \frac{1}{\sqrt{2 \pi} \sigma} \exp\left(\frac{-\gamma^2}{2 \sigma^2} \right) \, .$$

In this case we can do the integral, and the result is

$$P_{\gamma_2 - \gamma_1}(\gamma) = \frac{1}{\sqrt{2\pi} (\sqrt{2} \sigma)} \exp \left( \frac{-\gamma^2}{2 (\sqrt{2}\sigma)^2}\right) \, ,$$

which is just a Gaussian with standard deviation $\sqrt{2} \sigma$. Now remember we also divide by $\Delta t$, and doing this too modifies the distribution. The result is that the probability distribution is still a Gaussian where the standard deviation turns out to be $\sqrt{2}\sigma / \Delta t$.

So that's your answer: the error in the derivative is completely described by a Gaussian probability distribution with standard deviation $\sqrt{2} \sigma / \Delta t$.

Answer 2

I would avoid using a spline in data analysis in general. The spline draws smooth curves through arbitrary sets of points but destroys a lot of the information in the points and adds extraneous, meaningless "information". Use it for making things pretty, not for analysis, unless you really think your data should follow a power series (which would be rather unusual data).

I like @DanielSank's answer a lot (and I voted it up) as it leads to good ways to characterize the noise in your data. Here is a "rough and ready" answer that you might find easier to use with typical data. It is quite crude and will only be valid if the temperature varies slowly compared to your sampling frequency.

If your temperature is slowly varying compared to your sampling frequency, $1/\Delta t$, then instead of taking a spline you could just "bin" measurements so that every group of $n$ successive measurements is averaged into a single temperature measurement. So you are replacing your $N$ temperature measurements with $N/n$ temperature measurements. You can now estimate the uncertainty of each of these measurements as simply the standard deviation of the mean (SDOM, a.k.a. standard error) of the measurements in that bin. Now you will have much less noisy data and you can do the numerical derivative as you normally would.

Danger Will Robinson

1. Do this only if you are sure that your temperature is slowly varying. Even then it is very crude.

2. If you, later, need to take a Fourier transform because you want to know the spectrum of temperature variation then this binning will have "chopped off" the high frequency part of the spectrum.

Answer 3

Answer to the OP can be obtained easily using the Savitzky-Golay (SG) smoothing-differentiation filter. Suppose we have noisy $n$-point data such as the temperature ($T$) vs. time ($t$) as in the OP. As per the OP we want to smooth the data, find the time rate of change, and the uncertainty that the SG procedure would introduce in the time rate of change. The SG method works like this:

The Savitzky-Golay procedure fits a polynomial of degree $m$ to the data contained in a window of size $2w+1 << n$. Of interest in operating over a particular window is to determine the smooth temperature, the derivative, and the derivative uncertainty ONLY at the center of this window . The window is then moved across the data so that the quantities of interest may be calculated at every data point. To illustrate how this method works I've used $m$=2. Consider the procedure operating on a window centered at point at $j$. In this case the fitting polynomial is: $$T_j=b_{0,j}+b_{1,j}\bar t+b_{2,j}{\bar t}^2$$ where $\bar t=(t-t_j)/\Delta t$ , i.e., the abscissa of a data-point in the window relative to the abscissa of the window center-point , normalized by the spacing between adjacent data points . In simple terms, $\bar t$ represents the abscissa-distance of a data-point at $t$ in a window from the center of that window at $t_j$ in units of the spacing $\Delta t$. With this transformation $b_0$, $b_1/\Delta t$, and $\delta (b_1/\Delta t)$ respectively are the smooth temperature, the temperature derivative with time, and the uncertainty in derivative of temperature at $t_j$ (window center). The fit parameters for the window $j$ are given by:

$$b_{0,j}=\sum_{i=-w}^{w}C_{0i,j}T_{i,j}$$ $$b_{1,j}=\sum_{i=-w}^{w}C_{1i,j}T_{i,j}$$ $$b_{2,j}=\sum_{i=-w}^{w}C_{2i,j}T_{i,j}$$

In these expressions the coefficients $C_{0i,j}$, $C_{1i,j}$, and $C_{2i,j}$ are obtained by the procedure devised by Savitzky Golay. The fit uncertainty for window $j$ is then:

$$\delta_j=\sqrt{\frac{1}{2w-1}\sum_{i=-w}^{w}(T_{i,j}-T_j(t_i))^2}$$ which when propagated through the expression of $b_{1,j}$ yields $\delta {b_{1,j}}$, which the uncertainty in the time derivative of temperature. Note this uncertainty is purely comes from the smoothing procedure so that by choosing certain values of $m$ and $w$ the uncertainty can be made as small as one requires. The result of this procedure with $w$=10 is shown below: