40. Statistics

1

40. Statistics
Revised August 2023 by G. Cowan (RHUL).
This chapter gives an overview of statistical methods used in high-energy physics. In statistics,
we are interested in using a given sample of data to make inferences about a probabilistic model,e.g., to assess the model’s validity or to determine the values of its parameters. There are two main
approaches to statistical inference, which we may call frequentist and Bayesian.
In frequentist statistics, probability is interpreted as the limiting frequency of the outcome ofa repeatable experiment. The most important tools in this framework are parameter estimation,
covered in Section 40.2, statistical tests, discussed in Section 40.3, and confidence intervals, which
are constructed so as to cover the true value of a parameter with a specified probability, as describedin Section 40.4.2. Note that in frequentist statistics one does not define a probability for a hypothesis
or for the value of a parameter.
In Bayesian statistics, the subjective interpretation of probability is used to quantify one’s
degree of belief in a hypothesis. This allows one to define a probability density function (p.d.f.) fora parameter, which reflects one’s knowledge about where its true value lies.
Bayesian methods provide a natural means to include additional information, which in general
may be subjective; in fact they require prior probabilities for the hypotheses (or parameters) inquestion, i.e., the degree of belief about the parameters’ values, before carrying out the measurement. Using Bayes’ theorem (Eq. (39.4)), the prior degree of belief is updated by the data from theexperiment. Bayesian methods for interval estimation are discussed in Sections 40.4.1 and 40.4.2.4.
For many inference problems, the frequentist and Bayesian approaches give similar numericalvalues, even though they answer different questions and are based on fundamentally different interpretations of probability. In some important cases, however, the two approaches may yield very
different results. For a discussion of Bayesian vs. non-Bayesian methods, see references written bya statistician [1], by a physicist [2], or the detailed comparison in Ref. [3].

40.1 Fundamental concepts
Consider an experiment whose outcome is characterized by one or more data values, which
we can write as a vector x. A hypothesis H is a statement about the probability for the data,
often written P (x|H). (We will usually use a capital letter for a probability and lower case fora probability density. Often the term p.d.f. is used loosely to refer to either a probability or aprobability density.) This could, for example, define completely the p.d.f. for the data (a simple
hypothesis), or it could specify only the functional form of the p.d.f., with the values of one or moreparameters not determined (a composite hypothesis).
If the probability P (x|H) for data x is regarded as a function of the hypothesis H, then it iscalled the likelihood of H, usually written L(H). Often the hypothesis is characterized by one ormore parameters θ, in which case L(θ) = P (x|θ) is called the likelihood function.
In some cases one can obtain at least approximate frequentist results using the likelihood
evaluated only with the data obtained, for example, when determining confidence regions with
Eq. (40.79). In general, however, the frequentist approach requires a full specification of the probability model P (x|H) both as a function of the data x and hypothesis H.
In the Bayesian approach, inference is based on the posterior probability for H given the datax, which represents one’s degree of belief that H is true given the data. This is obtained from
Bayes’ theorem (39.4), which can be written
P (H|x) = R

P (x|H)π(H).
P (x|H 0 )π(H 0 ) dH 0

S. Navas et al. (Particle Data Group), Phys. Rev. D 110, 030001 (2024)
31st May, 2024 10:20am(40.1)

40. Statistics

2

Here P (x|H) is the likelihood for H, which depends only on the data actually obtained. The
quantity π(H) is the prior probability for H, which represents one’s degree of belief for H before
carrying out the measurement. The integral in the denominator (or sum, for discrete hypotheses)
serves as a normalization factor. If H is characterized by a continuous parameter θ then theposterior probability is a p.d.f. p(θ|x). Note that the likelihood function itself is not a p.d.f. for θ.

40.2 Parameter estimation
Here we review point estimation of parameters, first with an overview of the frequentist approach and its two most important methods, maximum likelihood and least squares, treated in
Sections 40.2.2 and 40.2.3. The Bayesian approach is outlined in Sec. 40.2.6.
An estimator θb (written with a hat) is a function of the data used to estimate the value of theparameter θ. Sometimes the word ‘estimate’ is used to denote the value of the estimator whenevaluated with given data. There is no fundamental rule dictating how an estimator must be
constructed. One tries, therefore, to choose that estimator which has the best properties. The
most important of these are (a) consistency, (b) bias, (c) efficiency, and (d) robustness.
(a) An estimator is said to be consistent if the estimate θb converges in probability (see Ref. [3]) to
the true value θ as the amount of data increases. This property is so important that it is possessedby all commonly used estimators.
(b) The bias, b = E[ θb ] − θ, is the difference between the expectation value of the estimator andthe true value of the parameter. The expectation value is taken over a hypothetical set of similar
experiments in which θb is constructed in the same way. When b = 0, the estimator is said to beunbiased. The bias depends on the chosen metric, i.e., if θb is an unbiased estimator of θ, then θb 2
is not in general an unbiased estimator for θ2 .
(c) Efficiency is the ratio of the minimum possible variance for any estimator of θ to the varianceb For the case of a single parameter, under rather general conditions theV [ θb ] of the estimator θ.
minimum variance is given by the Rao-Cramér-Fréchet bound,
2
σmin
=

where

"

I(θ) = E



∂b
1+
∂θ

∂ ln L
∂θ

2 #

2

/I(θ) ,"

∂ 2 ln L= −E∂θ2(40.2)#(40.3)

is the Fisher information, L is the likelihood, and the operator E[] in (40.3) is the expectation
value with respect to the data. For the final equality to hold, the range of allowed data values must
not depend on θ.
The mean-squared error,b + b2 ,MSE = E[(θb − θ)2 ] = V [θ](40.4)
is a measure of an estimator’s quality which combines bias and variance.
(d) Robustness is the property of being insensitive to departures from assumptions about the p.d.f.,e.g., owing to uncertainties in the distribution’s tails.
It is not in general possible to optimize simultaneously for all the measures of estimator quality
described above. For some common estimators, the properties above are known exactly. Moregenerally, it is possible to evaluate them by Monte Carlo simulation. Note that they will in generaldepend on the unknown θ.
40.2.1 Estimators for mean, variance, and medianSuppose we have a set of n independent measurements, x1 , . . . , xn , each assumed to follow ap.d.f. with unknown mean µ and unknown variance σ 2 (the measurements do not necessarily have
31st May, 2024

40. Statistics

3
to follow a Gaussian distribution). Then

b=
µ

n
1X
xi
n i=1

(40.5)

c2 =
σ

