Combining Experiments with Systematic Errors

Roger John Barlow

arXiv:1701.03701·physics.data-an·September 29, 2020

Combining Experiments with Systematic Errors

Roger John Barlow

PDF

TL;DR

This paper analyzes methods for combining multiple datasets with shared systematic uncertainties, demonstrating their equivalence or differences depending on error type and application, which is crucial for accurate data analysis at particle colliders.

Contribution

It clarifies the theoretical relationship between the full matrix and extra parameter techniques for systematic error treatment in data fitting, including conditions for their equivalence or bias.

Findings

01

For additive uncertainties, the methods are theoretically equivalent.

02

For multiplicative errors, the matrix fit can be biased if applied to data points.

03

Applying the multiplicative factor to the function avoids bias in the fit.

Abstract

We consider fits to two or more datasets for which results from the sa me experiment share a common systematic uncertainty in addition to their individ ual statistical errors. This is important in extracting the maximum information from a set of similar bu t different experiments (or the same experiment under different conditions) an alysing similar but different datasets, as happens at the LHC and other particle colliders. There are two techniques in use: using the full matrix and using extra paramneters, and we show, for a completely general fit, that for an addit ive uncertainty they are in principle equivalent even though in practice the det ails differ. For a multiplicative error the matrix fit is equivalent to the extra parameter fit if the factor is applied to the data points but not if it is applied to the function: the former leads to biassed estimates and the latter avoids…

Equations44

χ^{2} = i \sum (\frac{y _{i} - f ( x _{i} ; a )}{σ _{i}})^{2}

χ^{2} = i \sum (\frac{y _{i} - f ( x _{i} ; a )}{σ _{i}})^{2}

χ^{2} = i, j \sum (y_{i} - f (x_{i}; a)) V_{ij}^{- 1} (y_{j} - f (x_{j}; a))

χ^{2} = i, j \sum (y_{i} - f (x_{i}; a)) V_{ij}^{- 1} (y_{j} - f (x_{j}; a))

χ^{2} = r \sum i \in E_{r} \sum (\frac{y _{i} - f ( x _{i} ; a ) - z _{r}}{σ _{i}})^{2} + r \sum (\frac{z _{r}}{S _{r}})^{2}

χ^{2} = r \sum i \in E_{r} \sum (\frac{y _{i} - f ( x _{i} ; a ) - z _{r}}{σ _{i}})^{2} + r \sum (\frac{z _{r}}{S _{r}})^{2}

χ^{2} = r \sum i \in E_{r} \sum (\frac{y _{i} - z _{r} f ( x _{i} ; a )}{σ _{i}})^{2} + r \sum (\frac{z _{r} - 1}{S _{r}})^{2}

χ^{2} = r \sum i \in E_{r} \sum (\frac{y _{i} - z _{r} f ( x _{i} ; a )}{σ _{i}})^{2} + r \sum (\frac{z _{r} - 1}{S _{r}})^{2}

V_{r}=\left(\matrix{\sigma_{1}^{2}+S_{r}^{2}&S_{r}^{2}&S_{r}^{2}&\dots\cr S_{r}^{2}&\sigma_{2}^{2}+S_{r}^{2}&S_{r}^{2}&\dots\cr S_{r}^{2}&S_{r}^{2}&\sigma_{3}^{2}+S_{r}^{2}&\dots\cr\vdots&\vdots&\vdots\cr}\right)

V_{r}=\left(\matrix{\sigma_{1}^{2}+S_{r}^{2}&S_{r}^{2}&S_{r}^{2}&\dots\cr S_{r}^{2}&\sigma_{2}^{2}+S_{r}^{2}&S_{r}^{2}&\dots\cr S_{r}^{2}&S_{r}^{2}&\sigma_{3}^{2}+S_{r}^{2}&\dots\cr\vdots&\vdots&\vdots\cr}\right)

χ^{2} = r \sum i \in E_{r} \sum j \in E_{r} \sum Δ_{i} (V_{r}^{- 1})_{ij} Δ_{j}

