Combining parameter values or $p$-values

We assume that there are N determinations $x_{i}\pm\sigma_{i}$ of a physical quantity $x$. A way of estimating $x_{comb}$ is to minimise the weighted sum of squares

$$S(x)=\Sigma[(x-x_{i})/\sigma_{i}]^{2},$$

(1)

with the summation running over the number of observations¹¹ 1 Many people refer to this as $\chi^{2}$. We prefer to use a different symbol ($S$), as this makes more understandable the question of whether or not the distribution of $S$ is the mathematical $\chi^{2}$.. This yields

$$x_{comb}=\Sigma w_{i}x_{i}/\Sigma{w_{i}},\ \ \ w_{i}=1/\sigma_{i}^{2}$$

(2)

i.e. the best value of $x$ is given by the weighted average of the $x_{i}$, where the weights are equal to $1/\sigma_{i}^{2}$. Thus the smaller the uncertainty on a measurement, the larger the weight. In an informal sense, the weight of an experiment can be thought of as its information content.

The uncertainty $\sigma_{comb}$ on $x_{comb}$ is given by

$$1/\sigma_{comb}^{2}=\Sigma(1/\sigma_{i}^{2})$$

(3)

This ensures that $\sigma_{comb}$ is at least as small as the smallest of the individual uncertainties; this is the motivation for combination. It is also guaranteed to be not larger than the uncertainty on the unweighted average. In terms of the weights, equation 3 is $w_{comb}=\Sigma{w_{i}}$, i.e. the information content of the combined value is the sum of those for each of the individual measurements.

The uncertainties in this note depend only on the uncertainties of the individual measurements, and their possible correlations, but not on the degree of consistency of the separate measurements. Thus the uncertainty on the combination of uncorrelated measurements $0\pm 3$ and $z\pm 3$ will be $\sim 2$, independent of whether $z$ is $1$ or $7$.

An example demonstrates that care is needed in applying the formulae. Consider high energy cosmic rays being recorded by a large counter system for two consecutive one-week periods, with the number of counts being $100\pm 10$ and $1\pm 1$ ²² 2 It is a crime (punishable by a forcible transfer to doing a doctorate on Astrology) to combine such discrepant measurements. It seems likely that someone turned off the detector between the two runs; or there was a large background in the first measurement which was eliminated for the second; etc. The only reason for using such discrepant numbers is to produce a dramatically stupid result. The effect would have been present with measurements like $100\pm 10$ and $81\pm 9$.. Unthinking application of the formulae for the combined result gives the ridiculous $2\pm 1$. What has gone wrong?

The answer is that we are supposed to use the true accuracies of the individual measurements to assign the weights. Here we have used the estimated accuracies. Because the estimated uncertainty depends on the estimated rate³³ 3 The problem arises here because the standard deviation in a Poisson process is equal to the square root of the rate. Even worse, in determining the lifetime $\tau$ of a particle from a set of $N$ measured decay times, the uncertainty on $\tau$ is $\tau/\sqrt{N}$ i.e. the estimated uncertainty is proportional to the estimated lifetime., a downward fluctuation in the measurement results in an underestimated uncertainty, an overestimated weight, and a downward bias in the combination. In our example, the combination should assume that the true rate was the same in the two measurements which used the same detector and which lasted the same time as each other, and hence their true accuracies are (unknown but) equal. So the two measurements should each be given the same weight, which yields the more sensible combined result of $50.5\pm 5$ counts per week.

A general way of mitigating this problem within this approach is discussed in ref [1]. It incorporates the way the uncertainty for each result is expected to vary with the estimated parameter value. It is equivalent to an iterative approach, in which at each stage, the input uncertainties are recalculated assuming that they can be obtained using the parameter value as determined in the previous iteration.

A method of combining correlated results is the ‘Best Linear Unbiassed Estimate’ (BLUE) [2]. We look for the best linear unbiassed combination

$$x_{BLUE}=\Sigma w_{i}x_{i},$$

(4)