n
1 X
b )2
(xi − µn − 1 i=1

(40.6)

c2 is
b is σ 2 /n and the variance of σ
are unbiased estimators of µ and σ 2 . The variance of µh ic2 = 1 m − n − 3 σ 4
V σ
4
n
n−1




,

(40.7)

where m4 is the 4th central moment of x (see Eq. (39.8)). For Gaussian distributed
xi , this√b is σ/ 2n. For any
becomes 2σ 4 /(n − 1) for any n ≥ 2, and for large n the standard deviation of σc2 are uncorrelated.
b is an efficient estimator for µ, and the estimators µb and σn and Gaussian xi , µ
Otherwise the arithmetic mean (40.5) is not necessarily the most efficient estimator; this is discussed
further in Sec. 8.7 of Ref. [4].
b, as can be seen from Eq. (40.5); however, if
If σ 2 is known, it does not improve the estimate µb in Eq. (40.6) and replace n − 1 by n to obtain an estimatorµ is known, one can substitute it for µof σ 2 still with zero bias but smaller variance. If the xi have different, known variances σi2 , then
the weighted average
n
1X
b=
µ
wi x i ,(40.8)
w i=1
where wi = 1/σi2 and w = i wi , is an unbiased estimator for µ with a smaller variance than an√
b is 1/ w.
unweighted average. The standard deviation of µbmed such that half the xi are below
As an estimator for the median xmed , one can use the value x
and half above (the sample median). If there are an even number of observations and the sample
median lies between two observed values, the estimator is set by convention to their arithmetic
average. If the p.d.f. of x has the form f (x − µ) and µ is both mean and median, then for largen the variance of the sample median approaches 1/[4nf 2 (0)], provided f (0) > 0 [5]. Althoughestimating the median can often be more difficult computationally than the mean, the resulting
estimator is generally more robust, as it is insensitive to the exact shape of the tails of a distribution.
P

40.2.2 The method of maximum likelihood
Suppose we have a set of measured quantities x and the likelihood L(θ) = P (x|θ) for a set ofparameters θ = (θ1 , . . . , θM ). The maximum likelihood estimators (MLEs) for θ are defined as thevalues that give the maximum of L. Because of the properties of the logarithm, it is usually easier
to work with ln L, and since both are maximized for the same parameter values θ, the MLEs canbe found by solving the likelihood equations,∂ ln L=0,∂θii = 1, . . . , M .

(40.9)

Often the solution must be found numerically. Maximum likelihood estimators are important because they are asymptotically (i.e., for large data samples) unbiased, efficient and have a Gaussian
sampling distribution under quite general conditions, and the method has a wide range of applicability.
In general the likelihood function is obtained from the probability of the data under assumptionof the parameters. An important special case is when the data consist of i.i.d. (independent
31st May, 2024

440. Statistics

and identically distributed) values. Here one has a set of n statistically independent quantities
x = (x1 , . . . , xn ), where each component follows the same p.d.f. f (x; θ). In this case the joint p.d.f.
of the data sample factorizes and the likelihood function isL(θ) =nYf (xi ; θ) .

(40.10)i=1Here the number n of observations (usually individual “events” in particle physics) is regarded as
fixed. If, however, the probability to observe n events itself depends on the parameters θ, then thisdependence should be included in the likelihood. For example, if n follows a Poisson distribution
with mean µ and the independent x values all follow f (x; θ), then the likelihood becomesL(θ) =nµn −µ Yef (xi ; θ) .
n!
i=1

(40.11)Equation (40.11) is often called the extended likelihood (see, e.g., Refs. [6–8]). If µ is given as afunction of θ, then including the probability for n given θ in the likelihood provides additional information about the parameters. This therefore leads to a reduction in their statistical uncertainties
and in general changes their estimated values.
In evaluating the likelihood function, it is important that any normalization factors in the p.d.f.
that involve θ be included. However, we will only be interested in the maximum of L and in ratiosof L at different values of the parameters; hence any multiplicative factors that do not involve theparameters that we want to estimate may be dropped, including factors that depend on the data
but not on θ.
b That is,Under a one-to-one change of parameters from θ to η, the MLEs θb transform to η(θ).
the maximum-likelihood solution is invariant under change of parameter. However, other properties
of MLEs, in particular the bias, are not invariant under change of parameter.
The inverse V −1 of the covariance matrix Vij = cov[θbi , θbj ] for a set of MLEs can be estimatedby using∂ 2 ln L(Vb −1 )ij = −.
(40.12)∂θi ∂θj θb
Equation (40.12) holds for a sufficiently large data sample (or for a small sample only in special
cases, e.g., where the means of Gaussian distributed data are linear functions of the parameters);
otherwise it can result in a misestimation of the variances. Under the conditions where the equation
is valid, L has a Gaussian form and ln L is (hyper)parabolic. In this case, s times the standard
deviations σi of the estimators for the parameters can be obtained from the hypersurface definedby the θ such thatln L(θ) = ln Lmax − s2 /2 ,(40.13)where ln Lmax is the value of ln L at the solution point (compare with Eq. (40.79)). The minimum
and maximum values of θi on the hypersurface then give an approximate s-standard deviation
confidence interval for θi (see Section 40.4.2.2).
40.2.2.1 Maximum likelihood with binned data
If the total number of data values xi , (i = 1, . . . , ntot ), is small, the unbinned maximumlikelihood method, i.e., use of Equation (40.10) (or (40.11) for extended maximum likelihood), ispreferred since binning can only result in a loss of information, and hence larger statistical errors
for the parameter estimates. If the sample is large, it can be convenient to bin the values in ahistogram with N bins, so that one obtains a vector of data n = (n1 , . . . , nN ) with expectation
31st May, 2024

540. Statistics

values µ = E[n] and probabilities f (n; µ). Suppose the mean values µ can be determined as afunction of a set of parameters θ. Then one may maximize the likelihood function based on thecontents of the bins.
P
As mentioned in Sec. 40.2.2, the total number of events ntot = i ni can be regarded either as
fixed or as a random variable. If it is fixed, the histogram follows a multinomial distribution,fM (n; θ) =ntot !
pn1 · · · pnNN ,n1 ! · · · nN ! 1(40.14)

where we assume the probabilities pi are given functions of the parameters θ. The distribution canbe written equivalently in terms of the expected number of events in each bin, µi = ntot pi . If theni are regarded as independent and Poisson distributed, then the data are instead described by a
product of Poisson probabilities,fP (n; θ) =

N
Y
µni
i

n!
i=1 i

e−µi ,(40.15)where the mean values µi are given functions of θ. The total number of events ntot thus follows aP
Poisson distribution with mean µtot = i µi .
When using maximum likelihood with binned data, one can find the estimators and at thesame time obtain a statistic usable for a test of goodness of fit (see Sec. 40.3.3.1). For independent
Poisson distributed data, maximizing the likelihood L(θ) = fP (n; θ) is equivalent to maximizing
the likelihood ratio λ(θ) = fP (n; θ)/fsat (n; µ̂), where in the denominator fsat (n; µ) is the corresponding model with an adjustable parameter for each bin, µ = (µ1 , . . . , µN ), with estimatorsµ̂ = (n1 , . . . , nN ) (called the “saturated model”). Equivalently one often minimizes the quantity−2 ln λ(θ). For independent Poisson distributed ni this is [9]
N 
X

ni
µi (θ) − ni + ni ln− 2 ln λ(θ) = 2
µ
i (θ)
i=1



,

(40.16)

where for bins with ni = 0, the last term in (40.16) is zero. The expression (40.16) without theterms µi − ni also gives −2 ln λ(θ) for multinomially distributed ni , i.e., when the total numberof entries is regarded as fixed. In the limit of zero bin width, minimizing (40.16) is equivalent tomaximizing the unbinned extended likelihood function (40.11); in the corresponding multinomialcase without the µi − ni terms one obtains Eq. (40.10).
b corresponds to better agreement between the data and the hyA smaller value of −2 ln λ(θ)
b can thus be translated into a p-value as a measure
pothesized form of µ(θ). The value of −2 ln λ(θ)
of goodness-of-fit, as described in Sec. 40.3.3.1. Assuming the model is correct, then according to
Wilks’ theorem [10], for sufficiently large µi and provided certain regularity conditions are met, theminimum of −2 ln λ as defined by Eq. (40.16) follows a χ2 distribution (see, e.g., Ref. [9]). If there
are N bins and M fitted parameters, then the number of degrees of freedom for the χ2 distribution
is N − M if the data are treated as Poisson-distributed, and N − M − 1 if the ni are multinomiallydistributed.
P
Suppose the ni are Poisson-distributed and the overall normalization µtot = i µi is taken as an
adjustable parameter, so that µi = µtot pi (θ), where the probability to be in the ith bin, pi (θ), does
not depend on µtot . Then by minimizing Eq. (40.16), one obtains that the area under the fittedPPfunction is equal to the sum of the histogram contents, i.e., i µ̂i = i ni . This is a property notpossessed by the estimators from the method of least squares (see, e.g., Sec. 40.2.3 and Ref. [8]).
31st May, 2024

6

40. Statistics

40.2.2.2 Frequentist treatment of nuisance parameters
Suppose we want to determine the values of parameters θ using a set of measurements x
described by a probability model P (x|θ). In general the model is not perfect, which is to say it
cannot provide an accurate description of the data even at the most optimal point of its parameter
space. As a result, the estimated parameters can have a systematic bias.
One can improve the model by including in it additional parameters. That is, P (x|θ) is replaced
by a more general model P (x|θ, ν), which depends on parameters of interest θ and nuisanceparameters ν. The additional parameters are not of intrinsic interest but must be included for themodel to be sufficiently accurate for some point in the enlarged parameter space.
Although including additional parameters may eliminate or at least reduce the effect of systematic uncertainties, their presence will result in increased statistical uncertainties for the parametersof interest. This occurs because the estimators for the nuisance parameters and those of interest
will in general be correlated, which results in an enlargement of the contour defined by Eq. (40.13).
To reduce the impact of the nuisance parameters one often tries to constrain their values bymeans of control or calibration measurements, say, having data y. For example, some components ofy could represent estimates of the nuisance parameters, often from separate experiments. Suppose
the measurements y are statistically independent from x and are described by a model P (y|ν).
The joint model for both x and y is in this case therefore the product of the probabilities for x
and y, and thus the likelihood function for the full set of parameters isL(θ, ν) = P (x|θ, ν)P (y|ν) .

(40.17)Note that in this case if one wants to simulate the experiment by means of Monte Carlo, boththe primary and control measurements, x and y, must be generated for each repetition underassumption of fixed values for the parameters θ and ν.
Using all of the parameters (θ, ν) in Eq. (40.13) to find the statistical errors in the parametersof interest θ is equivalent to using the profile likelihood, which depends only on θ. It is defined asbb
Lp (θ) = L(θ, ν(θ)),(40.18)where the double-hat notation indicates the profiled values of the parameters ν, defined as the values
that maximize L for the specified θ. The profile likelihood is discussed further in Section 40.3.2.1
in connection with hypothesis tests.
40.2.3 The method of least squares
The method of least squares (LS) coincides with the method of maximum likelihood in thefollowing special case. Consider a set of N independent measurements yi at known points xi . The
measurement yi is assumed to be Gaussian distributed with mean µ(xi ; θ) and known variance σi2 .
The goal is to construct estimators for the unknown parameters θ. The log-likelihood function
contains the sum of squares
χ2 (θ) = −2 ln L(θ) + constant =NX(yi − µ(xi ; θ))2
i=1

σi2

.

(40.19)The parameter values that maximize L are the same as those which minimize χ2 .
The minimum of the chi-square function in Equation (40.19) defines the least-squares estimators θb for the more general case where the yi are not Gaussian distributed as long as they are
independent. If they are not independent but rather have a covariance matrix Vij = cov[yi , yj ],then the LS estimators are determined by the minimum ofχ2 (θ) = (y − µ(θ))T V −1 (y − µ(θ)) ,31st May, 2024

(40.20)40. Statistics

7

where y = (y1 , . . . , yN ) is the (column) vector of measurements, µ(θ) is the corresponding vector
of predicted values, and the superscript T denotes the transpose. If the yi are not Gaussian
distributed, then the least-squares and maximum-likelihood estimators will not in general coincide.
Often one further restricts the problem to the case where µ(xi ; θ) is a linear function of theparameters, i.e.,µ(xi ; θ) =mXθj hj (xi ) .

(40.21)j=1Here the hj (x) are m linearly independent functions, e.g., 1, x, x2 , . . . , xm−1 or Legendre polynomials. We require m < N and at least m of the xi must be distinct.
Minimizing χ2 in this case with m parameters reduces to solving a system of m linear equations.
Defining Hij = hj (xi ) and minimizing χ2 by setting its derivatives with respect to the θi equal tozero gives the LS estimators,
θb = (H T V −1 H)−1 H T V −1 y ≡ Dy .

(40.22)The covariance matrix for the estimators Uij = cov[θbi , θbj ] is given by
U = DV DT = (H T V −1 H)−1 ,(40.23)or equivalently, its inverse U −1 can be found from(U −1 )ij =NX
1 ∂ 2 χ2=hi (xk )(V −1 )kl hj (xl ) .
2 ∂θi ∂θj θ=θb k,l=1(40.24)

The LS estimators can also be found from the expression
θb = U g ,(40.25)

yj hi (xk )(V −1 )jk .

(40.26)

where the vector g is defined by
gi =

N
X
j,k=1

For the case of uncorrelated yi , for example, one can use (40.25) with
(U

−1

)ij =
gi =

N
X
hi (xk )hj (xk )
k=1
N
X

σk2

,

yk hi (xk )
.
σk2
k=1

(40.27)(40.28)b one finds that the contour in parameter space defined by
Expanding χ2 (θ) about θ,b + 1 = χ2 + 1
χ2 (θ) = χ2 (θ)
min(40.29)has tangent planes located at plus-or-minus-one standard deviation σθb from the LS estimates θb(the relation is approximate if the fit function µ(x; θ) is nonlinear in the parameters).
In constructing the quantity χ2 (θ) one requires the variances or, in the case of correlated
measurements, the covariance matrix. Often these quantities are not known a priori and must be
31st May, 2024

8

40. Statistics

estimated from the data. In this case the least-squares and maximum-likelihood methods are nolonger exactly equivalent even for Gaussian distributed measurements. An important example iswhere the measured value yi represents the event count in a histogram bin. If, for example, yi
represents a Poisson variable, for which the variance is equal to the mean, then one can either
estimate the variance from the predicted value, µ(xi ; θ), or from the observed number itself, yi . In
the first option, the variances become functions of the parameters, and as a result the estimators
may need to be found numerically. The second option can be undefined if yi is zero, and for small yi ,the variance will be poorly estimated. In either case, one should constrain the normalization of thefitted curve to the correct value, i.e., one should determine the area under the fitted curve directlyfrom the number of entries in the histogram (see Ref. [8], Section 7.4). As noted in Sec. 40.2.2.1,
this issue is avoided when using the method of extended maximum likelihood with binned data by
minimizing Eq. (40.16). In that case if the expected number of events µtot does not depend on theother fitted parameters θ, then its extended MLE is equal to the observed total number of events.
As the minimum value of the χ2 represents the level of agreement between the measurements
and the fitted function, it can be used for assessing the goodness-of-fit; this is discussed further in
Section 40.3.3.1.
40.2.4 Parameter estimation with constraints
In some applications one is interested in using a set of measured quantities y = (y1 , . . . , yN ) toestimate parameters θ = (θ1 , . . . , θM ) subject to a number of constraints. For example, one may
have measured coordinates from two tracks, and the goal is to estimate their momentum vectors
subject to the constraint that the tracks have a common vertex. The parameters can also include
momenta of undetected particles such as neutrinos, as long as the constraints from conservation ofenergy and momentum and from known masses of particles involved in the reaction chain provide
enough information for these quantities to be inferred.
A set of K constraints can be given in the form of equations
ck (θ) = 0 ,k = 1, . . . , K .

(40.30)In some problems it may be possible to define a new set of parameters η = (η1 , . . . , ηL ) with
L = M − K such that every point in η-space automatically satisfies the constraints. If this ispossible then the problem reduces to one of estimating η with, e.g., maximum likelihood or leastsquares and then transforming the estimators back into θ-space.
In many cases it may be difficult or impossible to find an appropriate transformation η(θ).
Suppose that the parameters are determined by minimizing an objective function such as χ2 (θ) in
the method of least squares. Here one may enforce the constraints by finding the stationary pointsof the Lagrange function
L(θ, λ, y) = χ2 (θ, y) +KX

λk ck (θ)(40.31)k=1

with respect to both the parameters θ and a set of Lagrange multipliers λ = (λ1 , . . . , λK ). Combining the parameters and Lagrange multipliers into an (M +K)-component vector γ = (θ1 , . . . , θM , λ1 ,. . . , λK ), the solutions for γ, i.e., the estimators γ̂, are found (e.g., numerically) from the systemof equations∂L
Fi (γ, y) ≡= 0 , i = 1, . . . , M + K .
(40.32)∂γi
To obtain the covariance matrix of the estimated parameters one can find solutions γ̃ corresponding
to the expectation values of the data hyi and expand Fi (γ̂, y) to first order about these values.
31st May, 2024

40. Statistics

9

This gives (see, e.g., Sec. 11.6 of Ref. [8]) linearized approximations for the estimators, γ̂(y) ≈γ̃ + C(y − hyi), where the matrix C = −A−1 B, and A and B are given by
"

#

∂Fi
Aij =
∂γj γ̃,hyi

"

#

∂Fi
and Bij =.
∂yj γ̃,hyi(40.33)

In practice the values hyi and corresponding solutions γ̃ are estimated using the data from theactual measurement. Using this approximation for γ̂(y), one can find the covariance matrix Uij =cov[γ̂i , γ̂j ] of the estimators for the γi in terms of that of the data Vij = cov[yi , yj ] using error
propagation (cf. Eqs. (39.17) and (39.18)),
U = CV C T .

(40.34)The upper-left M × M block of the matrix U gives the covariance matrix for the estimated parameters cov[θ̂i , θ̂j ]. One can show for linear constraints that cov[θ̂i , θ̂j ] is also given by the upper-left
M × M block of 2A−1 . If the parameters are estimated using the method of least squares, then thenumber of degrees of freedom for the distribution of the minimized χ2 is increased by the numberof constraints, i.e., it becomes N − M + K. Further details can be found in, e.g., Ch. 8 of Ref. [4]and Ch. 7 of Ref. [11].
40.2.5 Unfolding
An important class of parameter estimation problem involves measurement of the differential
distribution of a kinematic variable in the form of a histogram with N bins. The data thus consistof the vector of measured data values n = (n1 , . . . nN ), with expectation values νi = E[ni ]. The ni
are usually independent and often modeled as Poisson distributed, from which it follows that the
maximum-likelihood estimators are ν̂i = ni for all i.
Because of the limited acceptance and resolution of the experiment, however, the measured
values of the kinematic variable in question differ in general from their true values, and as aconsequence the form of the data histogram is distorted relative to what would be obtained withperfect resolution. The desired parameters are not, therefore, the νi but rather one wants toestimate the expected number of entries in a given bin that would be found with a perfect detector.
We call this the “true histogram” and denote it with µ = (µ1 , . . . , µM ), where the number of bins
M is not required to equal that of the observed histogram. The µi are related to the expected
numbers of events in the observed histogram by
νi =

M
X

Rij µj + βi ,(40.35)j=1

where βi here represents the expected number of events in bin i due to background processes andthe N × M response matrix Rij gives the probability for an event to be observed in bin i of ν
given that the true value of the variable was in bin j of µ. For purposes of this discussion let us
suppose that the response matrix and expected background values are known. There are two main
approaches to this type of problem, which we can call unfolding and folding.
In unfolding, one treats the true histogram µ as the parameters of interest. The result isthus given by estimators µ̂ and the corresponding covariance matrix Uij = cov[µ̂i , µ̂j ]. Provided
the response matrix can be inverted, the maximum-likelihood solution is easily found as µ̂ =R−1 (n − β). If the response matrix allows for significant migration of events between bins, thenthe variances of these estimators can be very large, sometimes to the point where the µ̂ bear
essentially no resemblance to the true µ. In such cases the estimators can be found by maximizing
31st May, 2024

1040. Statistics

a linear combination of the log-likelihood and a regularization function that imposes some degreeof smoothness on the unfolded distribution. In achieving a reduction in variance one inevitably
introduces some bias into the estimators.
In the approach of folding, by contrast, to test a given model prediction for the true histogramµ it is first “folded” with the response matrix and corrected for expected background to give thecorresponding ν according to Eq. (40.35), and these are compared to the corresponding n using, e.g.,a likelihood ratio as shown in Sec. 40.2.2.1. In this way, one avoids any bias due to regularization.
To account for systematic uncertainties, the response matrix and background values may not be
expressed as constants but rather as functions of nuisance parameters θ. These may be constrainedby auxiliary measurements u with p.d.f. p(u|θ) so that the full likelihood becomes (as in Eq. (40.17)
but with different notation)
L(µ, θ) = P (n|ν(µ, θ))p(u|θ) .

(40.36)

Here the mean values ν = R(θ)µ + β(θ) now depend on both the true histogram parameters µ
as well as the nuisance parameters θ. Thus to record enough information for a future analysis
using folding one must provide all of the ingredients used above: the primary data n, the auxiliary
measurements u, and also the response matrix R(θ) and expected backgrounds β(θ) as functionsof the nuisance parameters. If unfolding is used, then future model tests or comparisons with otherexperiments can be carried out directly using the estimators µ̂ and their covariance matrix.
If several distributions are unfolded, then to combine these in a test of a given model one should
know how estimators for bins of different distributions are correlated, as can arise, e.g., throughcommon systematic effects. If only the unfolded distributions and their separate covariance matrices
are reported, however, then information on such correlations is not retained. In folding, one can
include information on correlated systematic effects if the ingredients (R and β) are known in termsof nuisance parameters that are common to different distributions.
In unfolding, the estimators for µ can be constructed either by maximizing the log-likelihoodusing Eq. (40.36) or in the case where regularization is required one can maximize ϕ(µ, θ) =ln L(µ, θ) + τ S(µ), where the regularization function S(µ) reflects the smoothness of the true
histogram and serves to reduce the variances of the estimators. Possible functions are based on
the mean squared second derivative (a commonly used type of Tikhonov regularization) or theentropy of the true histogram. The parameter τ fixes the relative weighting of the log-likelihood
and the regularization function and thus determines the balance between the bias and varianceof the resulting estimators. Further discussion of the unfolding problem including methods forchoosing the regularization function and parameter, as well as techniques that employ other types ofregularization such as the iterative Bayes (Richardson-Lucy) method, can be found in Refs. [8,11–13]
and references therein.
40.2.6 The Bayesian approach
In the frequentist methods discussed above, probability is associated only with data, not withthe value of a parameter. This is no longer the case in Bayesian statistics, however, which we
introduce in this section. For general introductions to Bayesian statistics see, e.g., Refs. [14–17].
Suppose the outcome of an experiment is characterized by a vector of data x, whose probability
distribution depends on an unknown parameter (or parameters) θ that we wish to determine. In
Bayesian statistics, all knowledge about θ is summarized by the posterior p.d.f. p(θ|x), whose
integral over any given region gives the degree of belief for θ to take on values in that region, given
the data x. It is obtained by using Bayes’ theorem,p(θ|x) = R

P (x|θ)π(θ),
P (x|θ 0 )π(θ 0 ) dθ 0

31st May, 2024

(40.37)

11

40. Statistics

where P (x|θ) is the likelihood function, i.e., the joint p.d.f. for the data viewed as a function of
θ, evaluated with the data actually obtained in the experiment, and π(θ) is the prior p.d.f. for θ.
Note that the denominator in Eq. (40.37) serves to normalize the posterior p.d.f. to unity.
As it can be difficult to report the full posterior p.d.f. p(θ|x), one would usually summarize itwith statistics such as the mean, mode or median value and covariance matrix. In addition one
may construct intervals with a given probability content, as is discussed in Sec. 40.4.1 on Bayesianinterval estimation.
40.2.6.1 Priors
Bayesian statistics supplies no unique rule for determining the prior π(θ); this reflects theanalyst’s subjective degree of belief (or state of knowledge) about θ before the measurement was
carried out. For the result to be of value to the broader community, whose members may not sharethese beliefs, it is important to carry out a sensitivity analysis, that is, to show how the resultchanges under a reasonable variation of the prior probabilities.
One might like to construct π(θ) to represent complete ignorance about the parameters bysetting it equal to a constant. A problem here is that if the prior p.d.f. is flat in θ, then it is notflat for a nonlinear function of θ, and so a different parametrization of the problem would lead in
general to a non-equivalent posterior p.d.f.
For the special case of a constant prior, one can see from Bayes’ theorem (40.37) that theposterior is proportional to the likelihood, and therefore the mode (peak position) of the posterior
is equal to the maximum-likelihood estimator. The posterior mode, however, will change in general
upon a transformation of parameter. One may use as the Bayesian estimator a summary statistic
other than the mode, such as the median, which is invariant under parameter transformation. But
this will not in general coincide with the MLE.
The difficult and subjective nature of encoding personal knowledge into priors has led to what
is called objective Bayesian statistics, where prior probabilities are based not on an actual degreeof belief but rather derived from formal rules. These give, for example, priors which are invariant
under a transformation of parameters, or ones which result in a maximum gain in information fora given set of measurements. For an extensive review see, e.g., Ref. [18].
Objective priors do not in general reflect degree of belief, but they could in some cases betaken as possible, although perhaps extreme, subjective priors. The posterior probabilities as welltherefore do not necessarily reflect a degree of belief. However one may regard investigating avariety of objective priors to be an important part of the sensitivity analysis. Furthermore, use ofobjective priors with Bayes’ theorem can be viewed as a recipe for producing estimators or intervals
which have desirable frequentist properties.
An important procedure for deriving objective priors is due to Jeffreys. According to Jeffreys’
rule one takes the prior asqπ(θ) ∝ det(I(θ)) ,(40.38)where"∂ 2 ln P (x|θ)
Iij (θ) = −E
∂θi ∂θj#(40.39)is the Fisher information matrix. One can show that the Jeffreys prior leads to inference that isinvariant under a transformation of parameters. One should note that the Jeffreys prior does notin general correspond to one’s degree of belief about the value of a parameter. As examples, theJeffreys prior for the mean µ of a Gaussian distribution is a constant, and for the mean of a Poisson√
distribution one finds π(µ) ∝ 1/ µ.
√
Neither the constant nor 1/ µ priors can be normalized to unit area and are therefore said tobe improper. This can be allowed because the prior always appears multiplied by the likelihood
31st May, 2024

1240. Statistics

function, and if the likelihood falls to zero sufficiently quickly then one may have a normalizable
posterior density.
An important type of objective prior is the reference prior due to Bernardo and Berger [19]. To
find the reference prior for a given problem one considers the Kullback-Leibler divergence Dn [π, p]of the posterior p(θ|x) relative to a prior π(θ), obtained from a set of i.i.d. data x = (x1 , . . . , xn ):Z

Dn [π, p] =p(θ|x) lnp(θ|x)dθ .
π(θ)(40.40)This is effectively a measure of the gain in information provided by the data. The reference prior ischosen so that the expectation value of this information gain is maximized for the limiting case ofn → ∞, where the expectation is computed with respect to the marginal distribution of the data,Zp(x) =

p(x|θ)π(θ) dθ .

(40.41)

For a single, continuous parameter the reference prior is usually identical to the Jeffreys prior.
In the multiparameter case an iterative algorithm exists, which requires sorting the parameters byorder of inferential importance. Often the result does not depend on this order, but when it does,
this can be part of a sensitivity analysis. Further discussion and applications to particle physics
problems can be found in Ref. [20].
40.2.6.2 Bayesian treatment of nuisance parameters
As discussed in Sec. 40.2.2, a model may depend on parameters of interest θ as well as on
nuisance parameters ν, which must be included for an accurate description of the data. Knowledge about the values of ν may be supplied by control measurements, theoretical insights, physicalconstraints, etc. Suppose, for example, one has data y from a control measurement which is characterized by a probability p(y|ν). Suppose further that before carrying out the control measurement
one’s state of knowledge about ν is described by an initial prior π0 (ν), which in practice is often
taken to be a constant or in any case very broad. By using Bayes’ theorem (40.1) one obtains the
updated prior π(ν) (i.e., now π(ν) = π(ν|y), the probability for ν given y),π(ν|y) ∝ p(y|ν)π0 (ν) .

(40.42)In the absence of a model for P (y|ν) one may make some reasonable but ad hoc choices inorder to approximate π(ν). For a single nuisance parameter ν, for example, one might characterize
the uncertainty by a p.d.f. π(ν) centered about its nominal value with a certain standard deviationσν . Often a Gaussian p.d.f. provides a reasonable model for one’s degree of belief about a nuisanceparameter; in other cases, more complicated shapes may be appropriate. If, for example, theparameter represents a non-negative quantity then a log-normal or gamma p.d.f. can be a morenatural choice than a Gaussian truncated at zero. Note also that truncation of the prior of a
nuisance parameter ν at zero will in general make π(ν) nonzero at ν = 0, which can lead to an
unnormalizable posterior for a parameter of interest that appears multiplied by ν.
The likelihood function, prior, and posterior p.d.f.s all depend on both θ and ν, and are relatedby Bayes’ theorem, as usual. Note that the likelihood here only refers to the primary measurement
x. Once any control measurements y are used to find the updated prior π(ν) for the nuisanceparameters, this information is fully encapsulated in π(ν) and the control measurements do notappear further.
One can obtain the posterior p.d.f. for θ alone by integrating over the nuisance parameters, i.e.,Zp(θ|x) =p(θ, ν|x) dν .

31st May, 2024

(40.43)

13

40. Statistics

Such integrals can often not be carried out in closed form, and if the number of nuisance parameters
is large, then they can be difficult to compute with standard Monte Carlo methods. Markov Chain
Monte Carlo (MCMC) techniques are often used for computing integrals of this type (see Sec. 42.6).

40.3 Statistical tests
In addition to estimating parameters, one often wants to assess the validity of certain statements
concerning the data’s underlying distribution. Frequentist hypothesis tests, described in Sec. 40.3.1,
provide a rule for accepting or rejecting hypotheses depending on the outcome of a measurement. In
significance tests, covered in Sec. 40.3.2, one gives the probability to obtain a level of incompatibility
with a certain hypothesis that is greater than or equal to the level observed with the actual data.
Goodness-of-fit tests that quantify the general level of compatibility between data and a hypothesis
are described in Sec. 40.3.3. In the Bayesian approach, the corresponding procedure is based
fundamentally on the posterior probabilities of the competing hypotheses. In Sec. 40.3.4 we describea related construct called the Bayes factor, which can be used to quantify the degree to which thedata prefer one or another hypothesis.
40.3.1 Hypothesis tests
A frequentist test of a hypothesis (often called the null hypothesis, H0 ) is a rule that states
for which data values x the hypothesis is rejected. A region of x-space called the critical region,w, is specified such that there is no more than a given probability under H0 , α, called the size or
significance level of the test, to find x ∈ w. If the data are discrete, it may not be possible to finda critical region with exact probability content α, and thus we require P (x ∈ w|H0 ) ≤ α. If thedata are observed in the critical region, H0 is rejected.
The data x used to construct a test could be, for example, a set of values that characterizes anindividual event. In this case the test corresponds to classification as, e.g., signal or background.
Alternatively the data could represent a set of values from a collection of events. Often one isinterested in knowing whether all of the events are of a certain type (background), or whether the
sample contains at least some events of a new type (signal). Here the background-only hypothesis
plays the role of H0 , and in the alternative H1 both signal and background are present. Rejecting
H0 is, from the standpoint of frequentist statistics, the required step to establish discovery of thesignal process.
The critical region is not unique. Its choice should take into account the probabilities for thedata predicted by some alternative hypothesis (or set of alternatives) H1 . Rejecting H0 if it istrue is called a type-I error, and occurs by construction with probability no greater than α. Not
rejecting H0 if an alternative H1 is true is called a type-II error, and for a given test this will havea certain probability β = P (x ∈/ w|H1 ). The quantity 1 − β is called the power of the test of H0
with respect to the alternative H1 . A strategy for defining the critical region can therefore be to
maximize the power with respect to some alternative (or alternatives) given a fixed size α.
To maximize the power of a test of H0 with respect to the alternative H1 , the Neyman–Pearsonlemma states that the critical region w should be chosen such that for all data values x inside w,
the likelihood ratio
f (x|H1 )λ(x) =(40.44)f (x|H0 )is greater than or equal to a given constant cα , and everywhere outside the critical region one has
λ(x) < cα , where the value of cα is determined by the size of the test α. The hypotheses H0 andH1 must be simple, i.e., they should not contain undetermined parameters.
It is convenient to define the test using a scalar function of the data x called a test statistic,t(x), such that the boundary of the critical region is given by a surface of constant t(x). The
Neyman–Pearson lemma is equivalent to the statement that the likelihood ratio (40.44) represents
31st May, 2024

14

40. Statistics

the optimal test statistic. It can be difficult in practice, however, to determine λ(x), since thisrequires knowledge of the joint p.d.f.s f (x|H0 ) and f (x|H1 ). Often one does not have explicit
formulae for these, but rather Monte Carlo models that allow one to generate instances of x that
follow the p.d.f.s.
In the case where the likelihood ratio (40.44) cannot be used explicitly, there exist a variety
of other multivariate methods for constructing a test statistic that may approach its performance.
These are based on machine-learning algorithms that use samples of training data corresponding
to the hypotheses in question, often generated from Monte Carlo models. Further information on
Machine Learning can be found in Sec. 41 of this Review.
The multivariate algorithms designed to classify events into signal and background types also
form the basis of tests of the hypothesis that a sample of events consists of background only. Sucha test can be constructed using the distributions of the test statistic t(x) for event classification
obtained from a multivariate algorithm such as a Neural Network output. The distributions p(t|s)
and p(t|b) for signal and background events, respectively, are used to construct the likelihood ratioof the signal-plus-background hypothesis relative to that of background only. To the extent that thetest statistic t(x) approximates the likelihood ratio (or a monotonic function thereof) for individual
events given by (40.44), the resulting test of the background-only hypothesis for the event sample
will have maximum power with respect to the signal-plus-background alternative (see Ref. [21]).
40.3.2 Tests of significance, p-values
A frequentist significance test allows one to assign a numerical value, called the p-value, thatreflects the level of agreement between a hypothesis H0 and the observed data x. To carry out asignificance test one must specify a subset of the data space, whose boundary includes the observed
data point xobs , that is deemed to have equal or lesser compatibility with H0 relative to xobs . The
p-value of H0 is the probability, assuming data distributed according to H0 , to find x in the regionof equal or lesser compatibility.
For continuous data, one can show that the p-value of H0 , here called p0 , follows a uniform
distribution for data generated under H0 . Therefore one has P (p0 ≤ α|H0 ) = α for any constant
α in [0, 1]. In this way the region of data space where p0 ≤ α gives a critical region w of a test of
size α, as it satisfies P (x ∈ w|H0 ) ≤ α. Thus rejecting H0 if p0 ≤ α is equivalent to the hypothesis
test of Sec. 40.3.1 above. For the critical region w defined by p0 ≤ α one can find the power
P (x ∈ w|H1 ) with respect to an alternative hypothesis H1 , as defined in Sec. 40.3.1. In general theregion of low compatibility with H0 is defined so as to have a high probability for data generated
under H1 , so that the distribution of p0 , assuming data x that follows H1 , will be concentrated at
low values, and the critical region w derived from p0 ≤ α will have a high power with respect toH1 .
Similar to construction of a test’s critical region in Sec. 40.3.1, finding the p-value is generally
done by defining a statistic t(x) whose value corresponds to better or worse levels of agreementwith the hypothesis. The hypothesis being tested, H0 , will determine the p.d.f. f (t|H0 ) for thestatistic. If, for example, t is defined such that large values correspond to poor agreement with the
hypothesis, then the p-value isZ∞p=f (t|H0 ) dt ,(40.45)tobs

where tobs = t(xobs ) is the value of the statistic obtained in the actual experiment.
Note that the p-value is not the probability of the hypothesis, which in frequentist statistics is
not defined. The p-value should also not be confused with the size (significance level) of a test, or
the confidence level of a confidence interval (Section 40.4), both of which are pre-specified constants.
When searching for a new phenomenon, one tries to reject the hypothesis H0 that the data
31st May, 2024

15

40. Statistics

are consistent with known (e.g., Standard Model) processes. If the p-value of H0 is sufficiently
low, then one is willing to accept that some alternative hypothesis is true. Often one converts thep-value into an effective significance Z, defined as an equivalent number of standard deviations ofa Gaussian distributed random variable. In a search for an intrinsically positive signal, i.e., where
only upward fluctuations of the estimated rate appear signal like, this is defined as
Z = Φ−1 (1 − p) .

(40.46)

Here Φ is the cumulative distribution of the standard Gaussian, and Φ−1 is its inverse (quantile)function. In this way, a p-value of 1/2 gives Z = 0, i.e., having an estimated signal rate of zerocorresponds to zero significance. If either positive or negative data fluctuations would indicate
evidence of the signal, then one defines
Z = Φ−1 (1 − p/2) .

(40.47)In this case an estimated signal rate of zero gives p = 1 and Z = 0.
Often in particle physics the level of significance where an effect is said to qualify as a discovery
is Z = 5, i.e., a 5σ effect, corresponding to a p-value of 2.87 × 10−7 . One’s actual degree of belief
that a new process is present, however, will depend in general on other factors as well, such as
the plausibility of the new signal hypothesis and the degree to which it can describe the data,
one’s confidence in the model that led to the observed p-value, and possible corrections for multiple
observations out of which one focuses on the smallest p-value obtained (the “look-elsewhere effect”,
discussed in Section 40.3.2.2).
40.3.2.1 Frequentist treatment of nuisance parameters and asymptotic methods for tests
Suppose one wants to test hypothetical values of parameters θ, but the model also contains
nuisance parameters ν. To find a p-value for θ we can construct a test statistic qθ such that larger
values constitute increasing incompatibility between the data and the hypothesis. Then for an
observed value of the statistic qθ,obs , the p-value of θ isZ ∞pθ (ν) =qθ,obsf (qθ |θ, ν) dqθ ,(40.48)which depends in general on the nuisance parameters ν. In the strict frequentist approach, θ isrejected only if the p-value is less than α for all possible values of the nuisance parameters.
The difficulty described above is effectively solved if we can define the test statistic qθ in sucha way that its distribution f (qθ |θ) is independent of the nuisance parameters. Although exact
independence is only found in special cases, it can be achieved approximately by use of the profile
likelihood ratio. This is given by the profile likelihood from Eq.(40.18) divided by the value of theblikelihood at its maximum, i.e., when evaluated with the maximum-likelihood estimators θb and ν:
λp (θ) =bb
L(θ, ν(θ))
b ν)
b
L(θ,

.

(40.49)Wilks’ theorem [10] states that, providing certain general conditions are satisfied, the distributionof −2 ln λp (θ), under assumption of θ, approaches a χ2 distribution in the limit where the data
sample is very large, independent of the values of the nuisance parameters ν. Here the number ofdegrees of freedom is equal to the number of components of θ. More details on use of the profile
likelihood are given in Refs. [22,23] and in contributions to the PHYSTAT conferences [24]; explicit
formulae for special cases can be found in Ref. [25]. Further discussion on how to incorporate
systematic uncertainties into p-values can be found in Ref. [26].
31st May, 2024

16

40. Statistics

Even with use of the profile likelihood ratio, for a finite data sample the p-value of hypothesizedparameters θ will retain in general some dependence on the nuisance parameters ν. Ideally one
would find the the maximum of pθ (ν) from Eq. (40.48) explicitly, but that is often impractical. Anbbbb
are the profiled
where ν(θ)
approximate and computationally feasible technique is to use pθ (ν(θ)),
values of the nuisance parameters as defined in Section 40.2.2.2. The resulting p-value is correct if
the true values of the nuisance parameters are equal to the profiled values used; otherwise it could
be either too high or too low. This is discussed further in Section 40.4.2 on confidence intervals.
The methods above based on the profile likelihood ratio are useful in practice because, provided
one has a sufficiently large data sample, the distributions of the test statistics are related throughWilks’ theorem to the chi-square p.d.f. If the data sample is not large enough to justify use of
the asymptotic formulae, then some adjustments can be made to the test statistics so that theirdistributions are better approximated by the asymptotic ones; these improvements go under thegeneral name of “higher-order asymptotics”. Using these methods one may obtain p-values using
the asymptotic formulae even for smaller data samples, avoiding costly Monte Carlo simulation tofind distributions of test statistics.
Roughly speaking, the chi-square distribution for the profile likelihood ratio relies on having a
Gaussian distribution for estimators of the model’s parameters. An important type of corrected
statistic relies on an improved distribution for the maximum-likelihood estimator due to BarndorffNielsen [27] (the p∗ approximation). Another class of improved statistic due to Bartlett [28, 29]
starts with a statistic t which, asymptotically, should be chi-square distributed with nd degrees of
freedom. A new statistic is defined as t0 = tnd /E[t], which by construction has mean E[t0 ] = nd ,as for the chi-square distribution. Further details and applications of these methods are describedin Refs. [30, 31].
One may also treat model uncertainties in a Bayesian manner but then use the resulting modelin a frequentist test. Suppose the uncertainty in a set of nuisance parameters ν is characterized
by a Bayesian prior p.d.f. π(ν). This can be used to construct the marginal (also called the priorpredictive) model for the data x and parameters of interest θ,Z

Pm (x|θ) =P (x|θ, ν)π(ν) dν .

(40.50)The marginal model does not represent the probability of data that would be generated if one werereally to repeat the experiment, as in that case one would assume that the nuisance parameters
do not vary. Rather, the marginal model represents a situation in which every repetition of theexperiment is carried out with new values of ν, randomly sampled from π(ν). It is in effect anaverage of models each with a given ν, where the average is carried out with respect to the priorp.d.f. π(ν).
The marginal model for the data x can be used to determine the distribution of a test statistic
Q, which can be writtenZ
Pm (Q|θ) =P (Q|θ, ν)π(ν) dν .

(40.51)

In a search for a new signal process, the test statistic can be based on the ratio of likelihoods
corresponding to the experiments where signal and background events are both present, Ls+b , tothat of background only, Lb . Often the likelihoods are evaluated with the profiled values of thenuisance parameters, which may give improved performance. It is important to note, however, that
it is through use of the marginal model for the distribution of Q that the uncertainties related to
the nuisance parameters are incorporated into the result of the test. Different choices for the test
statistic itself only result in variations of the power of the test with respect to different alternatives.
31st May, 2024

17

40. Statistics

Studies of marginalization versus profiling of nuisance parameters for specific problems, e.g., related
to a Poisson counting experiment, can be found in Refs. [32, 33].
40.3.2.2 The look-elsewhere effect
The “look-elsewhere effect” relates to multiple measurements used to test a single hypothesis.
The classic example is when one searches in a distribution for a peak whose position is not predicted
in advance. Here the no-peak hypothesis is tested using data in a given range of the distribution. In
the frequentist approach the correct p-value of the no-peak hypothesis is the probability, assuming
background only, to find a signal as significant as the one found or more so anywhere in the searchregion. This can be substantially higher than the probability to find a peak of equal or greater
significance in the particular place where it appeared. There is in general some ambiguity as to what
constitutes the relevant search region or even the broader set of relevant measurements. Although
the desired p-value is well defined once the search region has been fixed, an exact treatment canrequire extensive computation.
The “brute-force” solution to this problem by Monte Carlo involves generating data under the
background-only hypothesis and for each data set, fitting a peak of unknown position and recording
a measure of its significance. To establish a discovery one often requires a p-value smaller than
2.87×10−7 , corresponding to a 5σ or larger effect. Determining this with Monte Carlo thus requires
generating and fitting a very large number of experiments, perhaps several times 107 . In contrast,
if the position of the peak is fixed, then the fit to the distribution is much easier, and furthermore
one can in many cases use formulae valid for sufficiently large samples that bypass completely theneed for Monte Carlo (see, e.g., Ref. [25]). However, this fixed-position or “local” p-value would
not be correct in general, as it assumes the position of the peak was known in advance.
A method that allows one to modify the local p-value computed under assumption of a fixed
position to obtain an approximation to the correct “global” value using a relatively simple calculation is described in Ref. [34]. Suppose a model contains a nuisance parameter such as thepeak position that is only defined under the signal model (there is no peak in the background-onlymodel). Furthermore, suppose a test statistic q0 is defined using the profile likelihood ratio, so thatby Wilk’s theorem it would be asymptotically chi-square distributed if the peak position were to
be fixed. The asymptotic distribution no longer holds if the peak position is adjustable, however,as this violates the regularity conditions of Wilks’ theorem. If the statistic has an observed valueu, then an approximation for the global p-value is found to be
pglobal ≈ plocal + hNu i ,(40.52)where hNu i, which is much smaller than one in cases of interest, is the mean number of “upcrossings” (as defined in Ref. [34]) of the statistic q0 above the level u in the range of the nuisanceparameter considered (e.g., the mass range). The value of hNu i can be estimated from the numberof upcrossings hNu0 i above some much lower value, u0 , by using a relation due to Davis [35],hNu i ≈ hNu0 ie−(u−u0 )/2 .

(40.53)

By choosing u0 sufficiently low, the value of hNu i can be estimated by simulating only a very small
number of experiments, or even from the observed data, rather than of the order 107 needed if one
is dealing with a 5σ effect.
40.3.3 Goodness of fit
At times one wants to quantify the level of agreement between the data x and a hypothesis H0
without explicit reference to alternative hypotheses. In the frequentist approach this is done using
a goodness-of-fit test. The result is quantified with a p-value, in general obtained from a statistic
31st May, 2024

40. Statistics

18

t(x) as done in Sec. 40.3.2 with a significance test. Here, however, the goodness-of-fit statistic isdefined using general considerations of what constitutes greater or lesser compatibility between
data and hypothesis and not with reference to alternative hypotheses or power. Nevertheless, fora given statistic t(x), the power of the resulting test with respect to an alternative can be found,so that types of alternatives to which one is sensitive (i.e., to which the test has high power) canbe investigated.
The Neyman-Pearson lemma (see Sec. 40.3.1) states that a test of H0 that has maximum powerwith respect to an alternative H1 should be based on a statistic that is monotonic in the likelihood
ratio P (x|H1 )/P (x|H0 ). Therefore a goodness-of-fit statistic derived without reference to a specific
H1 will not in general attain this maximum power. And even if a given t(x) happens to be optimalor nearly so for some alternative H1 , there may be other relevant alternatives for which its power islower. In practice, goodness-of-fit statistics are defined so as to provide sensitivity to a broad classof alternatives, e.g., those that differ in the location of the bulk of the distribution’s probabilityor the distribution’s width or the nature of the tails. It may therefore be useful to use several
goodness-of-fit tests that are sensitive to different aspects of the data distribution in question.
The most important application of goodness-of-fit tests is no doubt in conjunction with themethod of least squares, which we describe in Sec. 40.3.3.1 below. Other types of tests, particularly
those based on the empirical distribution function are described in Sec. 40.3.3.2. Further discussionof goodness of fit can be found in Refs. [4, 36] and references therein.
40.3.3.1 Goodness-of-fit with chi-square or likelihood ratio
An important type of goodness-of-fit test arises when a set of measurements y = (y1 , . . . , yN )are used to test the hypothesis that the values are Gaussian distributed with given mean values
E[yi ] = µi and variances V [yi ] = σi2 (or covariance Vij = cov[yi , yj ] if the values are correlated).
The general procedures also hold when the data are Poisson or multinomially distributed with
sufficiently large mean values so that the Gaussian represents an adequate approximation. The
Gaussian hypothesis may be simple (all parameters fully specified), or the mean values µi may be
given as functions of some other parameters θ = (θ1 , . . . , θM ), with M < N . The (co)variances aretreated as fixed.
The test statistic derived from the likelihood ratio has already been introduced in Sec. 40.2.3
in connection with the method of least squares, and is given by the minimum of the function χ2 (θ)
with respect to θ, i.e., t = χ2min = χ2 (θ̂), e.g., in the case of uncorrelated measurementst=NX(yi − µi (θ̂))2
i=1

σi2

.

(40.54)For Gaussian distributed y, the statistic t follows a chi-square distribution for N − M degrees
of freedom. From Eq. (40.54) one sees that greater t corresponds to increasing incompatibility
between the measured yi and predicted µi . The p-value of the Gaussian hypothesis is therefore
taken to be
Z ∞

p=

χ2min

f (t; nd ) dt ,(40.55)where f (t; nd ) is the χ2 p.d.f. and nd = N − M is the appropriate number of degrees of freedom.
Values are shown in Fig. 40.1 or obtained from standard computer libraries. If the asymptotic
conditions for using the χ2 p.d.f. do not hold, the statistic can still be defined as before, but its
p.d.f. must be determined by other means in order to obtain the p-value, e.g., using a Monte Carlocalculation.
31st May, 2024

40. Statistics

19

Often the data represent numbers of events in N bins of a histogram, i.e., n = (n1 , . . . , nN ). An
P
important case is when n follows a multinomial distribution with ntot = N
i=1 ni total entries andmean values µi , or equivalently bin probabilities pi = µi /ntot . Alternatively one may model the ni
as independent and Poisson distributed with means µi . If there are no adjustable parameters in
the hypotheses, then the goodness-of-fit can be quantified with Pearson’s χ2 statistic, defined for
multinomial data asNX(ni − ntot pi )2,(40.56)tM =ntot pi
i=1
or for the independent Poisson data as
tP =

N
X
(ni − µi )2
i=1

µi

.

(40.57)In the limit where the means µi are large, the statistics tM and tP are found to follow chi-square
distributions for N −1 and N degrees of freedom, respectively. Having one fewer degrees of freedom
in the multinomial case is related to the fact that ntot is fixed and thus the ni are correlated.
Alternatively one can use the test statistic based on the likelihood ratio t = −2 ln λ(θ̂) given
previously in Eq. (40.16). Here the likelihood ratio λ is chosen to correspond to the multinomial or
Poisson model as appropriate, as described in Sec. 40.2.2.1. One finds that the distribution of thelikelihood-ratio statistic approaches the asymptotic limit faster than does Pearson’s chi-square andthus when using the chi-square p.d.f. to compute a p-value one obtains in general a more accurate
result from −2 ln λ (see [9]).
Since the mean of the chi-square distribution is equal to nd , one expects in a “reasonable”
experiment to obtain χ2 ≈ nd (here χ2 refers to the minimized value). Hence the quantity χ2 /nd
is sometimes reported. Since the p.d.f. of χ2 /nd depends on nd , however, one must report nd as
well if one wishes to determine the p-value. The p-values obtained for different values of χ2 /nd areshown in Fig. 40.2.
If the minimized χ2 value indicates a low level of agreement between data and hypothesis, one
may be tempted to expect a high degree of uncertainty for any fitted parameters. Poor goodness-offit, however, does not mean that one will have large statistical errors for parameter estimates. If, forexample, the error bars (or covariance matrix) used in constructing the χ2 are underestimated, thenthis will lead to underestimated statistical errors for the fitted parameters and an increased valueof the minimized χ2 . The standard deviations of estimators that one finds from, say, Eq. (40.13)
reflect how widely the estimates would be distributed if one were to repeat the measurement manytimes, assuming that the hypothesis and measurement errors used in the χ2 are also correct. They
do not include the systematic error which may result from an incorrect hypothesis or incorrectly
estimated measurement errors in the χ2 .
40.3.3.2 Goodness-of-fit with the empirical distribution function
Suppose a measurement yields a sample of independent and identically distributed (i.i.d.) values
x = (x1 , . . . , xn ). Often the values may be summarized by creating a histogram, but this inevitably
results in a loss of information, since the position of the x values within each bin is not retained.
Particularly with small numbers of observations n one therefore may prefer to keep the values ofeach of the xi (i.e., “unbinned”). The sample may be displayed graphically using the empirical
distribution function (e.d.f.)Fn (x) =number of elements in sample with value ≤ x.
n31st May, 2024

(40.58)

40. Statistics

20

1.000

p-value for test
α for confidence intervals

0.500
0.200

n=1

2

3 4

15

6 8

25

40

0.100
10

0.050

20 30

50

0.020
0.010
0.005
0.002
0.001
1

2

3 4 5

7 10

20

30 40 50 70 100

χ2Figure 40.1: One minus the χ2 cumulative distribution, 1 − F (χ2 ; n), for n degrees of freedom.
This gives the p-value for the χ2 goodness-of-fit test as well as one minus the coverage probability
for confidence regions (see Sec. 40.4.2.2).
The e.d.f. Fn (x) begins at zero for x less than the lowest observed value, increases by 1/n at eachobserved xi , and saturates at unity for x equal to the highest value in the sample.
Given an i.i.d. data sample x = (x1 , . . . , xn ), suppose one wants to test the hypothesis H0 thatx follows the p.d.f. f (x) using a goodness-of-fit test, i.e., without reference to a specific alternative.
The hypothesisof the p.d.f. f (x) is equivalent to saying the cumulative distribution function (c.d.f.)
Rx
is F (x) = −∞ f (x0 ) dx0 . As the e.d.f. Fn (x) can be viewed as an estimator for the c.d.f. F (x), we canderive a goodness-of-fit statistic by some appropriately defined measure of the “distance” between
the two.
A widely used example is the Kolmogorov-Smirnov (K-S) test, based on
Dn = sup |Fn (x) − F (x)| ,(40.59)xwhere sup denotes the supremum. That is, Dn is the greatest vertical distance between Fn (x)and F (x) for any x. In the limit of large sample size n, one can show that the distribution of
Dn , assuming the data follows the c.d.f. F (x), approaches the Kolmogorov distribution, which is
independent of the specific form of F (x) and is provided in many software packages. A larger valueof Dn corresponds to greater incompatibility between data and hypothesis, so the p-value is theprobability P (Dn ≥ Dn,obs |H0 ). The K-S test is sensitive in particular to differences in location(i.e., in the mean of x), as this corresponds to a horizontal shift between Fn (x) and F (x), giving alarge maximum vertical distance Dn .
The Cramér-von Mises family of tests is based on the squared difference between Fn (x) andF (x) through an integral of the form
31st May, 2024

40. Statistics

21

2.5

2.0
1%
1.5

χ2/n
1.0

50%

0.5

90%
99%

0.0

5%

10%

0

10

32%
68%
95%

20
30
Degrees of freedom n

40

50

Figure 40.2: The ‘reduced’ χ2 , equal to χ2 /n, for n degrees of freedom. The curves show as afunction of n the χ2 /n that corresponds to a given p-value.

2

Z ∞

W =
−∞

(F (x) − Fn (x))2 w(x)f (x) dx ,(40.60)where f (x) = dF/dx is the p.d.f. of x and w(x) is a weight function that must be specified. Two
important choices are the classic Cramér-von Mises test with w(x) = 1, and the Anderson-Darling
test with w(x) = [F (x)(1 − F (x))]−1 , which gives increased sensitivity to departures in the tailsof the distribution. For cases such as these where the weight function depends on x only throughthe cumulative distribution F (x), one can show (see, e.g., Refs. [4, 36]) that the distribution of
W 2 , assuming data that follows F (x), is independent of the form of F (x). The distributions of
the Smirnov-Cramér-von Mises and Anderson-Darling statistics are available in standard computer
libraries.
40.3.4 Bayes factors
In Bayesian statistics, all of one’s knowledge about a model is contained in its posterior probability, which one obtains using Bayes’ theorem (Eq. (40.37)). Thus one could reject a hypothesis
H if its posterior probability P (H|x) is sufficiently small. The difficulty here is that P (H|x) is
proportional to the prior probability P (H), and there will not be a consensus about the prior probabilities for the existence of new phenomena. Nevertheless one can construct a quantity called theBayes factor (described below), which can be used to quantify the degree to which the data prefer
one hypothesis over another, and is independent of their prior probabilities.
Consider two models (hypotheses), Hi and Hj , described by vectors of parameters θi and θj ,31st May, 2024

40. Statistics

22

respectively. Some of the components will be common to both models and others may be distinct.
The full prior probability for each model can be written in the formπ(Hi , θi ) = P (Hi )π(θi |Hi ) .

(40.61)

Here P (Hi ) is the overall prior probability for Hi , and π(θi |Hi ) is the normalized p.d.f. of its
parameters. For each model, the posterior probability is found using Bayes’ theorem,R

P (Hi |x) =P (x|θi , Hi )P (Hi )π(θi |Hi ) dθi,
P (x)(40.62)where the integration is carried out over the internal parameters θi of the model. The ratio of
posterior probabilities for the models is thereforeR

P (Hi |x)P (x|θi , Hi )π(θi |Hi ) dθi P (Hi )=R.
P (Hj |x)P (x|θj , Hj )π(θj |Hj ) dθj P (Hj )(40.63)

The Bayes factor is defined asR
P (x|θi , Hi )π(θi |Hi ) dθi.
Bij = RP (x|θj , Hj )π(θj |Hj ) dθj(40.64)This gives what the ratio of posterior probabilities for models i and j would be if the overall priorprobabilities for the two models were equal. If the models have no nuisance parameters, i.e., nointernal parameters described by priors, then the Bayes factor is simply the likelihood ratio. The
Bayes factor therefore shows by how much the probability ratio of model i to model j changes in
the light of the data, and thus can be viewed as a numerical measure of evidence supplied by thedata in favor of one hypothesis over the other.
Although the Bayes factor is by construction independent of the overall prior probabilities
P (Hi ) and P (Hj ), it does require priors for all internal parameters of a model, i.e., one needs
the functions π(θi |Hi ) and π(θj |Hj ). In a Bayesian analysis where one is only interested in theposterior p.d.f. of a parameter, it may be acceptable to take an unnormalizable function for the prior(an improper prior) as long as the product of likelihood and prior can be normalized. Improper
priors are, however, only defined up to an arbitrary multiplicative constant, and so the Bayes factor
would depend on this constant. Furthermore, although the range of a constant normalized prior is
unimportant for parameter determination (provided it is wider than the likelihood), this is not sofor the Bayes factor when such a prior is used for only one of the hypotheses. So to compute a Bayes
factor, all internal parameters must be described by normalized priors that represent meaningfulprobabilities over the entire range where they are defined.
An exception to this rule may be considered when the identical parameter appears in themodels for both numerator and denominator of the Bayes factor. In this case one can argue thatthe arbitrary constants would cancel. One must exercise some caution, however, as parameterswith the same name and physical meaning may still play different roles in the two models.
Both integrals in Equation (40.64) are of the form
Zm=

P (x|θ)π(θ) dθ ,(40.65)which is similar to the marginal likelihood seen previously in Eq. (40.50) (in some fields this quantity
is called the evidence). Computing marginal likelihoods can be difficult; in many cases it can bedone with the nested sampling algorithm [37] as implemented, e.g., in the program MultiNest [38].
A review of Bayes factors can be found in Ref. [39].
31st May, 2024

40. Statistics

23

40.4 Intervals and limits
When the goal of an experiment is to determine a parameter θ, the result is usually expressedby quoting, in addition to the point estimate, some sort of interval which reflects the statisticalprecision of the measurement. In the simplest case, this can be given by the parameter’s estimatedb σ
bb. If, however, the p.d.f. of thevalue θb plus or minus an estimate of the standard deviation of θ,θestimator is not Gaussian or if there are physical boundaries on the possible values of the parameter,
then one usually quotes instead an interval according to one of the procedures described below.
In reporting an interval or limit, the experimenter may wish to• communicate as objectively as possible the result of the experiment;• provide an interval that is constructed to cover on average the true value of the parameterwith a specified probability;• provide the information needed by the consumer of the result to draw conclusions about theparameter or to make a particular decision;• draw conclusions about the parameter that incorporate stated prior beliefs.
With a sufficiently large data sample, the point estimate and standard deviation (or for themultiparameter case, the parameter estimates and covariance matrix) satisfy essentially all of thesegoals. For small data samples, no single method for quoting an interval will achieve all of them.
In addition to the goals listed above, the choice of method may be influenced by practicalconsiderations such as ease of producing an interval from the results of several measurements. Of
course the experimenter is not restricted to quoting a single interval or limit; one may choose,for example, first to communicate the result with a confidence interval having certain frequentist
properties, and then in addition to draw conclusions about a parameter using a judiciously chosen
subjective Bayesian prior. It is recommended, however, that there be a clear separation between
these two aspects of reporting a result. In the remainder of this section, we assess the extent towhich various types of intervals achieve the goals stated here.
40.4.1 Bayesian intervals
As described in Sec. 40.2.6, a Bayesian posterior probability may be used to determine regionsthat will have a given probability of containing the true value of a parameter. In the single parameter
case, for example, an interval (called a Bayesian or credible interval) [θlo , θup ] can be determinedwhich contains a given fraction 1 − α of the posterior probability, i.e.,
1−α=

Z θup

p(θ|x) dθ .

(40.66)θlo

Sometimes an upper or lower limit is desired, i.e., θlo or θup can be set to a physical boundary or
to plus or minus infinity. In other cases, one might be interested in the set of θ values for whichp(θ|x) is higher than for any θ not belonging to the set, which may constitute a single interval or aset of disjoint regions; these are called highest posterior density (HPD) intervals. Note that HPDintervals are not invariant under a nonlinear transformation of the parameter.
If a parameter is constrained to be non-negative, then the prior p.d.f. can simply be set to zerofor negative values. An important example is the case of a Poisson variable n, which counts signal
events with unknown mean s, as well as background with mean b, assumed known. For the signal
mean s, one often uses the prior(0s<0π(s) =.
(40.67)1s≥0This prior may be regarded as providing an interval whose frequentist properties can be studied,
31st May, 2024

40. Statistics

24

rather than as representing a degree of belief. For example, to obtain an upper limit on s, one mayproceed as follows. The likelihood for s is given by the Poisson distribution for n with mean s + b,
(s + b)n −(s+b)
e,(40.68)n!
along with the prior (40.67) in (40.37) gives the posterior density for s. An upper limit sup at
confidence level (or here, rather, credibility level) 1 − α can be obtained by requiring
P (n|s) =

1−α=

Z sup
−∞

R sup

P (n|s) π(s) ds,−∞ P (n|s) π(s) dsp(s|n)ds = R−∞∞(40.69)where the lower limit of integration is effectively zero because of the cut-off in π(s). By relating the
integrals in Eq. (40.69) to incomplete gamma functions, the solution for the upper limit is foundto be
1
sup = Fχ−1
(40.70)2 [p, 2(n + 1)] − b ,22where Fχ−1
2 is the quantile of the χ distribution (inverse of the cumulative distribution). Here thequantity p is
(40.71)
p = 1 − α 1 − Fχ2 [2b, 2(n + 1)] ,where Fχ2 is the cumulative χ2 distribution. For both Fχ2 and Fχ−1
2 above, the argument 2(n + 1)
gives the number of degrees of freedom. For the special case of b = 0, the limit reduces to
1
sup = Fχ−1
2 (1 − α; 2(n + 1)) .
2(40.72)It happens that for the case of b = 0, the upper limit from Eq. (40.72) coincides numerically withthe frequentist upper limit discussed in Section 40.4.2.3. Values for 1 − α = 0.9 and 0.95 are givenby the values µup in Table 40.3. The frequentist properties of confidence intervals for the Poissonmean found in this way are discussed in Refs. [2] and [40].
As in any Bayesian analysis, it is important to show how the result changes under assumptionof different prior probabilities. For example, one could consider the Jeffreys
prior as described in√Sec. 40.2.6. For this problem one finds the Jeffreys prior π(s) ∝ 1/ s + b for s ≥ 0 and zerootherwise. As with the constant prior, one would not regard this as representing one’s prior beliefs
about s, both because it is improper and also as it depends on b. Rather it is used with Bayes’
theorem to produce an interval whose frequentist properties can be studied.
If the model contains nuisance parameters then these are eliminated by marginalizing, as in
Eq. (40.43), to obtain the p.d.f. for the parameters of interest. For example, if the parameter b in the
Poisson counting problem above were to be characterized by a prior p.d.f. π(b), then
one would firstR
use Bayes’ theorem to find p(s, b|n). This is then marginalized to find p(s|n) = p(s, b|n)π(b) db,
from which one may determine an interval for s. One may not be certain whether to extend a modelby including more nuisance parameters. In this case, a Bayes factor may be used to determine towhat extent the data prefer a model with additional parameters, as described in Section 40.3.4.
40.4.2 Frequentist confidence intervals
The unqualified phrase “confidence intervals” refers to frequentist intervals obtained with aprocedure due to Neyman [41], described below. The boundary of the interval (or in the multiparameter case, region) is given by a specific function of the data, which would fluctuate if one wereto repeat the experiment many times. The coverage probability refers to the fraction of intervals insuch an ensemble that contain the true parameter value. Confidence intervals are constructed so
31st May, 2024

40. Statistics

25

as to have a coverage probability greater than or equal to a given confidence level, regardless of thetrue parameter’s value. It is important to note that in the frequentist approach, such a probability
is not meaningful for a fixed interval. In this section we discuss several techniques for producing
intervals that have, at least approximately, this property of coverage.
40.4.2.1 The Neyman construction for confidence intervals
Consider a p.d.f. f (x; θ) where x represents the outcome of the experiment and θ is the unknown
parameter for which we want to construct a confidence interval. The variable x could (and often
does) represent an estimator for θ. Using f (x; θ), we can find using a pre-defined rule and probability
1 − α for every value of θ, a set of values x1 (θ, α) and x2 (θ, α) such that
Z x2

P (x1 < x < x2 ; θ) =f (x; θ) dx ≥ 1 − α .

(40.73)x1

If x is discrete, the integral is replaced by the corresponding sum. In that case there may not exist
a range of x values whose summed probability is exactly equal to a given value of 1 − α, and one
requires by convention P (x1 < x < x2 ; θ) ≥ 1 − α.
This is illustrated for continuous x in Fig. 40.3: a horizontal line segment [x1 (θ, α), x2 (θ, α)] is
drawn for representative values of θ. The union of such intervals for all values of θ, designated in
the figure as D(α), is known as a confidence belt. Typically the curves x1 (θ, α) and x2 (θ, α) aremonotonic functions of θ, which we assume for this discussion.

parameter θD(α)x2(θ), θ2(x)θ0
x1(θ), θ1(x)x1(θ0)x2(θ0)

Possible experimental values x
Figure 40.3: Construction of the confidence belt (see text).
Upon performing an experiment to measure x and obtaining a value x0 , one draws a verticalline through x0 . The confidence interval for θ is the set of all values of θ for which the corresponding
31st May, 2024

26

40. Statistics

line segment [x1 (θ, α), x2 (θ, α)] is intercepted by this vertical line. Such confidence intervals are
said to have a confidence level (CL) equal to 1 − α.
Now suppose that the true value of θ is θ0 , indicated in the figure. We see from the figure that
θ0 lies between θ1 (x) and θ2 (x) if and only if x lies between x1 (θ0 ) and x2 (θ0 ). The two events thus
have the same probability, and since this is true for any value θ0 , we can drop the subscript 0 andobtain1 − α = P (x1 (θ) < x < x2 (θ)) = P (θ2 (x) < θ < θ1 (x)) .
(40.74)In this probability statement, θ1 (x) and θ2 (x), i.e., the endpoints of the interval, are the randomvariables and θ is an unknown constant. If the experiment were to be repeated a large numberof times, the interval [θ1 , θ2 ] would vary, covering the fixed value θ in a fraction 1 − α of theexperiments.
The condition of coverage in Eq. (40.73) does not determine x1 and x2 uniquely, and additional
criteria are needed. One possibility is to choose central intervals such that the probabilities to findx below x1 and above x2 are each α/2. In other cases, one may want to report only an upper orlower limit, in which case one of P (x ≤ x1 ) or P (x ≥ x2 ) can be set to α and the other to zero.
Another principle based on likelihood ratio ordering for determining which values of x should be
included in the confidence belt is discussed below.
When the observed random variable x is continuous, the coverage probability obtained with
the Neyman construction is 1 − α, regardless of the true value of the parameter. Because of therequirement P (x1 < x < x2 ) ≥ 1−α when x is discrete, one obtains in that case confidence intervals
that include the true parameter with a probability greater than or equal to 1 − α.
An equivalent method of constructing confidence intervals is to consider a test (see Sec. 40.3)of the hypothesis that the parameter’s true value is θ (assume one constructs a test for all physicalvalues of θ). One then excludes all values of θ where the hypothesis would be rejected in a test of
size α or less. The remaining values constitute the confidence interval at confidence level 1 − α. If
the critical region of the test is characterized by having a p-value pθ ≤ α, then the endpoints of theconfidence interval are found in practice by solving pθ = α for θ.
In the procedure outlined above, one is still free to choose the test to be used; this corresponds
to the freedom in the Neyman construction as to which values of the data are included in theconfidence belt. One possibility is to use a test statistic based on the likelihood ratio,λ(θ) =

f (x; θ)f (x; θb ),(40.75)where θb is the value of the parameter which, out of all allowed values, maximizes f (x; θ). The
confidence belt is taken to contain the values of x that give the greatest values of λ(θ). This resultsin the intervals described in Ref. [42] by Feldman and Cousins.
If the model contains nuisance parameters ν, then these can be incorporated into the test (orthe p-values) used to determine the limit by profiling as discussed in Section 40.3.2.1. As mentioned
there, the strict frequentist approach is to regard the parameter of interest θ as excluded only if it
is rejected for all possible values of ν. The resulting interval for θ will then cover the true value
with a probability greater than or equal to the nominal confidence level for all points in ν-space.
bb
If the p-value is based on the profiled values of the nuisance parameters, i.e., with ν = ν(θ)
usedin Eq. (40.48), then the resulting interval for the parameter of interest will have the correct coverage
if the true values of ν are equal to the profiled values. Otherwise the coverage probability may be
too high or too low. This procedure has been called profile construction in particle physics [32] (seealso Ref. [26]).
31st May, 2024

40. Statistics

27

40.4.2.2 Gaussian distributed measurements
An important example of constructing a confidence interval is when the data consists of a single
random variable x that follows a Gaussian distribution; this is often the case when x represents anestimator for a parameter and one has a sufficiently large data sample. If there is more than one
parameter being estimated, the multivariate Gaussian is used. For the univariate case with known
σ, the probability that the measured value x will fall within ±δ of the true value µ is1−α= √

1
2πσ

Z µ+δ



δ
2σ

= erf

√

2

2

δ
σ

−1,

e−(x−µ) /2σ dx

µ−δ



 

= 2Φ

(40.76)

where erf is the Gaussian error function, which is rewritten in the final equality using Φ, theGaussian cumulative distribution. Fig. 40.4 shows a δ = 1.64σ confidence interval unshaded. The
choice δ = σ gives an interval called the standard error which has 1 − α = 68.27% if σ is known.
Values of α for other frequently used choices of δ are given in Table 40.1.

f (x; µ,σ)

1− α

α /2
−3

α /2
−2

−1

0
1
(x−µ) /σ2

3Figure 40.4: Illustration of a symmetric 90% confidence interval (unshaded) for a Gaussiandistributed measurement of a single quantity. Integrated probabilities, defined by α = 0.1, are asshown.
We can set a one-sided (upper or lower) limit by excluding above x + δ (or below x − δ). The
values of α for such limits are half the values in Table 40.1.
31st May, 2024

40. Statistics

28

Table 40.1: Area of the tails α outside ±δ from the mean of a Gaussian
distribution.
α0.3173
4.55 ×10−2
2.7 ×10−3
6.3×10−5
5.7×10−7
2.0×10−9

δ
1σ
2σ
3σ
4σ
5σ
6σ

α
0.2
0.1
0.05
0.01
0.001
10−4

δ
1.28σ
1.64σ
1.96σ
2.58σ
3.29σ
3.89σ

The relation (40.76) can be re-expressed using the cumulative distribution function for the χ2
distribution as
α = 1 − F (χ2 ; n) ,(40.77)for χ2 = (δ/σ)2 and n = 1 degree of freedom. This can be seen as the n = 1 curve in Fig. 40.1 orobtained using standard computer libraries. For multivariate measurements of, say, M parameter
estimates θb = (θb1 , . . . , θbM ), construction of the confidence region requires the full covariance matrix
Vij = cov[θbi , θbj ], which can be estimated as described in Sections 40.2.2 and 40.2.3. Under fairly
general conditions with the methods of maximum-likelihood or least-squares in the large sample
limit, the estimators will be distributed according to a multivariate Gaussian centered about thetrue (unknown) values θ, and furthermore, the likelihood function itself will take on a Gaussian
shape.
The standard error ellipse for the pair (θbi , θbj ) is shown in Fig. 40.5, corresponding to a contourb andχ2 = χ2min + 1 or ln L = ln Lmax − 1/2. The ellipse is centered about the estimated values θ,
the tangents to the ellipse give the standard deviations of the estimators, σi and σj . The angle ofthe major axis of the ellipse is given bytan 2φ =2ρij σi σj,σj2 − σi2(40.78)where ρij = cov[θbi , θbj ]/σi σj is the correlation coefficient.
The correlation coefficient can be visualized as the fraction of the distance σi from the ellipse’s
horizontal center-line at which the ellipse becomes tangent to vertical, i.e., at the distance ρij σibelow the center-line as shown. As ρij goes to +1 or −1, the ellipse thins to a diagonal line.
It could happen that one of the parameters, say, θj , is known from previous measurements toa precision much better than σj , so that the current measurement contributes almost nothing to
the knowledge of θj . However, the current measurement of θi and its dependence on θj may still
be important. In this case, instead of quoting both parameter estimates and their correlation, one
sometimes reports the value of θi , which minimizes χ2 at a fixed value of θj , such as the PDG bestvalue. This θi value lies along the dotted line between the points where the ellipse becomes tangent
to vertical, and has statistical error σinner as shown on the figure, where σinner = (1 − ρ2ij )1/2 σi .
Instead of the correlation ρij , one reports the dependency dθbi /dθj , which is the slope of the dottedline. This slope is related to the correlation coefficient by dθbi /dθj = ρij × σσji .
As in the single-variable case, because of the symmetry of the Gaussian function between θb one finds that contours of constant ln L or χ2 cover the true values with a certain, fixedand θ,
probability. That is, the confidence region is determined byln L(θ) ≥ ln Lmax − ∆ ln L ,31st May, 2024

(40.79)40. Statistics

29

θi

σi

σ inner

^

θi

ρijσi

σi
σj

φ

σj
θ^j

θj

Figure 40.5: Standard error ellipse for the estimators θbi and θbj . In the case shown the correlationis negative.
or where a χ2 has been defined for use with the method of least-squares,
χ2 (θ) ≤ χ2min + ∆χ2 .

(40.80)

Values of ∆χ2 or 2∆ ln L are given in Table 40.2 for several values of the coverage probability
1 − α and number of fitted parameters M . For Gaussian distributed data, these are related by−1
is the chi-square quantile (inverse of the cumulative∆χ2 = 2∆ ln L = Fχ−1
2 (1 − α), where Fχ2
M

M

distribution) for M degrees of freedom.
Table 40.2: Values of ∆χ2 or 2∆ ln L corresponding to a coverage probability 1 − α in the large data sample limit, for joint estimation of M
parameters.
(1 − α) (%)
68.27
90.
95.
95.45
99.
99.73

M =1
1.00
2.71
3.84
4.00
6.63
9.00

M =2
2.30
4.61
5.99
6.18
9.21
11.83

M =3
3.53
6.25
7.82
8.03
11.34
14.16

For non-Gaussian data samples, the probability for the regions determined by Equations (40.79)
or (40.80) to cover the true value of θ becomes independent of θ only in the large-sample limit. So
for a finite data sample these are not exact confidence regions according to our previous definition.
Nevertheless, they can still have a coverage probability only weakly dependent on the true parameter, and approximately as given in Table 40.2. In any case, the coverage probability of the intervals
or regions obtained according to this procedure can in principle be determined as a function of thetrue parameter(s), for example, using a Monte Carlo calculation.
One of the practical advantages of intervals that can be constructed from the log-likelihood
function or χ2 is that it is relatively simple to produce the interval for the combination of several
31st May, 2024

40. Statistics

30

experiments. If N independent measurements result in log-likelihood functions ln Li (θ), then the
combined log-likelihood function is simply the sum,ln L(θ) =N
X

ln Li (θ) .

(40.81)i=1This can then be used to determine an approximate confidence interval or region with Eq. (40.79),
just as with a single experiment.
40.4.2.3 Poisson or binomial data
Another important class of measurements consists of counting a certain number of events, n.
In this section, we will assume these are all events of the desired type, i.e., there is no background.
If n represents the number of events produced in a reaction with cross section σ and selection
efficiency ε in a fixed integrated luminosity L, then it follows a Poisson distribution with mean
µ = σεL. If, on the other hand, one has selected a larger sample of N events and found n ofthem to have a particular property, then n follows a binomial distribution where the parameter p
gives the probability for the event to possess the property in question. This is appropriate, e.g.,
for estimates of branching ratios or selection efficiencies based on a given total number of events.
For the case of Poisson distributed n, limits on the mean value µ can be found from the Neyman
procedure as discussed in Section 40.4.2.1 with n used directly as the statistic x . The upper andlower limits are found to be1µlo = Fχ−1
(40.82a)2 (αlo ; 2n) ,21µup = Fχ−1
(40.82b)2 (1 − αup ; 2(n + 1)) ,2
where confidence levels of 1 − αlo and 1 − αup refer separately to the corresponding intervals µ ≥ µlo2
and µ ≤ µup , and Fχ−1
2 is the quantile of the χ distribution (inverse of the cumulative distribution).
For central confidence intervals at confidence level 1 − α, set αlo = αup = α/2.
Table 40.3: Lower and upper (one-sided) limits for the mean µ of a
Poisson variable given n observed events in the absence of background, forconfidence levels of 90% and 95%.

n
0
1
2
3
4
5
6
7
8
9
10

1 − α =90%
µlo
µup
–
2.30
0.105 3.89
0.532 5.32
1.10
6.68
1.74
7.99
2.43
9.27
3.15 10.53
3.89 11.77
4.66 12.99
5.43 14.21
6.22 15.41

1 − α =95%
µlo
µup
–
3.00
0.051 4.74
0.355 6.30
0.818 7.75
1.37
9.15
1.97 10.51
2.61 11.84
3.29 13.15
3.98 14.43
4.70 15.71
5.43 16.96

It happens that the upper limit from Eq. (40.82b) coincides numerically with the Bayesian
upper limit for a Poisson parameter, using a uniform prior p.d.f. for µ. Values for confidence levels
31st May, 2024

40. Statistics

31

of 90% and 95% are shown in Table 40.3. For the case of binomially distributed n successes out ofN trials with probability of success p, the upper and lower limits on p are found to beplo =nFF−1 [αlo ; 2n, 2(N − n + 1)],
N − n + 1 + nFF−1 [αlo ; 2n, 2(N − n + 1)](40.83a)pup =(n + 1)FF−1 [1 − αup ; 2(n + 1), 2(N − n)].
(N − n) + (n + 1)FF−1 [1 − αup ; 2(n + 1), 2(N − n)](40.83b)

Here FF−1 is the quantile of the F distribution (also called the Fisher–Snedecor distribution; see
Ref. [4]).
40.4.2.4 Parameter exclusion in cases of low sensitivity
An important example of a statistical test arises in the search for a new signal process. Suppose
the parameter µ is defined such that it is proportional to the signal cross section. A statistical test
may be carried out for hypothesized values of µ, which may be done by computing a p-value, pµ ,for all µ. Those values not rejected in a test of size α, i.e., for which one does not find pµ ≤ α,constitute a confidence interval with confidence level 1 − α.
In general one will find that for some regions in the parameter space of the signal model, thepredictions for data are almost indistinguishable from those of the background-only model. This
corresponds to the case where µ is very small, as would occur, e.g., in a search for a new particle
with a mass so high that its production rate in a given experiment is negligible. That is, one has
essentially no experimental sensitivity to such a model.
One would prefer that if the sensitivity to a model (or a point in a model’s parameter space) isvery low, then it should not be excluded. Even if the outcomes predicted with or without signal are
identical, however, the probability to reject the signal model will equal α, the type-I error rate. As
one often takes α to be 5%, this would mean that in a large number of searches covering a broad
range of a signal model’s parameter space, there would inevitably be excluded regions in which theexperimental sensitivity is very small, and thus one may question whether it is justified to regard
such parameter values as disfavored.
Exclusion of models to which one has little or no sensitivity occurs, for example, if the data
fluctuate very low relative to the expectation of the background-only hypothesis. In this case theresulting upper limit on µ may be anomalously low. As a means of controlling this effect one
often determines the mean or median limit under assumption of the background-only hypothesis,
as discussed in Sec. 40.5.
One way to mitigate the problem of excluding models to which one is not sensitive is the CLs
method, where the measure used to test a parameter is increased for decreasing sensitivity [43, 44].
The procedure is based on a statistic called CLs , which is defined as
CLs =

pµ
,
1 − pb(40.84)

where pb is the p-value of the background-only hypothesis. In the usual formulation of the method,both pµ and pb are defined using a single test statistic, and the definition of CLs above assumes
this statistic is continuous; more details can be found in Refs. [43, 44].
A point in a model’s parameter space is regarded as excluded if one finds CLs ≤ α. As thedenominator in Eq. (40.84) is always less than or equal to unity, the exclusion criterion based on
CLs is more stringent than the usual requirement pµ ≤ α. In this sense the CLs procedure isconservative, and the coverage probability of the corresponding intervals will exceed the nominal
confidence level 1 − α. If the experimental sensitivity to a given value of µ is very low, then one
31st May, 2024

32

40. Statistics

finds that as pµ decreases, so does the denominator 1 − pb , and thus the condition CLs ≤ α iseffectively prevented from being satisfied. In this way the exclusion of parameters in the case of
low sensitivity is suppressed.
The CLs procedure has the attractive feature that the resulting intervals coincide with thoseobtained from the Bayesian method in two important cases: the mean value of a Poisson or Gaussian
distributed measurement with a constant prior. The CLs intervals overcover for all values of theparameter µ, however, by an amount that depends on µ.
The problem of excluding parameter values to which one has little sensitivity is particularly
acute when one wants to set a one-sided limit, e.g., an upper limit on a cross section. Here one
tests a value of a rate parameter µ against the alternative of a lower rate, and therefore the criticalregion of the test is taken to correspond to data outcomes with a low event yield. If the number of
events found in the search region fluctuates low enough, however, it can happen that all physicallymeaningful signal parameter values, including those to which one has very little sensitivity, arerejected by the test.
Another solution to this problem, therefore, is to replace the one-sided test by one based on
the likelihood ratio, where the critical region is not restricted to low rates. This is the approach
followed in the Feldman-Cousins procedure described in Section 40.4.2.1. The critical region for
the test of a given value of µ contains data values characteristic of both higher and lower rates.
As a result, for a given observed rate one can in general obtain a two-sided interval. If, however,the parameter estimate µ̂ is sufficiently close to the lower limit of zero, then only high values ofµ are rejected, and the lower edge of the confidence interval is at zero. Note, however, that thecoverage property of 1 − α pertains to the entire interval, not to the probability for the upper edge
µup to be greater than the true value µ. For parameter estimates increasingly far away from theboundary, i.e., for increasing signal significance, the point µ = 0 is excluded and the interval hasnonzero upper and lower edges.
An additional difficulty arises when a parameter estimate is not significantly far away from theboundary, in which case it is natural to report a one-sided confidence interval (often an upper limit).
It is straightforward to force the Neyman prescription to produce only an upper limit by setting
x2 = ∞ in Eq. (40.73). Then x1 is uniquely determined and the upper limit can be obtained. If,
however, the data come out such that the parameter estimate is not so close to the boundary, one
might wish to report a central confidence interval (i.e., an interval based on a two-sided test withequal upper and lower tail areas). As pointed out by Feldman and Cousins [42], if the decision toreport an upper limit or two-sided interval is made by looking at the data (“flip-flopping”), then ingeneral there will be parameter values for which the resulting intervals have a coverage probability
less than 1 − α. With the confidence intervals suggested in Ref. [42], the prescription determines
whether the interval is one- or two-sided in a way which preserves the coverage probability (andare thus said to be unified).
The intervals according to this method for the mean of Poisson variable in the absence of
background are given in Table 40.4. (Note that α in Ref. [42] is defined following Neyman [41] as
the coverage probability; this is opposite the modern convention used here in which the coverage
probability is 1 − α.) The values of 1 − α given here refer to the coverage of the true parameter by
the whole interval [µ1 , µ2 ]. In Table 40.3 for the one-sided upper limit, however, 1 − α refers to theprobability to have µup ≥ µ (or µlo ≤ µ for lower limits).
A potential difficulty with unified intervals arises if, for example, one constructs such an interval
for a Poisson parameter s of some yet to be discovered signal process with, say, 1 − α = 0.9. If
the true signal parameter is zero, or in any case much less than the expected background, one willusually obtain a one-sided upper limit on s. In a certain fraction of the experiments, however, atwo-sided interval for s will result. Since, however, one typically chooses 1 − α to be only 0.9 or 0.95
31st May, 2024

40. Statistics

33

Table 40.4: Unified confidence intervals [µ1 , µ2 ] for a the mean of a
Poisson variable given n observed events in the absence of background, forconfidence levels of 90% and 95%.
1 − α =90%n µ1
µ2
0 0.002.44
1 0.114.36
2 0.535.91
3 1.107.42
4 1.478.60
5 1.849.99
6 2.21 11.47
7 3.56 12.53
8 3.96 13.99
9 4.36 15.30
10 5.50 16.50

1 − α =95%
µ1
µ2
0.00
3.09
0.05
5.14
0.36
6.72
0.82
8.25
1.37
9.76
1.84 11.26
2.21 12.75
2.58 13.81
2.94 15.29
4.36 16.77
4.75 17.82

when setting limits, the value s = 0 may be found below the lower edge of the interval before theexistence of the effect is well established. It must then be communicated carefully that in excludings = 0 at, say, 90% or 95% confidence level from the interval, one is not necessarily claiming to havediscovered the effect, for which one would usually require a higher level of significance (e.g., 5 σ).
Another possibility is to construct a Bayesian interval as described in Section 40.4.1. The
presence of the boundary can be incorporated simply by setting the prior density to zero in theunphysical region. More specifically, the prior may be chosen using formal rules such as the reference
prior or Jeffreys prior mentioned in Sec. 40.2.6.
In particle physics a widely used prior for the mean µ of a Poisson distributed measurement
has been the uniform distribution for µ ≥ 0. This prior does not follow from any fundamental
rule nor can it be regarded as reflecting a reasonable degree of belief, since the prior probability
for µ to lie between any two finite values is zero. The procedure above can be more appropriately
regarded as a way for obtaining intervals with frequentist properties that can be investigated. The
resulting upper limits have a coverage probability that depends on the true value of the Poissonparameter, and is nowhere smaller than the stated probability content. Lower limits and two-sided
intervals for the Poisson mean based on uniform priors undercover, however, for some values of theparameter, although to an extent that in practical cases may not be too severe [2, 40].
In any case, it is important to always report sufficient information so that the result can be
combined with other measurements. Often this means giving an unbiased estimator and its standard
deviation, even if the estimated value is in the unphysical region.
It can also be useful with a frequentist interval to calculate its subjective probability content
using the posterior p.d.f. based on one or several reasonable guesses for the prior p.d.f. If it turns
out to be significantly less than the stated confidence level, this warns that it would be particularly
misleading to draw conclusions about the parameter’s value from the interval alone.

40.5 Experimental sensitivity
In this section we describe methods for characterizing the sensitivity of a search for a newphysics signal. As discussed in Sec. 40.3, an experimental analysis can often be formulated as a testof hypothetical model parameters. Therefore we may quantify the sensitivity by giving the results
31st May, 2024

34

40. Statistics

that we expect from such a test under specific assumptions about the signal process.
Here to be concrete we will consider a parameter µ proportional to the rate of a signal process,although the concepts described in this section may be easily generalized to other parameters. One
may wish to establish discovery of the signal process by testing and rejecting the hypothesis thatµ = 0, and in addition one often wants to test nonzero values of µ to construct a confidence interval(e.g., limits) as described in Sec. 40.4. In the frequentist framework, the result of each tested valueof µ is the p-value pµ or equivalently the significance Zµ = Φ−1 (1 − pµ ), where as usual Φ is thestandard Gaussian cumulative distribution and its inverse Φ−1 is the standard Gaussian quantile.
Prior to carrying out the experiment, one generally wants to quantify what significance Zµ is
expected under given assumptions for the presence or absence of the signal process. Specifically, forthe significance of a test of µ = 0 (the discovery significance) one usually quotes the Z0 one would
expect if the signal is present at a given nominal rate, which we can define in general to correspondto µ = 1. For limits, one often gives the expected limit under assumption of the background-only(µ = 0) model. These quantities are used to optimize the analysis and to quantify the experimental
sensitivity, that is, to characterize how likely it is to make a discovery if the signal is present, andto say what values of µ one may be able to exclude if the signal is in fact absent.
First we clarify the notion of expected significance. Because the significance Zµ is a function of
the data, it is itself a random quantity characterized by a certain sampling distribution. This distribution depends on the assumed value of µ, which is not necessarily the same as the hypothesizedvalue of µ being tested. We may therefore consider the distribution f (Zµ |µ0 ), i.e., the distributionof Zµ that would be obtained by considering data samples generated under assumption of µ0 . In asimilar way one can talk about the sampling distribution of an upper limit for µ, f (µup |µ0 ).
One can identify the expected significance or limit with either the mean or median of these distributions, but the median may be preferred since it is invariant under monotonic transformations.
For example, the monotonic relation between p-value and significance, p = 1 − Φ(Z), then gives
med[pµ |µ0 ] = 1 − Φ(med[Zµ |µ0 ]), whereas the corresponding relation does not hold in general forthe mean.
In some cases one may be able to write down approximate formulae for the distributions ofZµ and for limits, but more generally they must be determined from Monte Carlo calculations. In
many cases of interest, the significance Zµ and the limits on µ will have approximate Gaussian
distributions.
As an example, consider a Poisson counting experiment, where the result consists of an observed
number n of events, modeled as a Poisson distributed variable with a mean of µs + b. Here s andb, the expected numbers of events from signal and background processes, are taken to be known.
If we are interested in discovering the signal process we test and try to reject the hypothesis µ = 0.
To characterize the experimental sensitivity, we want to give the discovery significance expectedunder the assumption of µ = 1.
In the limit where its mean value is large, the Poisson variable n can be approximatedas an√almost continuous Gaussian variable with mean µs + b and standard deviation σ = µs + b. In
the usual case where a physical signal model corresponds to µ > 0, the p-value of µ = 0 is theprobability to find n greater than or equal to the value observed,n−bp0 = Φ √b

,(40.85)√and the corresponding significance is Z0 = Φ−1 (1 − p0 ) = (n − b)/ b. The median (here equal tothe mean) of n assuming µ = 1 is s + b, and therefore the median discovery significance iss
med[Z0 |µ = 1] = √ .
(40.86)b
31st May, 2024

40. Statistics

35√
The figure of merit “s/ b” has been widely used in particle physics as a measure of expecteddiscovery significance. A better approximation for the Poisson counting experiment, however, may
be obtained by testing µ = 0 using the likelihood ratio (40.49) λ(0) = L(0)/L(µ̂), whereL(µ) =(µs + b)n −(µs+b)en!

(40.87)is the likelihood function, and µ̂ = (n − b)/s is the estimator of µ. In this example there are nonuisance parameters, as s and b are taken to be known. For the case where the relevant signal models
correspond to positive µ, one may test the µ = 0 hypothesis with the statistic q0 = −2 ln λ(0) whenµ̂ > 0, i.e., an excess is observed, and q0 = 0 otherwise. One can show (see, e.g., Ref. [25]) that in√
the large-sample limit, the discovery significance is then Z0 = q0 , for which one finds
s 

Z0 =

n
2 n ln + b − nb(40.88)for n > b and Z0 = 0 otherwise. To approximate the expected discovery significance assumingµ = 1, one may simply replace n with the expected value E[n|µ = 1] = s + b (the so-called “Asimov
data set”), giving
s 



s
−s .
(40.89)med[Z0 |µ = 1] = 2 (s + b) ln 1 +bThis has been shown in Ref. [25] to provide a good approximation to the median discovery significance for valuesof s of several and for b well below unity. The right-hand side of Eq. (40.89)√
reduces to s/ b in the limit s  b.
Beyond the simple Poisson counting experiment, in general one may test values of a parameterµ with more complicated functions of the measured data to obtain a p-value pµ , and from thisone can quote the equivalent significance Zµ or find, e.g., an upper limit µup . In this case as well
one may quantify the experimental sensitivity by giving the significance Zµ expected if the data
are generated with a different value of the parameter µ0 . In some problems, finding the sampling
distribution of the significance or limits may be possible using large-sample formulae as described,
e.g., in Ref. [25]. In other cases a Monte Carlo study may be needed. Using whatever methodof calculation is most appropriate, one usually quotes the expected (mean or, preferably, median)discovery significance or exclusion limit as the primary measure of experimental sensitivity.
Even if the true signal is present at its nominal rate, the actual discovery significance Z0
obtained from the real data is subject to statistical fluctuations and will not in general be equal toits expected value. In an analogous way, the observed limit will differ from the expected limit even
if the signal is absent. Upon observing such a difference one would like to know how large this iscompared to expected statistical fluctuations. Therefore, in addition to the observed significanceand limits it is useful to communicate not only their expected values but also a measure of thewidth of their distributions.
As the distributions of significance and limits are often well approximated by a Gaussian, one
may indicate the intervals corresponding to plus-or-minus one and/or two standard deviations.
If the distributions are significantly non-Gaussian, one may use instead the quantiles that givethe same probability content, i.e., [0.1587, 0.8413] for ±1σ, [0.02275, 0.97725] for ±2σ. An upper
limit found significantly below the background-only expectation may indicate a strong downward
fluctuation of the data, or perhaps as well an incorrect estimate of the background rate.
The procedures described above pertain to frequentist hypothesis tests and limits. Bayesian
limits, just like those found from a frequentist procedure, are functions of the data and one may
31st May, 2024

40. Statistics

36

therefore find, usually with approximations or Monte Carlo studies, their sampling distribution andcorresponding mean (or, preferably, median) and standard deviation.
When trying to establish discovery of a signal process, the Bayesian approach may employ
a Bayes factor as described in Sec. 40.3.4. In the case of the Poisson counting experiment withthe likelihood from Eq. (40.87), the log of the Bayes factor that compares µ = 1 to µ = 0 isln B10 = ln(L(1)/L(0)) = n ln(1 + s/b) − s. That is, the expectation value, assuming µ = 1, ofln B10 for this problem isE[ln B10 |µ = 1] = (s + b) ln 1 +

s
b



−s.

(40.90)pComparing this to Eq. (40.89), one finds med[Z0 |1] = 2E[ln B10 |1]. Thus for this particular
problem the frequentist median discovery significance can be related to the corresponding Bayes
factor in a simple way.
In some analyses, the goal may not be to establish discovery of a signal process but rather
to measure, as accurately as possible, the signal rate. If we consider again the Poisson counting
experiment described by the likelihood function of Eq. (40.87), the maximum-likelihood estimatorµ̂ = (n − b)/s has a variance, assuming µ = 1, of
n−b
1
s+b
V [µ̂] = V= 2 V [n] = 2 ,(40.91)sss√√so that the standard deviation of µ̂ is σµ̂ = s + b/s. One may therefore use s/ s + b as a figure of
merit to be √maximized in order to obtain the best measurement accuracy of a rate parameter. The
quantity s/ s + b is also the expected significance with which one rejects s assuming the signal isabsent, and thus can be used to optimize the expected upper limit on s.


References[1] B. Efron, Am. Stat. 40, 11 (1986).
[2] R. D. Cousins, Am. J. Phys. 63, 398 (1995).
[3] A. Stuart, J.K. Ord, and S. Arnold, Kendall’s Advanced Theory of Statistics, Vol. 2A: ClassicalInference and the Linear Model, 6th ed., Oxford Univ. Press (1999), and earlier editions by
Kendall and Stuart. The likelihood-ratio ordering principle is described at the beginning of
Ch. 23. Chapter 26 compares different schools of statistical inference.
[4] F. James, Statistical methods in experimental physics (2006), ISBN 9789812567956.
[5] P. R. Rider, Journal of the American Statistical Association 55, 289, 148 (1960).
[6] L. Lyons, Statistics for Nuclear and Particle Physicists (1986), ISBN 9780521379342.
[7] R. J. Barlow, Nucl. Instrum. Meth. A297, 496 (1990).
[8] G. Cowan, Statistical data analysis (1998), ISBN 9780198501565.
[9] S. Baker and R. D. Cousins, Nucl. Instrum. Meth. 221, 437 (1984).
[10] S. S. Wilks, Annals Math. Statist. 9, 1, 60 (1938).
[11] O. Behnke et al., editors, Data analysis in high energy physics, Wiley-VCH, Weinheim,
Germany (2013), ISBN 9783527410583, 9783527653447, 9783527653430, URL http://www.
wiley-vch.de/publish/dt/books/ISBN3-527-41058-9.
[12] S. Schmitt, EPJ Web of Conferences 137, 11008 (2017), ISSN 2100-014X, URL http://dx.
doi.org/10.1051/epjconf/201713711008.
[13] L. Brenner et al., Int. J. Mod. Phys. A 35, 24, 2050145 (2020), [arXiv:1910.14654].
31st May, 2024

37

40. Statistics[14] A. O’Hagan and J.J. Forster, Bayesian Inference, (2nd edition, volume 2B of Kendall’s Advanced Theory of Statistics, Arnold, London, 2004).
[15] D. Sivia and J. Skilling, Data Analysis: A Bayesian Tutorial, (Oxford University Press, 2006).
[16] P.C. Gregory, Bayesian Logical Data Analysis for the Physical Sciences, (Cambridge University
Press, 2005).
[17] J.M. Bernardo and A.F.M. Smith, Bayesian Theory, (Wiley, 2000).
[18] Robert E. Kass and Larry Wasserman, J. Am. Stat. Assoc. 91, 1343 (1996).
[19] J.M. Bernardo, J. R. Statist. Soc. B41, 113 (1979); J.M. Bernardo and J.O. Berger, J. Am.
Stat. Assoc. 84, 200 (1989). See also J.M. Bernardo, Reference Analysis, in Handbook of Statistics, 25 (D.K. Dey and C.R. Rao, eds.), 17-90, Elsevier (2005) and references therein.
[20] L. Demortier, S. Jain and H. B. Prosper, Phys. Rev. D82, 034002 (2010), [arXiv:1002.1111].
[21] K. Cranmer, J. Pavez and G. Louppe (2015), [arXiv:1506.02169].
[22] N. Reid, Likelihood Inference in the Presence of Nuisance Parameters, Proceedings of PHYSTAT2003, L. Lyons, R. Mount, and R. Reitmeyer, eds., eConf C030908, Stanford, 2003.
[23] W. A. Rolke, A. M. Lopez and J. Conrad, Nucl. Instrum. Meth. A551, 493 (2005),[arXiv:physics/0403059].
[24] Links to the Proceedings of the PHYSTAT conference series (Durham 2002, Stanford 2003,
Oxford 2005, and Geneva 2007, 2011) can be found at https://espace.cern.ch/phystat.
[25] G. Cowan et al., Eur. Phys. J. C71, 1554 (2011), [arXiv:1007.1727]; G. Cowan et al., Eur.
Phys. J. C73, 1434 (2013).
[26] L. Demortier, P-Values and Nuisance Parameters, Proceedings of PHYSTAT 2007, CERN2008-001, p. 23.
[27] O. Barndorff-Nielsen, Biometrika 67, 2, 293 (1980), ISSN 00063444, URL http://www.jstor.
org/stable/2335474.
[28] M. S. Bartlett, Proceedings of the Royal Society of London. Series A - Mathematical andPhysical Sciences 160, 901, 268 (1937), URL https://royalsocietypublishing.org/doi/
abs/10.1098/rspa.1937.0109.
[29] D. N. Lawley, Biometrika 43, 3-4, 295 (1956), ISSN 0006-3444, URL https://doi.org/10.
1093/biomet/43.3-4.295.
[30] A. R. Brazzale, A. C. Davison and N. Reid, Applied Asymptotics: Case Studies in Small-Sample
Statistics, Cambridge University Press, Cambridge (2007), URL http://infoscience.epfl.
ch/record/104219.
[31] G. Cordeiro and F. Cribari-Neto, An Introduction to Bartlett Correction and Bias Reduction,
SpringerBriefs in Statistics, Springer Berlin Heidelberg (2014), ISBN 9783642552557, URL
https://books.google.co.uk/books?id=MJyKAwAAQBAJ.
[32] K. Cranmer, in “Statistical Problems in Particle Physics, Astrophysics and Cosmology(PHYSTAT 05): Proceedings, Oxford, UK, September 12-15, 2005,” 112–123 (2005),[arXiv:physics/0511028].
[33] R. D. Cousins, J. T. Linnemann and J. Tucker, Nucl. Instrum. Meth. A 595, 2, 480 (2008),
URL https://doi.org/10.1016%2Fj.nima.2008.07.086.
[34] E. Gross and O. Vitells, Eur. Phys. J. C70, 525 (2010), [arXiv:1005.1891].
[35] R. B. Davies, Biometrika 74, 33 (1987).
31st May, 2024

38

40. Statistics[36] R. D. Cousins, Lectures on Statistics in Theory: Prelude to Statistics in Practice (2023),[arXiv:1807.05996].
[37] J. Skilling, Nested Sampling, AIP Conference Proceedings, 735, 395–405 (2004).
[38] F. Feroz, M. P. Hobson and M. Bridges, Mon. Not. Roy. Astron. Soc. 398, 1601 (2009),[arXiv:0809.3437].
[39] R. E. Kass and A. E. Raftery, J. Am. Statist. Assoc. 90, 430, 773 (1995).
[40] B. P. Roe and M. B. Woodroofe, Phys. Rev. D63, 013009 (2001), [hep-ex/0007048].
[41] J. Neyman, Phil. Trans. Roy. Soc. Lond. A236, 767, 333 (1937); Reprinted in A Selection of
Early Statistical Papers on J. Neyman, (University of California Press, Berkeley, 1967).
[42] G. J. Feldman and R. D. Cousins, Phys. Rev. D57, 3873 (1998), [arXiv:physics/9711021]; This
paper does not specify what to do if the ordering principle gives equal rank to some values ofx. Eq. 21.6 of Ref. [3] gives the rule: all such points are included in the acceptance region (thedomain D(α)). Some authors have assumed the contrary, and shown that one can then obtainnull intervals.
[43] A. L. Read, in “Workshop on confidence limits, CERN, Geneva, Switzerland, 17-18 Jan 2000:
Proceedings,” 81–101 (2000), URL http://weblib.cern.ch/record/451614.
[44] T. Junk, Nucl. Instrum. Meth. A434, 435 (1999), [hep-ex/9902006].

31st May, 2024