χ^{2} = r \sum i \in E_{r} \sum j \in E_{r} \sum Δ_{i} (V_{r}^{- 1})_{ij} Δ_{j}

V_{r}^{- 1} = \frac{1}{σ ^{2}} I - (\frac{S _{r}^{2}}{σ ^{2} ( σ ^{2} + n _{r} S _{r}^{2} )}) U

V_{r}^{- 1} = \frac{1}{σ ^{2}} I - (\frac{S _{r}^{2}}{σ ^{2} ( σ ^{2} + n _{r} S _{r}^{2} )}) U

(V_{r}^{- 1})_{ij} = \frac{δ _{ij}}{σ _{i}^{2}} - \frac{S _{r}^{2}}{1 + \sum _{k \in E_{r}} \frac{S _{r}^{2}}{σ _{k}^{2}}} \frac{1}{σ _{i}^{2} σ _{j}^{2}}

(V_{r}^{- 1})_{ij} = \frac{δ _{ij}}{σ _{i}^{2}} - \frac{S _{r}^{2}}{1 + \sum _{k \in E_{r}} \frac{S _{r}^{2}}{σ _{k}^{2}}} \frac{1}{σ _{i}^{2} σ _{j}^{2}}

χ^{2} = r \sum i \in E_{r} \sum \frac{Δ _{i}^{2}}{σ _{i}^{2}} - r \sum \frac{S _{r}^{2}}{1 + \sum _{k \in E_{r}} \frac{S _{r}^{2}}{σ _{k}^{2}}} i \in E_{r} \sum j \in E_{r} \sum \frac{Δ _{i} Δ _{j}}{σ _{i}^{2} σ _{j}^{2}} .

χ^{2} = r \sum i \in E_{r} \sum \frac{Δ _{i}^{2}}{σ _{i}^{2}} - r \sum \frac{S _{r}^{2}}{1 + \sum _{k \in E_{r}} \frac{S _{r}^{2}}{σ _{k}^{2}}} i \in E_{r} \sum j \in E_{r} \sum \frac{Δ _{i} Δ _{j}}{σ _{i}^{2} σ _{j}^{2}} .

\frac{\partial χ ^{2}}{\partial z _{r}} = - 2 j \sum \frac{Δ _{j} - z _{r}}{σ _{j}^{2}} + 2 \frac{z _{r}}{S _{r}^{2}}

\frac{\partial χ ^{2}}{\partial z _{r}} = - 2 j \sum \frac{Δ _{j} - z _{r}}{σ _{j}^{2}} + 2 \frac{z _{r}}{S _{r}^{2}}

\overset{z}{^}_{r} = ζ_{r} j \sum \frac{Δ _{j}}{σ _{j}^{2}} .

\overset{z}{^}_{r} = ζ_{r} j \sum \frac{Δ _{j}}{σ _{j}^{2}} .

\eqalign{\chi^{2}=\sum_{r}\Biggl{(}\Biggr{.}&\sum_{i\in E_{r}}{\Delta_{i}^{2}\over\sigma_{i}^{2}}-2\sum_{i\in E_{r}}{\Delta_{i}\over\sigma_{i}^{2}}\zeta_{r}\sum_{j\in E_{r}}{\Delta_{j}\over\sigma_{j}^{2}}\cr&+\left(\sum_{k\in E_{r}}{1\over\sigma_{k}^{2}}+{1\over S^{2}}\right)\zeta_{r}^{2}\sum_{i\in E_{r}}\sum_{j\in E_{r}}{\Delta_{i}\over\sigma_{i}^{2}}{\Delta_{j}\over\sigma_{j}^{2}}\Biggl{.}\Biggr{)}}