where the weights are chosen to give the smallest uncertainty $\sigma_{BLUE}$ on $x_{BLUE}$. Also for the combination to be unbiassed, the weights must add up to unity. They are thus determined by minimising $\Sigma\Sigma w_{i}w_{j}E^{-1}_{ij}$, subject to the constraint $\Sigma w_{i}=1$; here $E$ is the covariance matrix for the correlated measurements. This gives

$$w_{i}=\Sigma e_{ij}/\Sigma\Sigma e_{ij},$$

(5)

where $e_{ij}$ is an element of $E^{-1}$, and the summation in the numerator is over the index j, while the double summation in the denominator is over $i$ and $j$.

The $BLUE$ procedure just described is equivalent to the $\chi^{2}$ approach for checking whether a correlated set of measurements are consistent with a common value. The advantage of $BLUE$ is that it provides the weights for each measurement in the combination. It thus enables us to calculate the contribution of various sources of uncertainty in the individual measurements to the uncertainty on the combined result.

When the correlation is so strong that the correlation coefficient⁴⁴ 4 Here $\sigma_{1}\leq\sigma_{2}$, and $\rho$ is the covariance divided by the product $\sigma_{1}\sigma_{2}$. $\rho>\sigma_{1}/\sigma_{2}$, the best estimate of $x$ falls outside the range of $x_{1}$ and $x_{2}$. This is in fact reasonable. If the correlation is strongly positive, it is likely that $x_{1}$ and $x_{2}$ lie on the same side of the true value $x_{true}$, with $x_{1}$ (the measurement with the smaller uncertainty) lying closer to $x_{true}$ than $x_{2}$ does. Thus it is entirely sensible that the best estimate should involve extrapolating from $x_{2}$ to beyond $x_{1}$.

However, the resulting $x_{BLUE}$ is sensitive to the values of the uncertainties and the correlations, so combining highly correlated values may not be sensible. This situation can arise, for example, when there is more than one group within a collaboration, analysing more or less the same data but with slightly different analyses, and with the same physics aim, e.g. measuring the top quark mass. This is likely to produce a set of answers that will be highly correlated. Rather than trying to combine the different results, it is better to decide which procedure should be used as the published result of the Collaboration (with the others being simply confirmatory). This choice should be based not on the results of the different methods, but rather on the expected sensitivity of each method; etc.

Another feature of large correlations is that the uncertainty on the combined value tends to zero as the correlation coefficient tends to +1 or -1. (Remember that even with complete correlation, the uncertainties do not have to be equal.)

When the individual measurements are uncorrelated, ${\bf BLUE}$ simplifies to the method described in Section 2.1.

Consider a remote island whose inhabitants are very conservative, and no-one leaves or arrives except for some anthropologists who wish to determine the number of married people there. Because the islanders are very traditional, it is necessary to send two teams of anthropologists, one consisting of males to interview the men, and the other of females for the women. There are too many islanders to interview them all, so each team interviews a sample and then extrapolates. The first team estimates the number of married men as $10,000\pm 300$. The second, who unfortunately have less funding and so can interview only a smaller sample, have a larger statistical uncertainty; they estimate $9,000\pm 900$ married women. Then how many married people are there on the island?

The simple approach is to add the numbers of married men and women, to give $19,000\pm 950$ married people. But if we use some theoretical input, maybe we can improve the accuracy of our estimate. So if we assume that the islanders are monogamous, the true numbers of married men and women should be equal. The weighted average is $9,900\pm 285$ married couples and hence $19,800\pm 570$ married people.

The contrast in these results is not so much the difference in the estimates, but that incorporating the assumption of monogamy and hence using the weighted average gives a smaller uncertainty on the answer. Of course, if our assumption is incorrect, this answer will be biassed.

A Particle Physics example incorporating the same idea of theoretical input reducing the uncertainty of a measurement is ‘Kinematic Fitting’. There the uncertainties on the measured momenta and energies of the objects produced in high energy interactions are reduced by assuming that energy and momentum conservation applies between the initial state collision particles and the final state objects measured in the detector.

There are situations where analyses determine two or more parameters. For example:

• •

  When fitting a peak plus a smooth background to a mass spectrum, the parameters include the location and strength of the signal, and perhaps its width.

