9
$\begingroup$

enter image description hereI 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?

$\endgroup$
  • 1
    $\begingroup$ One very enlightening and clear way to approach this is to model the measured value as the real value plus a noise which is drawn from a statistical distribution. $\endgroup$ – DanielSank Aug 12 '15 at 20:24
  • $\begingroup$ I see that you added a picture, but this doesn't help understand the spline issue. I think you're asking how to compute the uncertainty in a split fit to some data. I would have to know the formula for how the spline is being fit. $\endgroup$ – DanielSank Aug 15 '15 at 18:12
  • $\begingroup$ I am using Matlab to produce a spline, the smoothing parameter is 0.9. That is all I know about the spline... $\endgroup$ – cryonole Aug 15 '15 at 18:14
  • $\begingroup$ I'm not going to look up Matlab's spline fitting algorithm for you; you can do some leg work here :-) $\endgroup$ – DanielSank Aug 15 '15 at 18:16
  • $\begingroup$ What's the status of this question? If you want to know specifically how to find the uncertainty in the slope of the fitted spline I would move this to the math site. $\endgroup$ – DanielSank Aug 29 '15 at 15:07
8
$\begingroup$

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$.

$\endgroup$
  • 2
    $\begingroup$ Also, while this is good (+1), I think the OP's problem is that $\Delta t$ is small so the uncertainty in the derivative is large (as it should be, since noise will add more to the slope). But this is all assuming just 2 points. Since the OP must have more points to be fitting a spline anyway, perhaps they should be fitting a (differentiable) model to the full dataset rather than to every adjacent pair of points. $\endgroup$ – user10851 Aug 12 '15 at 22:04
  • $\begingroup$ @ChrisWhite I agree with your comment 100%. The goal of my answer was to answer OP's specific question which was how to compute the uncertainty in the derivative given the uncertainty in the measurements. I try to answer the question asked as much as possible, rather than venture into suggesting how OP can/should address the problem which inspired the Physics.SE post :-) $\endgroup$ – DanielSank Aug 12 '15 at 22:10
  • $\begingroup$ @DanielSank, in your post you show how the uncertainty in the derivative relates so the uncertainty in the measured signal. Uncertainty in $dT_{spline}/dt$ would then be given by $\sqrt(2){\sigma}_{spline}/\Delta t$. Am I thinking it right? $\endgroup$ – cryonole Aug 14 '15 at 1:07
  • $\begingroup$ @cryonole The really correct answer is that you have to take the formula used to compute the spline and transform the distributions in the way illustrated in my answer. However, people usually don't do this because it's laborious. Instead, they assume that the distributions are Gaussian all the way through and just use simple scaling arguments to compute the transformation of the standard deviation (you have have leaned of this as "error propagation" in school). I'm not sure how to comment on your proposed formula until I know what you mean by $\sigma_\text{spline}$. $\endgroup$ – DanielSank Aug 14 '15 at 4:42
  • $\begingroup$ @DanielSank You fit a spline and then evaluate the 'smoothed' signal using the spline coefficients. This signal still has noise, although lesser than the actual signal has. By $\sigma _{spline}$, I meant the standard deviation of this spline approximated/evaluated signal. $\endgroup$ – cryonole Aug 14 '15 at 13:27
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Is to 'bin' to downsample the data? $\endgroup$ – cryonole Aug 13 '15 at 19:47
  • 1
    $\begingroup$ I'm not using "bin" as a technical term. I'm sure there is a technical term for it. It is not downsampling. To downsample by a factor n you just throw away all of the data except for every $n^\textrm{th}$ point. What I am talking about is replacing each $n$ data points with the average of that group of $n$. This takes you from $N$ points to $N/n$ smoothed points. You could also do a "running average" where every data point is replaced by the average of the $n$ points centered on that point so you go from $N$ points to $N$ smoothed points. $\endgroup$ – gleedadswell Aug 13 '15 at 23:04
2
$\begingroup$

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:enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.