\eqalign{\chi^{2}=\sum_{r}\Biggl{(}\Biggr{.}&\sum_{i\in E_{r}}{\Delta_{i}^{2}\over\sigma_{i}^{2}}-2\sum_{i\in E_{r}}{\Delta_{i}\over\sigma_{i}^{2}}\zeta_{r}\sum_{j\in E_{r}}{\Delta_{j}\over\sigma_{j}^{2}}\cr&+\left(\sum_{k\in E_{r}}{1\over\sigma_{k}^{2}}+{1\over S^{2}}\right)\zeta_{r}^{2}\sum_{i\in E_{r}}\sum_{j\in E_{r}}{\Delta_{i}\over\sigma_{i}^{2}}{\Delta_{j}\over\sigma_{j}^{2}}\Biggl{.}\Biggr{)}}

(V_{r})_{ij}^{- 1} = \frac{δ _{ij}}{σ _{i}^{2}} - \frac{ξ _{r}^{2} y _{i} y _{j}}{1 + ξ _{r}^{2} \sum _{k \in E_{r}} \frac{y _{k}^{2}}{σ _{k}^{2}}} (\frac{1}{σ _{i}^{2} σ _{j}^{2}})

(V_{r})_{ij}^{- 1} = \frac{δ _{ij}}{σ _{i}^{2}} - \frac{ξ _{r}^{2} y _{i} y _{j}}{1 + ξ _{r}^{2} \sum _{k \in E_{r}} \frac{y _{k}^{2}}{σ _{k}^{2}}} (\frac{1}{σ _{i}^{2} σ _{j}^{2}})

χ^{2} = r \sum i \in E_{r} \sum (\frac{z _{r} y _{i} - f _{i}}{σ _{i}})^{2} + r \sum (\frac{z _{r} - 1}{ξ _{r}})^{2}

χ^{2} = r \sum i \in E_{r} \sum (\frac{z _{r} y _{i} - f _{i}}{σ _{i}})^{2} + r \sum (\frac{z _{r} - 1}{ξ _{r}})^{2}

χ^{2} = r \sum i \in E_{r} \sum (\frac{y _{i} - z _{r} f _{i}}{σ _{i}})^{2} + r \sum (\frac{z _{r} - 1}{ξ _{r}})^{2}

χ^{2} = r \sum i \in E_{r} \sum (\frac{y _{i} - z _{r} f _{i}}{σ _{i}})^{2} + r \sum (\frac{z _{r} - 1}{ξ _{r}})^{2}

\overset{z}{^}_{r} = \frac{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{f _{i} y _{i}}{σ _{i}^{2}}}{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{y _{i}^{2}}{σ _{i}^{2}}}

\overset{z}{^}_{r} = \frac{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{f _{i} y _{i}}{σ _{i}^{2}}}{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{y _{i}^{2}}{σ _{i}^{2}}}

χ^{2} = i \sum \frac{Δ _{i}^{2}}{σ _{i}^{2}} - r \sum \frac{ξ _{r}^{2}}{1 + ξ _{r}^{2} \sum _{k \in E_{r}} \frac{y _{k}^{2}}{σ _{k}^{2}}} (i \in E_{r} \sum \frac{Δ _{i} y _{i}}{σ _{i}^{2}})^{2}

χ^{2} = i \sum \frac{Δ _{i}^{2}}{σ _{i}^{2}} - r \sum \frac{ξ _{r}^{2}}{1 + ξ _{r}^{2} \sum _{k \in E_{r}} \frac{y _{k}^{2}}{σ _{k}^{2}}} (i \in E_{r} \sum \frac{Δ _{i} y _{i}}{σ _{i}^{2}})^{2}

\overset{z}{^}_{r} = \frac{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{f _{i} y _{i}}{σ _{i}^{2}}}{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{f _{i}^{2}}{σ _{i}^{2}}}

\overset{z}{^}_{r} = \frac{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{f _{i} y _{i}}{σ _{i}^{2}}}{1 + ξ _{r}^{2} \sum _{i \in E_{r}} \frac{f _{i}^{2}}{σ _{i}^{2}}}

χ^{2} = i \sum \frac{Δ _{i}^{2}}{σ _{i}^{2}} - r \sum \frac{ξ _{r}^{2}}{1 + ξ _{r}^{2} \sum _{k \in E_{r}} \frac{f _{k}^{2}}{σ _{k}^{2}}} (i \in E_{r} \sum \frac{Δ _{i} f _{i}}{σ _{i}^{2}})^{2}

