# Standard deviation of several measurements with uncertainties

I have two 2 hours of GPS data with a sampling rate of 1 Hz (7200 measurements). The data are given in the form $(X, X_\sigma, Y, Y_\sigma, Z, Z_\sigma)$, where $N_\sigma$ is the measurement uncertainty.

When I take the mean of all measurements (e.g. the average Z value of those two hours), what is its standard deviation? I can of course calculate the standard deviation from the Z values, but then I neglect the fact that there are known measurement uncertainties...

Edit: The data is all from the same station, and all coordinates are remeasured every second. Due to satellite constellations etc., every measurement has a different uncertainity. The purpose of my analysis is to find the displacement due to the an external event, (i.e. an earthquake). I would like to take the mean for 7200 measurements (2h) before the earthquake and another mean for 2h after the earthquake, and then calculate the resulting difference (in height for example). In order to specifiy the standard deviation of this difference, I need to know the standard deviation of the two means.

• Good question. Even more importantly, the data will be strongly positively correlated over time: that will have a more profound effect on the answer than the variation in the measurement uncertainties.
– whuber
Apr 17, 2014 at 15:21
• Picking up on whuber's comment and Deathkill14's answer, you have not given us enough information to answer properly. It's important to know how the errors in measuring $X,Y,Z$ "work." For example, if the error measuring $X$ was positive at 3 seconds, it it more/less likely to be positive at 4 seconds---i.e. is there serial correlation? Second, if the error in $X$ was positive at 3 seconds, it is more/less likely for the error in $Y$ and/or $Z$ to be positive at 3 seconds? At 2 seconds? At 4 seconds?
– Bill
Apr 24, 2014 at 18:19
• A related by slightly different question is: how systematic is the measurement error? Suppose I said "Yeah, $X$ was measured a little high on my front lawn. $X$ is almost always measured a little high on my front lawn." Would that be a crazy statement? Does the measurement error work in such a way that a particular place might be very often too high while another particular place might be very often too low, etc" Or is all the error transitory?
– Bill
Apr 24, 2014 at 18:22
• @Bill: There is definitely serial correlation. The measurement errors are pretty much constant over the two hours. However, they are generally larger than the standard devtiation calculated from the data, which led me to this question. Apr 25, 2014 at 16:24
• Your question still does not clearly spell out the existence of serial correlation. Unfortunately, you have three carefully constructed answers not being nearly as useful to you as they might have been. Apr 27, 2014 at 9:13

I suspect that the previous responses to this question may be a bit off the mark. It seems to me that what the original poster is really asking here could be rephrased as, "given a series of vector measurements: $$\vec{\theta}_{i} = \left( \begin{array}{c} X_{i} \\ Y_{i} \\ Z_{i} \end{array}\right)$$ with $i=1, 2, 3,...,7200$, and measurement covariance: $$C_{i} = \left( \begin{array}{ccc} X_{\sigma,i}^{2} & 0 & 0 \\ 0 & Y_{\sigma,i}^{2} & 0 \\ 0 & 0 & Z_{\sigma,i}^{2} \end{array} \right)$$ how would I correctly calculate the covariance-weighted mean for this series of vector measurements, and afterward, how would I correctly calculate its standard deviation?" The answer to this question can be found in a lot of textbooks specializing in statistics for the physical sciences. One example that I like in particular is Frederick James, "Statistical Methods in Experimental Physics", 2nd edition, World Scientific, 2006, Section 11.5.2, "Combining independent estimates", pg. 323-324. Another very good, but more introductory-level text, which describes the variance-weighted mean calculation for scalar values (as opposed to full vector quantities as presented above) is Philip R. Bevington and D. Keith Robinson, "Data Reduction and Error Analysis for the Physical Sciences", 3rd edition, McGraw-Hill, 2003, Section 4.1.x, "Weighting the Data--Nonuniform Uncertainties". Because the poster's question happened to have a diagonalized covariance matrix in this case (i.e., all of the off-diagonal elements are zero), the problem is actually separable into three individual (i.e., X, Y, Z) scalar weighted mean problems, so the Bevington and Robinson analysis applies equally well here too.

In general, when responding to stackexchange.com questions, I don't normally find it useful to repackage long derivations that have already been presented before in numerous textbooks--if you want to truly understand the material, and understand why the answers look the way they do, then you really should just go and read the explanations which have already been published by the textbook authors. With that in mind, I'll simply jump directly to re-stating the answers that others have already provided. From Frederick James, setting $N=7200$, the weighted mean is: $$\vec{\theta}_{mean} = \left( \sum_{i=1}^{N} C_{i}^{-1} \right)^{-1} \left( \sum_{i=1}^{N} C_{i}^{-1} \vec{\theta_{i}} \right)$$ and the covariance of the weighted mean is: $$C_{mean} = \left( \sum_{i=1}^{N} C_{i}^{-1} \right)^{-1}$$ This answer is completely general, and will be valid no matter what the form of the $C_{i}$, even for non-diagonal measurement covariance matrices.

