How to combine two measurements of the same quantity with different confidences in order to obtain a single value and confidence

Question

Back in the lab at university, we were taught to measure the quantity of interest some number of times (call this N), and then calculate the standard error. The underlying assumption here is that you are attempting to measure the true value of some quantity x, but your experimental errors are random and probably distributed Gaussianly about the true value (although I don't think the Gaussian assumption is necessary). Thus, standard statistical inference theory tells us that we infer the variance of this distribution from the sample variance (including the $\frac{1}{N-1}$ Bessel correction) and consequently, the standard error is given by $\frac{\sigma}{\sqrt{N}}$.

Furthermore, if one measures two different quantities (say x and y) and then wishes to combine then by some function $f(x,y)$ in order to obtain the actual quantity of interest, we take the mean values of x and y calculated and our estimator for $f(x,y)$ (f is evaluated with our mean values for x and y). In order to calculate the standard error on f, one uses the formula for propagation of errors: $\sigma _{f}=\sqrt{\left(\frac{\partial f}{\partial x}\sigma_{x}\right)^{2}+\left(\frac{\partial f}{\partial y}\sigma_{y}\right)^{2}}$ in which the partial derivatives are evaluated at the mean values of (x,y) and I have somewhat clumsily used $\sigma_{x}$ to denote the standard error on x.

A question never addressed, and the one I'm currently interested in, is what if I measure some quantity x and my friend measures the SAME quantity x. How do I combine these values into one best estimate of the true mean and assign a single standard error?

The key point here is that I'm assuming our measurements to no longer be coming from the same distribution, we are using different apparatus which has some fundamental bias in it, so our distributions will be centered on different values which, if you like, is our "apparatus' true mean", and it makes sense to average over the two to get a better estimate of the true mean.

A simplistic approach would be to take the mean of our two measurements. Then, the formula for propagation of errors says that the new standard error will be $\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}$. Furthermore, one might pragmatically take a weighted mean, so that the person whose standard error is lower has their measurement weighted more highly and again, we can evaluate the formula for standard error.

I find this approach dis-satisfactory however for the following reason. If I measure my quantity to be equal to $10 \pm 0.5$ and my friend gets $10 \pm 0.7$, intuitively the new estimate should be $10 \pm a$ in which $a<0.5$. If conversely, my friend had obtained $20 \pm 0.5$, intuitively my new best estimate is $15 \pm b$ with $b \gg 0.5$. This has a very Bayesian feel to it, like I should take the first measurement into account and then do a Bayesian update, but I'm struggling to connect the dots (/have run out of ideas for what I need to be Googling).

Any tips would be highly appreciated.

Answer 1

For exponential families, you should assume that the carrier measure is part of the prior, and so if you have two beliefs $X$ and $Y$, then combining them should not double-count the carrier measure. If their densities are $$f_X(x) \propto e^{\theta_X^T T(x) + h(x)}$$ and similarly for $Y$ using natural parameter $\theta_Y$, then by independence, their combination should have density $$f_{Z}(x) \propto e^{(\theta_X+\theta_Y)^T T(x) + h(x)}.$$

Therefore, for an exponential family, the Bayesian way of combining evidence is to add natural parameters. This is close to a pointwise product of densities (also called a product of experts), but doesn't make the mistake of double-counting the carrier measure. For normal distributions, the carrier measure is zero and so this distinction is irrelevant.

This “Bayesian evidence combination” does almost exactly what you describe except in the second case, $b<0.5$ since when you have more evidence, the variance has to go down:

If I measure my quantity to be equal to $10 \pm 0.5$ and my friend gets $10 \pm 0.7$, intuitively the new estimate should be $10 \pm a$ in which $a<0.5$. If conversely, my friend had obtained $20 \pm 0.5$, intuitively my new best estimate is $15 \pm b$ with $b \gg 0.5$.

If you work this out for a normal distribution in your chosen parametrization (mean and variance), you will find that the combined distribution's mean is the weighted average of means (weighted by precision) $$\frac{\frac{\mu_X}{\sigma^2_X} + \frac{\mu_Y}{\sigma^2_Y}}{\frac{1}{\sigma^2_X} + \frac{1}{\sigma^2_Y}},$$ and its variance is the “harmonic sum” of variances $$\frac{1}{\frac{1}{\sigma^2_X} + \frac{1}{\sigma^2_Y}}.$$

Answer 2

I don't know how to solve your question theoretically, but we can Monte Carlo it and check the results.

First, let's create some data. Here I assume measures have a normal shape where the error bars include 95% of the measurements:

set.seed(101)
measure1 <- rnorm(1e4,10,.25)  # 10 +-0.5
measure2 <- rnorm(1e4,10,.35)  # 10 +-0.7

library(coda)                  # just checking:
HPDinterval(as.mcmc(measure1)) #  95% HPD: 9.513229 10.48211
HPDinterval(as.mcmc(measure2)) #  95% HPD: 9.345024 10.69419

y <- c(measure1, measure2)

First let's see what a MLE fit would give us:

library(fitdistrplus)
fitdist(measures, "norm")
# Fitting of the distribution ' norm ' by maximum likelihood 
# Parameters:
#        estimate  Std. Error
# mean 10.0013619 0.002135679
# sd    0.3020306 0.001510079

Herein, I specified the simple BUGS model:

model {
  for (i in 1:n) {
    y[i] ~ dnorm(mu, tau)
  }

  mu  ~ dunif(9, 11)      #      mean prior
  tau ~ dgamma(0.01,0.01) # precision prior
  sigma <- 1 / sqrt(tau)
}

The result is:

         mean       sd  MC_error val2.5pc  median val97.5pc start sample
mu    10.0000 0.002125 1.652e-05   9.9970 10.0000    10.010  1501  15000
sigma  0.3021 0.001507 1.137e-05   0.2991  0.3021     0.305  1501  15000

which is quite similar to the MLE fit.

Considering the case where the 2nd measure was around 20, and still assuming the model of a single unimodal distribution:

measure3 <- rnorm(1e4,20,.25) 
HPDinterval(as.mcmc(measure3)) # 95% HPD: 19.48301 20.46369
y2       <- c(measure1, measure3)

fitdist(y2, "norm")
# Fitting of the distribution ' norm ' by maximum likelihood 
# Parameters:
#       estimate Std. Error
# mean 15.001414 0.03539963
# sd    5.006263 0.02503131

The BUGS model is adapted to allow a wider prior for $\mu$

  model {
    for (i in 1:n) {
      y[i] ~ dnorm(mu, tau)
    }

    mu  ~ dunif(9, 21)      #      mean prior
    tau ~ dgamma(0.01,0.01) # precision prior
    sigma <- 1 / sqrt(tau)
  }

Result:

        mean      sd  MC_error val2.5pc median val97.5pc start sample
mu    15.000 0.03522 0.0002739   14.930 15.000    15.070  1501  15000
sigma  5.007 0.02498 0.0001884    4.957  5.007     5.056  1501  15000

Again quite similar to the MLE fit.

This second case, imho, should not be modeled like this. Even assuming that both measurements were ok (no strange measurement problem occurred), it would be more reasonable to assume a model with a bimodal distribution and, by some reason, both experiences censored the other 'half' of the measurement space. This could also be modeled by BUGS.

Btw, for arbitrary arithmetic expressions dealing with error measurements, there is a nice online tool Guesstimate which functions as a spreadsheet (it uses Monte Carlo to produce the resulting error bars):

Not sure if I helped you :-)