χ^{2} = i \sum \frac{Δ _{i}^{2}}{σ _{i}^{2}} - r \sum \frac{ξ _{r}^{2}}{1 + ξ _{r}^{2} \sum _{k \in E_{r}} \frac{f _{k}^{2}}{σ _{k}^{2}}} (i \in E_{r} \sum \frac{Δ _{i} f _{i}}{σ _{i}^{2}})^{2}

\chi^{2}=(x_{1}-\hat{x},x_{2}-\hat{x})\left(\matrix{x_{1}^{2}\xi^{2}+\sigma_{1}^{2}&x_{1}x_{2}\xi^{2}\cr x_{1}x_{2}\xi^{2}&x_{2}^{2}\xi^{2}+\sigma_{2}^{2})\cr}\right)^{-1}\left(\matrix{x_{1}-\hat{x}\cr x_{2}-\hat{x}}\right)

\chi^{2}=(x_{1}-\hat{x},x_{2}-\hat{x})\left(\matrix{x_{1}^{2}\xi^{2}+\sigma_{1}^{2}&x_{1}x_{2}\xi^{2}\cr x_{1}x_{2}\xi^{2}&x_{2}^{2}\xi^{2}+\sigma_{2}^{2})\cr}\right)^{-1}\left(\matrix{x_{1}-\hat{x}\cr x_{2}-\hat{x}}\right)

\overset{x}{^} = \frac{x _{1} σ _{2}^{2} + x _{2} σ _{1}^{2}}{σ _{1}^{2} + σ _{2}^{2} + ( x _{1} - x _{2} ) ^{2} ξ ^{2}}

\overset{x}{^} = \frac{x _{1} σ _{2}^{2} + x _{2} σ _{1}^{2}}{σ _{1}^{2} + σ _{2}^{2} + ( x _{1} - x _{2} ) ^{2} ξ ^{2}}

\chi^{2}=(x_{1}-\hat{x},x_{2}-\hat{x})\left(\matrix{\hat{x}^{2}\xi^{2}+\sigma_{1}^{2}&\hat{x}^{2}\xi^{2}\cr\hat{x}^{2}\xi^{2}&\hat{x}^{2}\xi^{2}+\sigma_{2}^{2})\cr}\right)^{-1}\left(\matrix{x_{1}-\hat{x}\cr x_{2}-\hat{x}}\right)

\chi^{2}=(x_{1}-\hat{x},x_{2}-\hat{x})\left(\matrix{\hat{x}^{2}\xi^{2}+\sigma_{1}^{2}&\hat{x}^{2}\xi^{2}\cr\hat{x}^{2}\xi^{2}&\hat{x}^{2}\xi^{2}+\sigma_{2}^{2})\cr}\right)^{-1}\left(\matrix{x_{1}-\hat{x}\cr x_{2}-\hat{x}}\right)

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Combining Experiments with Systematic Errors

Roger Barlow

The University of Huddersfield, Huddersfield HD1 3DH, UK

$26^{th}$ * September 2020*

Abstract

We consider fits to two or more datasets for which results from the same experiment share a common systematic uncertainty in addition to their individual statistical errors. This is important in extracting the maximum information from a set of similar but different experiments (or the same experiment under different conditions) analysing similar but different datasets, as happens at the LHC and other particle colliders. There are two techniques in use: using the full matrix and using extra paramneters, and we show, for a completely general fit, that for an additive uncertainty they are in principle equivalent even though in practice the details differ. For a multiplicative error the matrix fit is equivalent to the extra parameter fit if the factor is applied to the data points but not if it is applied to the function: the former leads to biassed estimates and the latter avoids them.

1 The problem

It is frequently required to fit the parameter(s) $a$ of a function $f(x;a)$ to values $\{x_{i},y_{i}\}$ where the measurements come from several different experiments; a measurement $i$ in experiment $E_{r}$ is subject not only to an independent error $\sigma_{i}$ but to a systematic uncertainty $S_{r}$ which is shared with all other measurements of that experiment. This uncertainty may be additive (for example, in a common background), or multiplicative (for example, in a common efficiency).