• •

  In a search for neutrino oscillations where only two flavours are relevant, the parameters are the amplitude of the oscillations $\sin^{2}2\theta$; and $\Delta m^{2}$, the difference of the mass-squared of the two neutrinos, which determines the oscillation frequency.

• •

  One or more physics parameters $\phi$ and nuisance parameter(s) $\nu$ for systematic effects.

• •

  Straight line fitting. The parameters are the gradient and the intercept of the line.

In many cases the uncertainties on the parameters will be correlated.

When several independent⁵⁵ 5 If there are correlations among the $N$ separate measurements of the two parameters, the covariance matrix $E$ is expanded to be of size $2N\times 2N$, and equation (6) is readily modified to include all correlations – see Section 3.2. measurements of the correlated parameters exist, we may want to combine the results. For a pair of parameters, the weighted sum of squares for consistency with a particular $(a,b)$ is

$$S(a,b)=\Sigma[(a_{i}-a)^{2}e_{i,11}+(b_{i}-b)^{2}e_{i,22}+2(a_{i}-a)(b_{i}-b)e_{i,12}]$$

(6)

where the summation is over the $N$ independent measurements $(a_{i},b_{i})$, and $e_{i,kl}$ is the $(k,l)$ element of the inverse covariance matrix $E^{-1}_{i}$ of the $i^{th}$ measurement. Then $S(a,b)$ is simply minimised with respect to the parameters $a$ and $b$. The uncertainties and correlation for the combined values are given by the covariance matrix $M$; the elements of its inverse are

$$\begin{split}M^{-1}_{11}&=0.5*\frac{\partial^{2}S}{\partial a^{2}}\\ M^{-1}_{22}&=0.5*\frac{\partial^{2}S}{\partial b^{2}}\\ M^{-1}_{12}&=0.5*\frac{\partial^{2}S}{\partial a\partial b}\end{split}$$

(7)

The extension to the case where each analysis measures more than two parameters is straightforward.

An extension of the above example is where we have $N$ observables, each of which is measured in $p$ different experiments, and there are possible correlations in all $n=N\times p$ variables. Valassi has extended BLUE, using the criterion of minimising the uncertainty on each of the $N$ combined values [3]. For Gaussian uncertainties, this is shown to be equivalent to minimising the weighted sum of squares. An example where this might be used would be a measurement of the differential cross section for some process, using several different decay channels; there could be correlations across the bins of the cross-section for a given channel, and also among the different channels. The output is the differential cross-section, with each bin being the combination of the different channels, taking all correlations into account.

The Valassi procedure is as follows. Let us assume that we are interested in $N$ observables, $X_{\alpha}=\{X_{1},...,X_{N}\}$, and that we have $n$ experimental results $y_{i}=\{y_{1},...,y_{n}\}$, such that each of the measurements $y_{i}$ corresponds to one of the observables $x_{\alpha}$ (and all observables are measured at least once: $n\geq N$). The $(n\times N)$ matrix $\mathscr{U}$ is defined by

