Mean or average of two partially correlated measurements

Fra mn/fys/epf
Revisjon per 2. des. 2016 kl. 09:02 av Read@uio.no (diskusjon | bidrag) (The derivation of the terms in the correlation matrix)

(diff) ← Eldre revisjon | Nåværende revisjon (diff) | Nyere revisjon → (diff)
Hopp til: navigasjon, søk

See this article for background material or a similar article with an example application from particle physics.

Use the root macro to perform an average of two measurements:
   x1 +- dx1s (statistical) +- dx1u (uncorrelated systematic) +- dx1c (correlated systematic)
   x2 +- dx2s (statistical) +- dx2u (uncorrelated systematic) +- dx2c (correlated systematic)

resulting in mean m +- dmstat +- dmsyst. The systematic errors in each channel are decomposed in uncorrelated and (100%) correlated components. The correlation matrix which appears in the solution is composed of the sum of a diagonal covariance matrix with elements given by the uncorrelated uncertainties summed in quadrature

dx1s2 + dx1u2
0
0
dx2s2 + dx2u2

and a non-diagonal covariance matrix for the correlated uncertainties with correlation coefficient rho=1

dx1c2
dx1c * dx2c
dx1c * dx2c
dx2c2

. From the form of the covariance matrix C

dx12
rho * dx1 * dx2
rho * dx1 * dx2
dx22

we can identify dxi2 = dxis2 + dxiu2 + dxic2 and rho = dx1c * dx2c / (dx1 * dx2).

Minimizing the generalized chi-squared XTC-1X, where X is a column vector

x1 - m
x2 - m

, we get for the minimum variance estimate of m

m = (x1/dx12 + x2/dx22 - rho * (x1 + x2) /(dx1 * dx2)) / (1/dx12 + 1/dx22 - 2 * rho/(dx1 * dx2))

and the variance of m

dm2 = (1-rho2) / ( 1/dx12 + 1/dx22 - 2 * rho/(dx1 * dx2)). We decompose the variance into statistical and systematic components by subtraction in quadrature of the statistical uncertainty dmstat2 = 1/(1/dx1s2 + 1/dx2s2), dmsyst2 = dm2 - dmstat 2.


There is protection in the code against the very special case that there is only a 100% correlated uncertainty. If this is truely the case then the 2 measurements must be equal by construction and the uncertainty may be taken as the smaller of the 2.


There is also code to show the results of a popular approximation that does not properly take into account the correlation between the 2 measurements.


Here is a sample output from the macro for some (almost) randomly chosen measurement results:

*******************************************
* *
* W E L C O M E to R O O T *
* *
* Version 5.28/00 14 December 2010 *
* *
* You are welcome to visit our Web site *
* http://root.cern.ch *
* *
*******************************************

ROOT 5.28/00 (trunk@37585, Dec 14 2010, 15:20:27 on linuxx8664gcc)

CINT/ROOT C/C++ Interpreter version 5.18.00, July 2, 2010
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
root [0]
Processing AverageMeasurements.C...
The measurements being averaged:
-------------------------------
x1 = 58.9 +- 3.4 (stat) +- 1.5 (uncorr syst) +- 2.4 (corr syst)
   = 58.9 +- 4.4238 (total)

x2 = 68.7 +- 2.8 (stat) +- 0.3 (uncorr syst) +- 3.9 (corr syst)
   = 68.7 +- 4.81041 (total)

Results for the generalized weighted average
--------------------------------------------
Correlation coefficient (rho) = 0.439844

m = 63.0708 +- 2.1614 (stat) +- 3.24854 (syst)
  = 63.0708 +- 3.90188 (total)

Generalized chi-squared = 4.00333

Approximate, simple formulae
----------------------------
m = 65.1253 +- 2.1614 (stat) +- 3.20753 (syst)
  = 65.1253 +- 3.8678 (total)

Generalized chi-squared = 4.28057 (for the approximate minimum)
Delta chi-squared with respect to exact minimum = 0.277241

The derivation of the terms in the correlation matrix

The covariance matrix is
[math]C_{ij} = \lt (x_i-\bar{x_i})(x_j-\bar{x_j})\gt [/math].
Suppose [math]\delta_1,\ \delta_2,\ \delta_3[/math] are three independent sources of normally-distributed unit fluctuations (with < δi > = 0 and < δi * δj > = δij and where [math]\delta_{ij}[/math] is the kroneker delta function (1 for [math]i=j[/math] and zero for [math]i\neq j[/math]).


Pseudo-measurements of (x1,x2) can be generated from the expressions (x1,x2) = (x10 + α * δ1 + β * δ2,x20 + γ * δ3 + λ * δ2).

Expanding the covariance matrix one finds [math]C_{ij} = \lt {x_i*x_j}\gt -\bar{x_i}\bar{x_j}[/math].

Using the properties of [math]\lt \delta_i*\delta_j\gt [/math] above and [math]\lt {x_i}\gt =\bar{x_i}[/math], one finds [math]\lt {x_1*x_2}\gt = x_{10}*x_{20}+\beta*\lambda[/math]; so that C12 = C21 = β * λ. Similar substitution gives for the diagonal matrix elements
[math]C_{11}=x_{10}^2+\alpha^2+\beta^2-x_{10}^2[/math]
and
[math]C_{22}=x_{20}^2+\gamma^2+\lambda^2-x_{20}^2[/math]
so that C=

α2 + β2 β * λ
[math]\beta*\lambda[/math] [math]\gamma^2+\lambda^2[/math]