For example the $e^{+}e^{-}$ hadronic cross section ratio $R$ has been measured by many experiments over a wide range of energies and fitting parameters, for example the $Z$ width [1], is best done by a global fit to the different experiments. Many determinations of neutrino mixing parameters, for example [2], involve global fits to several different neutrino experiments. Describing proton-proton elastic and diffractive scattering using Pomeron exchange models [3] fits parameters using data from many diffierent experiments at the ISR, the Tevatron and the LHC with each experiment sharing a systematic uncertainty.

Figure 1: Data from two different experiments

An illustration is shown in Figure 1, where straight-line-fit data from two different experiments is depicted by circles and squares. There is a clear systematic shift between the two datasets. Each on its own is well described by a straight line, but if the two are combined without regard for the systematic the dashed line is obtained, which is clearly wrong. The dotted line is produced by the techniques discussed below.

The problem is not that there is no obvious way to perform such a fit: the problem is that there are two such methods: by full matrix inversion and by introducing additional parameters.

Suppose that there are $n$ experiments, $E_{1}\dots E_{r}\dots E_{n}$ . (In what follows we use the index $r$ to denote the experiment, and indices $i,j,k$ to denote individual data points.) For uncorrelated data points one determines the value(s) of $a$ by minimising the $\chi^{2}$

[TABLE]

For correlated data the $1/\sigma^{2}$ terms are replaced by the inverse of the error matrix

[TABLE]

and the value(s) of $a$ are determined by minimising this.

However alternatively one can introduce parameters $z_{r}$ to describe the shared deviations of each experiment. For an additive uncertainty $z_{r}$ has mean zero and standard deviation $S_{r}$ , and the $\chi^{2}$ is written as

[TABLE]

If the uncertainty is multiplicative then the $\chi^{2}$ can be written

[TABLE]

where the $z_{r}$ scale factors have mean 1.0.

We shall show that these apparently different fitting procedures lead to the same result.

2 Past Approaches

D’Agostini [4] discusses the structure of the covariance matrix for muliplicative systematic uncertainties. In considering the special case of forming the average of two quantities (which corresponds to $f(x;a)=a)$ ) they assert that the methods give the same result, but do not generalise this statement. Demortier, in an unpublished note [5], shows the methods are equivalent for additive uncertainties, giving quite a complicated proof. Mo [6] considers multplicative uncertainties but only considers linear fits, for which they assert that the methods are equivalent (they refer to a paper which is unfortunately only in Chinese.) Blobel [7] asserts that the two approaches are equivalent “for simple problems” and “for additive uncertainties” but does not offer a proof. Fogli et al.[8] give a proof of equivalence for the additive case.

3 Additive Systematic Errors

The inversion of the matrix $V$ can be performed algebraically. First, note that if the data are grouped together according to their experiment, $V$ is block-diagonal, each block having the form

[TABLE]

so the inverse of $V$ is also block diagonal in the individual inverse matrices $V_{r}^{-1}$ , and Equation (2) has the form

[TABLE]

where the sums over $i$ and $j$ run over measurements of the experiment $r$ and $\Delta_{i}=y_{i}-f(x_{i};a)$ .

The individual submatrices can also be inverted. To see this, we first consider the case where all $\sigma_{i}$ are the same. Introduce the matrix $U$ for which all members are unity, and note that for an $n\times n$ square matrix, $U^{2}=nU$ . The submatrix may be written as $V_{r}=\sigma^{2}I+S_{r}^{2}U$ , and $V^{-1}_{r}$ is just:

[TABLE]

If the $\sigma_{i}$ are different then the situation is not quite so elegant, but the inverse is still writeable in closed form as

[TABLE]

as can be seen by multiplying out this expression and $(V_{r})_{ij}=\delta_{ij}\sigma_{i}^{2}+S_{r}^{2}$ .