$$\mathscr{U}_{i\alpha}=\left\{\begin{matrix}[l]1&\textrm{if }y_{i}\textrm{ is a measurement of }X_{\alpha},\\ 0&\textrm{if }y_{i}\textrm{ is not a measurement of }X_{\alpha}.\end{matrix}\right.$$

(8)

Each of the $n$ rows of $\mathscr{U}$ has one and only one element equal to 1. For instance, if we would combine a 3-bin differential cross section measurement between two channels, e.g. muon and electron, then $n=6$, $N=3$, and the $\mathscr{U}$ matrix would be:

$$\mathscr{U}_{i\alpha}=\begin{pmatrix}1&0&0\\ 0&1&0\\ 0&0&1\\ 1&0&0\\ 0&1&0\\ 0&0&1\end{pmatrix}$$

(9)

Let us also define the $(n\times n)$ covariance matrix of the measurements,

$$\mathscr{M}_{ij}=\textrm{cov}(y_{i},y_{j}).$$

(10)

Then the best linear estimate of each observable $X_{\alpha}$ is

$$\hat{x}_{\alpha}=\sum_{i=1}^{n}\lambda_{\alpha i}y_{i},$$

(11)

where the weights $\lambda_{\alpha i}$ are

$$\lambda_{\alpha i}=\sum_{\beta=1}^{N}\left(\mathscr{U}^{t}\mathscr{M}^{-1}\mathscr{U}\right)^{-1}_{\alpha\beta}\left(\mathscr{U}^{t}\mathscr{M}^{-1}\right)_{\beta i},$$

(12)

and $\mathscr{U}^{t}$ is the transpose matrix of $\mathscr{U}$. The covariance matrix for the estimates is

$$\textrm{cov}(\hat{x}_{\alpha},\hat{x}_{\beta})=\left(\mathscr{U}^{t}\mathscr{M}^{-1}\mathscr{U}\right)^{-1}_{\alpha\beta}.$$

(13)

Here we discuss a very simple example of combining straight line fits, where the answer can be appreciated intuitively. It consists of a simplified tracking situation, in which a particle passes through 6 detector planes, of which 3 are closely spaced and separated by some distance from another 3 closely spaced planes (See Fig. 1). The data consist of independent measurements $y_{i}\pm\sigma_{i}$ at well-defined $x_{i}$ values. There is no magnetic field and we consider only the $x$ and $y$ coordinates so the track is parametrised by the straight line $y=a+bx$. A straight line $L_{1}$ is fitted to the hits in the 3 left-most planes, and $L_{2}$ to the 3 right-most planes. Finally the results $(a_{1},b_{1})$ and $(a_{2},b_{2})$ are combined to give $(a_{comb},b_{comb})$.

A straight line fit to a set of closely separated points will determine the line’s gradient $b$ with a large uncertainty. Furthermore, if these points are centred away from $x=0$, there will be a strong correlation between $a$, the intercept at $x=0$, and $b$. The large uncertainty on the gradient $b$ then results in a large uncertainty on the intercept $a$. The covariance of $a$ and $b$ is obtained from equations 7, and is proportional to $-\langle x\rangle$, where $\langle x\rangle$ is the weighted average of the $x$-positions of the fitted data points (i.e. $\Sigma(x_{i}/\sigma_{i}^{2})/\Sigma(1/\sigma_{i}^{2})$).

Thus we expect that the lines $L_{1}$ and $L_{2}$ will have large uncertainties on $a$ and $b$, and strong correlations but of opposite signs. However, because of the larger range of $x$-values, $L_{comb}$ will have very much smaller uncertainties. The covariance ellipses for $L_{1}$ and $L_{2}$, as well as for $L_{comb}$ are shown in Fig. 1, and bear out these expectations. In terms of the covariance ellipses shown there, it is because they have different orientations that the combination results in vastly reduced uncertainties.

When the two sets of sub-detector planes are centred on the same side of the origin, as is usually the case in tracking, the best values of both the gradient and of the intercept of the combined line can be outside the ranges of the corresponding quantities for lines $L_{1}$ and $L_{2}$ (see Fig. 2).

For the straight line fits, the results for $(a,b)$ and for their covariance matrix $M$ are the same whether we combine the $L_{1}$ and $L_{2}$ results, or whether we do a single fit to all 6 data points.

A physical example of the big reduction in uncertainty when combining results of pairs of parameters with large internal correlations is the determination of the fraction of dark energy in the Universe $\Omega_{\Lambda}$. There are several different methods that provide information on this and the fraction of dark matter $\Omega_{M}$, but each on its own has a large uncertainty. Because of their different correlations, however, their combination provides a precise measurement of $\Omega_{\Lambda}$ (see Fig. 3).

For situations where we have several parameters, it is common and sometimes natural to choose one as the parameter of interest ($\phi$) and to profile or to marginalise over all the others ($\nu$); this is especially common when we have one physics parameter, and several nuisance parameters related to systematic effects. The profile likelihood is

$$\mathcal{L}_{prof}(\phi)={\mathcal{L}}(\phi,\nu_{best}(\phi))$$

(14)

where $\nu_{best}(\phi)$ is the value for $\nu$ that maximise the likelihood at that particular $\phi$. The profile likelihood is thus a function just of $\phi$. Fig. 2(c) shows the profile likelihood for the intercept $a$ (i.e. profiled over the gradient $b$) for all three lines. The widths of the curves for $\ln\mathcal{L}_{prof}$ correctly give the uncertainties on $a$. However it is important to note that, in contrast to the situation with the full likelihoods $\mathcal{L}(a,b)$, combining the profile likelihoods for $L_{1}$ and $L_{2}$ would not give the profile likelihood for $\mathcal{L}_{comb}$, even though the best values of $b$ for the two lines happen to be the same. However, except when $a=a_{best}$, the values of $b_{best}(a)$ for $L_{1}$ and for $L_{2}$ at the same $a$ are different - see Fig. 2(d)); that is a reason that the combination of profile likelihoods is not a sensible procedure.