Since it so happens that the measurement covariances are diagonal in this particular case, the Bevington and Robinson analysis can also be used to calculate variance-weighted means for the individual $X_{i}$, $Y_{i}$, and $Z_{i}$. The form of the scalar answer is similar the form of the vector answer: $$X_{mean} = \frac{\sum_{i=1}^{N} \frac{X_{i}}{X_{\sigma,i}^{2}}}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}}$$ and the variance is $$X_{\sigma,mean}^{2} = \frac{1}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}}$$ or equivalently, $$X_{\sigma,mean} = \sqrt{\frac{1}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}}}$$ and similarly for $Y_{mean}, Y_{\sigma, mean}$ and $Z_{mean}, Z_{\sigma, mean}$. A brief wikipedia entry which also arrives at this same answer for the scalar-valued case is available here.

• Maybe I was a bit unclear, so I have added some more info. I don't think that I need to weight my measurements. Apr 24, 2014 at 16:34
• Yes you do. Consider an extreme case, just as a thought experiment: suppose you have only 2 GPS measurements, instead of 7200. Suppose furthermore that one of the GPS measurements has an uncertainty of +/- 5 feet, while the other has an uncertainty of +/- 5 miles. The uncertainty number literally tells you how potentially inaccurate the measurement is. That means the +/- 5 miles value is likely to be several miles off, at least. Do you really want to include this number in your average, in any meaningful way? Weighted averaging allows you to discount values that shouldn't be trusted as much. Apr 24, 2014 at 17:21
• BTW, my answer has another thing going for it: in your original post, you mention that the reason you don't want to simply use the sample standard deviation, calculated directly from the Z values, is that in that case, you would, in your own words, "neglect the fact that there are known measurement uncertainties". My answer (well, really, the obscure textbook answer, which I'm simply sharing with you) uses the known measurement uncertainties, exactly as you asked for. It's just that it uses the information in more places (mean result as well as standard deviation) than you had been expecting. Apr 24, 2014 at 18:33
• You convinced me. Apr 28, 2014 at 20:59

This should be easily solved using bayesian inference. You know the measurement properties of the individual points with respect to their true value and want to infer the population mean and SD that generated the true values. This is a hierarchical model.

# Rephrasing the problem (Bayes basics)

Note that whereas orthodox statistics give you a single mean, in the bayesian framework you get a distribution of credible values of the mean. E.g. the observations (1, 2, 3) with SDs (2, 2, 3) could have been generated by the Maximum Likelihood Estimate of 2 but also by a mean of 2.1 or 1.8, though slightly less likely (given the data) than the MLE. So in addition to the SD, we also infer the mean.

Another conceptual difference is that you have to define your knowledge state before making the observations. We call this priors. You might know in advance that a certain area was scanned and in a certain height range. The complete absence of knowledge would be to have uniform(-90, 90) degrees as the prior in X and Y and maybe uniform(0, 10000) meters on height (above the ocean, below the highest point on earth). You have to define priors distributions for all parameters that you want to estimate, i.e. get posterior distributions for. This is true for the standard deviation as well.

So rephrasing your problem, I assume that you want to infer credible values for three means (X.mean, Y.mean, X.mean) and three standard deviations (X.sd, Y.sd, X.sd) which could have generated your data.

# The model

Using standard BUGS syntax (use WinBUGS, OpenBUGS, JAGS, stan or other packages to run this), your model would look something like this:

  model {
# Set priors on population parameters
X.mean ~ dunif(-90, 90)
Y.mean ~ dunif(-90, 90)
Z.mean ~ dunif(0, 10000)
X.sd ~ dunif(0, 10)  # use something with better properties, i.e. Jeffreys prior.
Y.sd ~ dunif(0, 10)
Z.sd ~ dunif(0, 100)

# Loop through data (or: set up plates)
# assuming observed(x, sd(x), y, sd(y) z, sd(z)) = d[i, 1:6]
for(i in 1:n.obs) {
# The true value was generated from population parameters
X[i] ~ dnorm(X.mean, X.sd^-2)  #^-2 converts from SD to precision
Y[i] ~ dnorm(Y.mean, Y.sd^-2)
Z[i] ~ dnorm(Z.mean, Z.sd^-2)

# The observation was generated from the true value and a known measurement error
d[i, 1] ~ dnorm(X[i], d[i, 2]^-2)  #^-2 converts from SD to precision
d[i, 3] ~ dnorm(Y[i], d[i, 4]^-2)
d[i, 5] ~ dnorm(Z[i], d[i, 6]^-2)
}
}


Naturally, you monitor the .mean and .sd parameters and use their posteriors for inference.

# Simulation

I simulated some data like this:

# Simulate 500 data points
x = rnorm(500, -10, 5)  # mean -10, sd 5
y = rnorm(500, 20, 5)  # mean 20, sd 4
z = rnorm(500, 2000, 10)  # mean 2000, sd 10
d = cbind(x, 0.1, y, 0.1, z, 3)  # added constant measurement errors of 0.1 deg, 0.1 deg and 3 meters
n.obs = dim(d)


Then ran the model using JAGS for 2000 iterations after a burnin of 500 iterations. Here's the result for X.sd. The blue range indicates the 95% Highest Posterior Density or Credible interval (where you believe the parameter is after having observed the data. notice that an orthodox confidence interval does not give you this).

The red vertical line is the MLE estimate of the raw data. It is usually the case that the most likely parameter in Bayesian estimation is also the most likely (maximum likelihood) parameter in orthodox stats. But you should not care too much about the top of the posterior. The mean or median is better if you want to boil it down to a single number.

Notice that MLE/top is not at 5 because the data were randomly generated, not because of wrong stats.

# Limitiations

This is a simple model which has several flaws currently.

1. It doesn't handle the identity of -90 and 90 degrees. This can be done, however, by making some intermediate variable which shifts extreme values of estimated parameters into the (-90, 90) range.
2. X, Y and Z are currently modeled as independent though they are probably correlated and this should be taken into account to get the most out of the data. It depends on whether the measurement device was moving (serial correlation and joint distribution of X, Y and Z will give you a lot of information) or standing still (independence is ok). I can expand the answer to approach this, if requested.

I should mention that there's a lot of literature on spatial Bayesian models which I am not knowledgeable about.

• Thanks for this answer. It is data from a fixed station, but does this imply that the data are independent? Apr 22, 2014 at 7:30
• @traindriver You need to provide some more info about the inference problem you face in order for us to help you. You could expand your question with an "update" section specifying at least (1) is it the same quantity that is measured repeatedly? I.e. the same coordinate. Or is an area scanned or ... (2) why do you want to infer the mean and sd? If it's an area, it might be that you want to use SD as an estimate of bumpiness or something like that. Apr 22, 2014 at 12:49
• I have added some more info in the original post. Apr 24, 2014 at 16:31

I first introduce some notation and set up the problem using the simple approach you mentioned. Then go further. I will use $\mathbf{z}$ to refer to the vector Z you gave.

Consider the following model, which lacks the explicit mention measurement error: $\bar{Z} = \frac{\sum_{i=1}^n\mu_{Z} + \epsilon_i}{n}$, where $\bar{Z}$ is the estimated average value of $\mathbf{z}$, and $\mu_Z$ is the true average value of Z. Here, $\mathbf{\epsilon}$ is a vector of the errors in your data, and you expect that if your sample is large $\bar{Z}$ will converge to $\mu_Z$. If you simply take the observed $Z$ values and average them, you get $\bar{Z}$ and if you compute the sample standard deviation you get $\hat{\sigma}$, the estimate of the true population standard deviation $\sigma$. What if you would like to make use of some knowledge about the measurement error?

First, note that we can reformulate the initial model as: $\mathbf{z} = \mathbf{1}\beta + \epsilon$, where $\mathbf{1}$ is a vector of ones, and $\beta$ will end up being $\bar{Z}$. Now this really looks like regression, but we are still basically just getting an estimate of $\mu_Z$. If we perform a regression like this, we will also get an estimate for the standard error of $\epsilon$, which is almost what we want - this is nothing but the standard error of $\mathbf{z}$ (but we still want to account for measurement error).

We can augment our initial model to obtain a mixed effects model. $\mathbf{z} = \mathbf{1}\beta + \mathbf{Q}u +\epsilon$, where $u$ is a vector of random effects, and $\mathbf{Q}$ is the regressor relating $\mathbf{z}$ to $u$. As with any random effect, you will need to make an assumption about the distribution of $u$. Is it correct that $Z_\sigma$ is the distribution of the measurement error for $\mathbf{z}$? If yes, this can be used to provide the distribution of the random effects. Typically, software to perform basic mixed effects modeling will assume the random effects have a normal distribution (with mean 0...) and estimate the variance for you. Perhaps you can try this to test the concept. If you wish to use your prior information about the distribution of the measurement error, a Bayesian mixed effects model is in order. You can use R2OpenBUGS.

After estimating this model, the standard error you get for the residuals $\epsilon$ is the standard error you express interest in. Intuitively, the random effects component of the model is soaking up some of the variation that you can explain because you know there is measurement error. This allows you to obtain a more relevant estimate of the variation of $\epsilon$

See this paper for a deeper discussion on this approach of random effects to account for measurement error. Your situation is similar to the one the authors introduce for $\mathbf{D}$ and its measurement error corrupted version $\mathbf{W}$. The example in Section 4 may offer some insights into your situation.

As mentioned by whuber, you may wish to account for autocorrelation in your data. Using random effects will not solve that problem.