The $\chi^{2}$ to be minimised in the first method is thus

[TABLE]

For the second method the expression for $\chi^{2}$ is given by Equation 3. To find the minimum the differentials are set to to zero, and the resulting system of simultaneous equations partially decouples. The differential for one of the $z_{r}$ is

[TABLE]

which gives the estimate of $z_{r}$ for a particular $a$ , writing, for convenience, $\zeta_{r}={1\over{1\over S_{r}^{2}}+\sum_{k}{1\over\sigma_{k}^{2}}}$

[TABLE]

Inserting Equation 11 into Equation 3 gives

[TABLE]

which reduces to the expression in Equation 9.

The two methods are thus equivalent, minimising the same $\chi^{2}$ . (This is used to obtain the dotted line in Figure 1.) From the mathematical point of view this may be surprising as the formulæ are very different; it is gratifying from the statistical viewpoint as both approaches make the same assumptions.

Either method can therefore be chosen in a particular problem. The matrix method is simpler, particularly if the model $f(x;a)$ is linear in $a$ as the normal equations can be solved directly, without the need for iteration. However the second method gives the $z_{r}$ offsets as a by-product, and these may be useful to investigate the quality of the different datasets.

The second method increases the difficulty of the minimisation by more than just having extra parameters to fit. If the function includes an additive constant term (which it generally will), the search space includes a direction in which a change in this term can be balanced by an opposite shift in the mean $z$ . In this direction $\chi^{2}$ changes slowly, through the $1/S^{2}$ terms, whereas the change in any orthogonal direction is, assuming that there are many more datapoints than datasets, much larger. This causes a problem of the Rosenbrock’s Valley type. Even if the minimiser reports ‘successful convergence’, the result should be confirmed with care. It may be preferable, rather than setting $a$ and all the $z_{r}$ as adjustable parameters in the $\chi^{2}$ of Equation 3, to work solely using the $a$ parameters, with the $z_{r}$ at every step being given explicitly by Equation 11.

4 Multiplicative Systematic Errors

Correlated multiplicative errors are more common than correlated additive errors, as uncertainties in acceptance, efficiency and overall normalisation are of this type.

Again one can use a matrix method or introduce additional parameters. The submatrix for individual experiment $E_{r}$ which has a fractional uncertainty $\xi_{r}$ is $(V_{r})_{ij}=\delta_{ij}\sigma_{i}^{2}+\xi^{2}_{r}y_{i}y_{j}$ , for which the block diagonal inverse, equivalent to Equation 8, is

[TABLE]

.

For the separate fit, there is an ambiguity about whether to apply the error to the value or the function. In comparing some measured $y_{i}$ values with function values $f(x_{i},a)$ with a nominal fractional error $\xi$ one can apply scale factors $z_{r}$ and write the $\chi^{2}$ as

[TABLE]

or as

[TABLE]

Leaving aside for now the arguments for the choice - which amounts to the difference between whether a quantity or its reciprocal is better described by a Gaussian, a question that does not arise for an additive correction - one can estimate the $z_{r}$ for a given $a$ in the first case by setting the differential of Equation 14 to zero

[TABLE]

and if these are inserted into Equation 14 then this gives

[TABLE]

which exactly matches the $\chi^{2}$ given by the matrix elements in Equation 13. So the methods are equivalent.

For the second case one differentiates Equation 15 to obtain

[TABLE]

Inserting this into Equation 15 gives

[TABLE]

which is similar to Equation 13 except that the measured $y_{i}$ have been replaced by the predicted $f_{i}$ everywhere except in the differences $\Delta_{i}=y_{i}-f_{i}$

Now consider which of Equations 14 and 15 should be used. In general the second choice, applying the error to the function, will be correct. If one measures a value $x=100.0$ with 10% precision one writes it as $100.0\pm 10.0$ for convenience, but this conceals the small difference between a true value of 110.0, for which the deviation is slightly under one sigma, and a true value of 90.0 for which it is slightly over. The situation is analogous to the choice between minimising $\sum{(n_{i}-f_{i})^{2}/f_{i}}$ or $\sum{(n_{i}-f_{i})^{2}/n_{i}}$ in a least squares fit to a histogram: the latter is quicker but leads to a bias which is noticeable for small $n_{i}$ .