If the only information available is the set of values $\mu_{i}$ from the separate measurements and their covariance matrix , then the only possibilities for combining are the methods described above. With a little more information, however, it might be possible to reconstruct approximately the likelihood functions and then to combine them rather than just the results.

For example, when the uncertainties are asymmetric, Barlow has suggested various ways of modifying the Gaussian shape by prescriptions for how the width varies [6].

A note by Cousins [7] points out that in some simplified circumstances, the shape of the likelihood is defined. Thus if the result is obtained by multiplying several factors, the equivalent of the central limit theorem causes the distribution of the product to be approximately log-normal. Some systematic uncertainties might also result in such a distribution. For example, theorists might say that their predicted cross-section for some process was accurate to within a factor of 2.

Alternatively for a search for a signal involving Poisson counting in the signal and in background regions (the ‘on-off’ problem), the distribution is a gamma function. This also applies to lifetime determinations using individually observed decay times with an expected exponential distribution. In all cases the parameters of the expected distributions are determined from the numerical values of $\mu$ and $\sigma$.

Of course the ideal situations for these distributions to be relevant are rarely realised in practice. For example, for the gamma distribution to apply to the lifetime measurement, we require the expected decay distribution to follow a perfect exponential. This means that we have a constant efficiency for observing decays over the full range of decay times from zero to infinity, and can ignore backgrounds, time resolution, etc.

Because correlations can lead to extrapolation and to small uncertainties, it is sometimes suggested that it would be a good idea to set the correlation coefficient $\rho$ to zero. This is thought to be conservative, but it throws away information and is against the spirit of BLUE - ‘B’ stands for ‘Best’, which means ‘smallest uncertainty’, rather than ‘most conservative’. Furthermore $\rho=0$ is not the most conservative choice. For example, if we have analysed a large data set $D$ and also a subset $S$ of this data, and we combine them using the covariance matrix with elements

$$E_{11}=\sigma_{D}^{2},\ \ \ E_{22}=\sigma_{S}^{2},\ \ \ E_{12}=E_{21}=\sigma_{S}^{2},$$

(15)

the weight ascribed to $S$ turns out to be zero, i.e. the subsample $S$ is ignored and the ‘combined’ result is simply that from $D$ (as is sensible). But if we set the covariance to zero, our combined result will have an incorrect ‘improvement’ with reduced uncertainty $\sigma_{C}$ given by

$$1/\sigma_{C}^{2}=1/\sigma_{D}^{2}+1/\sigma_{S}^{2}$$

(16)

While it might be sensible to choose the most conservative value for $\rho$ when the value of $\rho$ is unknown, otherwise its actual value seems a better choice.

Another problem is that in ignoring $\rho$, we will obtain an incorrect contribution to the weighted sum of squares $S$. Thus if the $pdf$ in $(x,y)$ is a 2-dimensional Gaussian centred on the origin with $\sigma_{x}=\sigma_{y}=1$ and $\rho=+0.9$, the correct contribution to $S$ from a measurement at (+1.0, +1.0) is 1, while for (+1.0, -1.0) it is 20; setting $\rho=0$ would incorrectly result in a contribution of 2 for both of them.