How do I propagate the error in the elements of a set through to the mean of the elements?

Question

I need to understand how to propagate the errors in the elements of a set through to the mean of the set and correctly determine the estimate for the error value of the mean.

I am designing an experiment to determine the time response of a sensor. My experiment will repeatably measure the time response of the same sensor and will generate a set of independent measurements of the time response using the ramped input method. Each element in the set will have the same random error which is RSS of the precision of the input to and output of the sensor expressed as a standard deviation which has been predetermined in another experiment $\bigl( \sigma_\epsilon \bigr) $. My thoughts are as follows:

Due to the error in the elements of the time data there will be an error in the calculation of the mean. Applying the standard error propagation formula to the equation for the mean: $$ \frac{\delta z}{z}\bigl( w,x,y,... \bigr)= \sqrt{\bigl( {\frac{\partial z}{\partial w}\delta w} \bigr)^2+\bigl( {\frac{\partial z}{\partial x}\delta x} \bigr)^2+\bigl( {\frac{\partial z}{\partial y}\delta y} \bigr)^2+...}$$

$$ \mu =\frac 1N \sum_{i=1}^N x_i $$ Yields the result: $$ \delta \mu =\frac 1N \sqrt {\sum_{i=1}^N \sigma _i^2} $$ As $ \sigma_i $ is constant for all i and normally distributed, the equation for the error in the mean simplifies to: $$ \sigma_{\mu} =\frac{\sigma_\epsilon}{\sqrt N}$$ The current estimate for the mean is therefore: $$ \mu =\frac 1N \sum_{i=1}^N x_i \pm \frac{n\sigma_\epsilon}{\sqrt N}$$ Where n is the number of standard deviations required.

In my case n = 3, so there is a 99.79% chance that the true mean is within $\pm $ three standard deviations of the estimated mean.

Next we need to consider the scatter of the elements around the estimated mean. This is the sample standard deviation of the elements ($\sigma_x$). So we now have two normally distributed independent contributions to the error; the error in the mean due to $ \sigma_\epsilon $, and the scatter around the mean $ \sigma_x $. The RSS of these values should be used and therefore the estimate for the time response would be: $$ \tau = \bar x \pm n\sqrt {\bigl( \sigma_\mu \bigr)^2 + \bigl( \sigma_x \bigr)^2} $$ I realise there is an error in $\sigma_x$ due to the error in the mean but this just adds an error to an error and I assume we can ignore this because it is symmetric?

This is probably all completely wrong and I would really appreciate someone correcting me.

Thanks.

Answer 1

Using physics notation, when adding or subtracting measurements $x,\ldots,w$ which have uncertainties $\delta x, \ldots, \delta w$, and the values are used to compute

$q=x + \cdots + z - (u + \cdots + w)$,

then the uncertainty in $q$ is

$\delta q = \sqrt{(\delta x)^2 + \cdots + (\delta z)^2 + (\delta u)^2 + \cdots + (\delta w)^2}$.

In statistical notation this is:

$\sigma_q = \sqrt{\sigma_x^2 + \cdots + \sigma_z^2 + \sigma_u^2 + \cdots + \sigma_w^2}$ .

Answer 2

The easiest way to think about this is as a chain of two Normal random variables.

You've said your sensor has a known error variance (or at least estimated from a previous experiment) and known no bias (implied it's symmetric about zero, I'll show you how to incorporate non-zero bias either way though). So you have basically defined the distribution of errors as $$e_i \sim N(b, \Delta^2)$$ where $b$ is the known bias of the sensor and $\Delta$ is the known standard deviation of the errors.

If you think of your measurements using this sensor as Normal, then, assuming you know the error for each observation $e_i$, your measurements are simply $$x_i | e_i \sim N(\mu + e_i, \sigma^2).$$ If we remember that the joint distribution is just the product of the conditional distribution and prior then $P(x_i, e_i, \mu) = P(x_i | e_i) P(e_i)$. We can get the marginal distribution for the $x_i$ variables by integrating out $e_i$. Using a common result about Normal distributions (called conjugacy if you want to look it up), we can see $$x_i \sim N(\mu + b, \sigma^2 + \Delta^2)$$.

The maximum likelihood estimate of $\mu$ and $\sigma^2$ are derived from this distribution using standard techniques. In short, $$\hat{\mu} = \bar{x}-b$$ and $$\hat{\sigma^2} = \frac{1}{N}\sum_i (x_i - \bar{x})^2 - \Delta^2.$$

So, if you're interested in the interval for $\mu$ (or $\tau$ in your derivation, I believe), then you just apply these MLEs and the results about averages of normal random variables to get $$ \bar{x} \sim N(\mu, \frac{\sigma^2 + \Delta^2}{N}) $$ and apply this to the typical Normal interval calculation: $$\mu = \bar{x} - b \pm z_\alpha \sqrt{\frac{\hat{\sigma^2} + \Delta^2}{N}}$$ which is simply $$\mu = \bar{x} - b \pm z_\alpha \sqrt{\frac{\frac{1}{N}\sum_i (x_i - \bar{x})^2 - \Delta^2 + \Delta^2}{N}}$$ which simplifies easily enough to $$\mu = \bar{x} - b \pm z_\alpha \sqrt{\frac{\frac{1}{N}\sum_i (x_i - \bar{x})^2}{N}}$$

In your case, you've stated $b = 0$, so that's easy enough to replace!