5 Remarks on bias

That the matrix method, applying the multiplicative errors to the values, leads to a bias was pointed out by D’Agostini[4].

He gives a simple example in which the effect can be seen. Consider just two values with a shared multiplicative error $\xi$ and one just seeks a combined best value. A covariance matrix is constructed according to $V_{ij}=\delta_{ij}\sigma^{2}+\xi^{2}x_{1}x_{2}$ and inverted to form the $\chi^{2}$ to be minimised.

[TABLE]

The $2\times 2$ matrix can be inverted explicitly. Finding the maximum by differentiating and setting to zero gives

[TABLE]

The third term in the denominator clearly leads to bias in the estimator as it is positive definite and makes $\hat{x}$ smaller than the unbiassed weighted sum. This is illustrated in Figure 2 which shows results from 100,000 simulations of averaging two numbers, both with mean value 1 and individual errors 0.1, with various values for the multiplicative error $\xi$ . As the multiplicative error increases the spread in the results increases, as one would expect, but also a definite bias emerges. the mean is less than 1 and the distribution of results has a negative skew.

Figure 2: Bias for the average of two numbers with different values of the shared systematic $\xi$ . The results on the left are using Equation 21 and on the right using Equation 22.

This bias arises from the scaling, rather than the correlation. Indeed, if it were not for the correlation then the term that produces the bias would be larger: $(x_{1}^{2}+x_{2}^{2})\xi^{2}$ . It is a familiar effect when downward and upward fluctuations are (wrongly) assigned errors that are smaller and larger, respectively. The correlation changes this to $(x_{1}-x_{2})^{2}\xi^{2}$ : if the measurements happen to be the same the effect vanishes.

The alternative is to apply the factor to the result and minimise

[TABLE]

which gives an unbiassed result, as also shown in Figure 2. This takes noticeably longer to solve as $\hat{x}$ has to be found by iteration rather than a neat closed form solution, but this is usually a price worth paying.

6 Conclusions

The matrix inversion method and the additional-parameter method have been shown to be equivalent, for both additive and multiplicative systematic errors. Although this has been asserted and proved in some cases in the past, we have shown that it holds in general. This result deserves to be widely known, as means that those performing such fits can adopt whichever method they prefer without deliberating over which is more ‘correct’ in their case.

With either method, introducing multiplicative errors will lead to a bias in the result if they are applied to the measurements, but this can readily be avoided by applying the factors to the fitted values. This will make a linear problem become non-linear, but this is unlikely to present a real drawback.

References

[1] V. Branchina et al., Hadronic Width of the $Z$ from a Global Fit to $e^{+}e^{-}\to$ hadrons data, Phys. Rev. Lett. 65 3237 (1990)

[2] G.L. Fogli et al., Hints of $\theta_{13}>0$ from Global Neutrino Data Analysis, Phys. Rev. Lett 101 141801 (2008)

[3] R. B. Appleby et al., The Practical Pomeron for High Energy Proton Collimation, Eur Phys. J. C76 520 (2016)

[4] G. D’Agostini, On the use of the covariance matrix to fit correlated data, Nucl. Instr. & Meth. A 346 306 (1994)

[5] L Demortier Equivalence of the best-fit and covariance-matrix methods for comparing binned data with a model in the presence of correlated systematic uncertainties, note CDF/MEMO/STATISTICS/PUBLIC/8661 (1999)

[6] X-H Mo, Unbiased $\chi^{2}$ Estimator for linear function fit involving correlated data, HEP & NP 31 745 (2007)

[7] V Blobel, Data Combination in Particle Physics, Terascale School on data combination and Limit Setting, DESY, unpublished (2011)

[8] G L Fogli et al.,Getting the most from the statistical analysis of solar neutrino oscillations Phys. Rev. D66 05310 (2002)