8.2 Frequentist statistics at the LHC

8.2.1 The likelihood function and test statistics

The data model

Let us take the simplest possible case of a (one bin) counting experiment, where in our “signal region” we expect $s$ signal events and $b$ background events. The probability to observe $n$ events in our signal region is distributed as a Poisson with mean $s+b$:

Since we have only one observation but two free parameters, this experiment is underconstrained. So, let’s also add a “control region” where we expect no signal and $b$ background events. The probability of observing $m$ events in our control region is therefore:

Combining the two, the joint probability distribution for $n$ and $m$ is:

This is also called the model for the data and is plotted for sample $s,b$ values in Figure 8.1.

The likelihood function

In the frequentist philosophy, however, all our parameters $n$, $m$ etc. are simply fixed values of nature and, hence, don’t have a probability distribution. Instead, we work with the likelihood function, which is a function only of our parameters of interest (POIs), $s$ in our example, and “nuisance parameters” ($b$), given fixed values for $n$ and $m$:

Importantly, this is not a probability distribution on $s$ and $b$! To derive that, we would have to use Bayes’ rule to go from $P(n,m;s,b)\to P(s,b;n,m)$; however, such probability distributions don’t make sense in our frequentist world view, so we’re stuck with this likelihood formulation. Often, it’s more convenient to consider the negative log-likelihood:

The profile likelihood ratio

Fundamentally, the goal of any experiment is to test the compatibility of the observed data ($n,m$ here) with a certain hypothesis $H$. We do this by mapping the data to a “test statistic” $t$, which is just a number, and comparing it against its distribution under $H$, $P(t|H)$. Our problem, thus, boils down to 1) choosing the most effective $t$ for testing $H$, and 2) obtaining $P(t|H)$.

In the case of testing a particular signal strength, we use the “profile likelihood ratio”:

where $ŝ,\widehat{b}$ are the maximum-likelihood estimates (MLEs) for $s$ and $b$, given the observations $n,m$, and $\widehat{\widehat{b}}(s)$ is the MLE for $b$ given $n,m$, and $s$. The MLE for a parameter is simply the value of it for which the likelihood is maximized, and will be discussed in the next section. The numerator of $\lambda (s)$ can be thought of as a way to “marginalize” over the nuisance parameters by simply values that maximize the likelihood for any given $s$, while the denominator is effectively a normalization factor, such that $\lambda (s)\le 1$.

Again, it’s often more convenient to use the (negative) logarithm:

Note that $\mathrm{Max}[\lambda (s)]=1\Rightarrow \mathrm{Min}[t_s]=0$. $\lambda (s)$ and $t_s$ are plotted for sample $n,m$ values with $n-m=10$ in Figure 8.2. The maximum (minimum) of the profile likelihood ratio ($t_s$) is at $s=n-m=10$, as we expect; however, as the ratio between $n$ and $m$ decreases — i.e., the experiment becomes more noisy — the distributions broaden, representing the reduction in sensitivity, or the higher uncertainty on the true value of $s$.

Note that the likelihood ratio and $t_s$ are also broadened due to the nuisance parameter; i.e., because we are missing information about $b$. This can be demonstrated by plotting them with $b=m$, emulating perfect information of $b$ (Figure 8.3), and indeed, we see the functions are narrower than in Figure 8.2. More generally, increasing (decreasing) the uncertainties on the nuisance parameters will broaden (narrow) the test statistic distribution. This is which is why experimentally we want to constrain them through auxiliary measurements as much as possible.

Maximum-likelihood estimates

MLEs for $s$ and $b$ can be found for this example by setting the derivative of the negative log-likelihood to 0 (more generally, this would require numerical minimization):

Solving simultaneously yields, as you might expect:

Or just for $\widehat{\widehat{b}}(s)$ from Eq. 8.2.8:

Plugging this back in, we can get $\lambda (s)$ and $t_s$ for any given $s$.

Alternative test statistic

So far, our construction allows for $s<0$; however, physically the number of signal events can’t be negative. Rather than incorporating this constraint in the model, it’s more convenient to impose this in the test statistic, by defining:

and

The difference between the nominal and alternative test statistics is highlighted in Figure 8.4. For $n<m$, the $\tilde{\lambda }(s)=1$ and ${\tilde{t}}_s=0$ values are at $s=0$, since physically that is what fits best with our data (even though the math says otherwise).

Next, we want to translate this to a probability distribution of ${\tilde{t}}_s$ under a particular signal hypothesis ($H_s$) (i.e., an assumed value of $s$): $p({\tilde{t}}_s|H_s)$, or just $p({\tilde{t}}_s|s)$ for simplicity.

8.2.2 Hypothesis testing

The goal of any experiment is to test whether our data support or exclude a particular hypothesis $H$, and quantify the (dis)agreement. For example, to what degree did our search for the Higgs boson agree or disagree with the standard model hypothesis?

We have already discussed the process of mapping data to a scalar test statistic $t$ that we can use to test $H$. However, we need to know the probability distribution of $t$ under $H$ to quantify the (in)consistency of the observed data with $H$ and decide whether or not to exclude $H$.

We must also recognize that there’s always a chance that we will exclude $H$ even if it’s true (called a Type I error, or a false positive), or not exclude $H$ when it’s false (Type II error, or false negative). The probability of each is referred to as $\alpha$ and $\beta$, respectively. This is summarized handily in Table 8.1.

Before the test, we should decide on a probability of making a Type I error, $\alpha$, that we are comfortable with, called a “significance level”. Typical values are 5% and 1%, although if we’re claiming something crazy like a new particle, we better be very sure this isn’t a false positive; hence, we set a much lower significance level for these tests of $3\times 10^{-7}$. (The significance of this value will be explained in Section 8.2.2. below.)

Deriving $p({\tilde{t}}_s|s)$

We can approximate $p({\tilde{t}}_s|s)$ by generating several pseudo- or “toy” datasets assuming $s$ expected signal events. In this case, this means sampling possible values for $n,m$ from our probability model. We will continue with our simple counting experiment (Section 8.2.1.), for which such toy datasets are generated and then used to create histograms for $p({\tilde{t}}_s|s)$ in Figure 8.5. Note that one complication in generating these toys is that the $n,m$ distributions from which we want to sample (Eq. 8.2.3) also depend on the nuisance parameter $b$. However, we see from the figure that this does not matter as much as we might expect.

We make two important observations:

1.

$p({\tilde{t}}_s|s)$ does not depend on nuisance parameters as long as we have sufficiently large statistics (in this case, when $b$ is sufficiently large). This is a key reason for basing our test statistic on the profile likelihood.

2.

In fact, $p({\tilde{t}}_s|s)$ doesn’t even depend on the POI $s$! (Again, as long as $s$ is large.)

Reference [301] shows that, asymptotically, this distribution follows a ${\chi }^2$ distribution with degrees of freedom equal to the number of POIs, as illustrated in Figure 8.6.¹ We can see that the asymptotic form looks accurate even for $s,b$ as low as $\sim 5$. Note that for cases where we can’t use the asymptotic form, Ref. [58] recommends using $b=\widehat{\widehat{b}}(s)$ when generating toys, so that we (approximately) maximize the agreement with the hypothesis.

$p$-values and significance

Now that we know the distribution of the test statistic $p({\tilde{t}}_s|H_s)\equiv p({\tilde{t}}_s|s)$, we can finally test $H_s$ with our experiment. We just need to calculate the “observed” test statistic ${\tilde{t}}_s^{\mathrm{obs}}$ from our observations, and compare it to the $p({\tilde{t}}_s|s)$.

Ultimately, we care about, given $p({\tilde{t}}_s|s)$, the probability of obtaining ${\tilde{t}}_s^{\mathrm{obs}}$ or a value more inconsistent with $H_s$; i.e., the green shaded region above. This is referred to as the $p$-value of the observation:

which is $0.30$ for this example, where

is the cumulative distribution function (CDF) of ${\tilde{t}}_s$. We reject the hypothesis if this $p$-value is less than our chosen significance level $\alpha$; the idea being that if $H_s$ were true and we repeated this measurement many times, then the probability of a false-positive ($p$-value $\le \alpha$) is exactly $\alpha$, as we intended.

The $p$-value is typically converted into a significance ($Z$), which is the corresponding number of standard deviations away from the mean in a Gaussian distribution:

where $\mathrm{\Phi }$ is the CDF of the standard Gaussian. This is more easily illustrated in Figure 8.8, where $\phi$ is the standard Gaussian distribution:

The significance in Example 8.2.1 is, therefore, ${\mathrm{\Phi }}^{-1}(1-0.30)=0.53$. We sometimes say that our measurement is (in)consistent or (in)compatible with $H$ at the $0.53\sigma$ level, or within $1\sigma$, etc.

Signal discovery

So far, we have been testing the signal hypothesis, but usually when searching for a particle, we instead test the “background-only” hypothesis $H_0$ and decide whether or not to reject it. This means we want ${\tilde{t}}_0^{\mathrm{obs}}$ and $p({\tilde{t}}_0|0)$ (Figure 8.9).²

We could say for this experiment, therefore, that we exclude the background-only hypothesis at the “3 sigma” level. However, for an actual search for a new particle at the LHC, this is insufficient to claim a discovery, as the probability of a false positive at $3\sigma$,  1/1000, is too high. The standard is instead set at $5\sigma$ for discovering new signals, corresponding to the $3\times 10^{-7}$ significance level quoted earlier, as we really don’t want to be making a mistake if we’re claiming to have discovered a new particle! $3\sigma$, $4\sigma$, and $5\sigma$ are commonly referred to as “evidence”, “observation”, and “discovery”, respectively, of the signals we’re searching for.

In summary, the framework for hypothesis testing comprises:

1.

Defining a test statistic $t$ to map data $x$ (in our example, $x=(n,m)$) to a single number.

2.

Deriving the distribution of $t$ under the hypothesis being tested $p(t|H)$ by sampling from “toy” datasets assuming $H$.

3.

Quantifying the compatibility of the observed data $x_{\mathrm{obs}}$ with $H$ with the $p$-value or significance $Z$ of $t_{\mathrm{obs}}$ relative to $p(t|H)$.

This $p$-value / significance is what we then use to decide whether or not to exclude $H$. A particularly important special case of this, as discussed above, is testing the background-only hypothesis when trying to discover a signal.

8.2.3 Confidence intervals and limits

Confidence intervals using the Neyman construction

Next, we discuss going beyond hypothesis testing to setting intervals and limits for parameters of interest. The machinery from Section 8.2.2 can be extended straightforwardly to extracting “confidence intervals” for our parameters of interest (POIs): a range of values of the POIs that are allowed, based on the experiment, at a certain “confidence level” (CL), e.g. 68% or 95%. Very similar to the idea of the significance level, the CL is defined such that if we were to repeat the experiment many times, a 95%-confidence-interval must contain, or cover, the true value of the parameter 95% of the time.

This can be ensured for any given CL by solving Eq. 8.2.14 for a $p$-value of $1-\mathrm{CL}$:

where $s_{-}$ and $s_{+}$ are the lower and upper limits on $s$, respectively.

This can be solved by scanning $s$ and finding the values of $s$ for which the RHS $=1-\mathrm{CL}$, as demonstrated in Figure 8.10 for the experiment in Example 8.2.1 ($n_{\mathrm{obs}}=20,m_{\mathrm{obs}}=5$). This procedure of inverting the hypothesis test by scanning along the values of the POIs is called the “Neyman construction”.

One subtlety to remember is that, in principle, we should also be scanning over the nuisance parameters ($b$) when estimating the $p$-values. However, this would be very computationally expensive so in practice, we continue to use $b=\widehat{\widehat{b}}(s)$, to always (approximately) maximize the agreement with the hypothesis. Ref. [58] calls this trick “profile construction”.

Upper limits

Typically if a search does not have enough sensitivity to directly observe a new signal, we instead quote an upper limit on the signal strength. This is similar in practice to the Neyman construction for confidence intervals, solving Eq. 8.2.17 only for the upper boundary. However, an important difference is that when setting upper limits, we have to modify the test statistic so that a best-fit signal strength greater than the expected signal ($ŝ>s$) does not lower the compatibility with $H_s$:

The upper limit test statistic $\tilde{q}(s)$ is set to $0$ for $ŝ>s$ so that this situation does not contribute to the $p$-value integral in Eq. 8.2.14. Figure 8.11 demonstrates this, and the difference between ${\tilde{t}}_s$ and ${\tilde{q}}_s$, for different sample observations.

Note that (as one may expect from Figure 8.11) the distribution $p({\tilde{q}}_s|s)$ no longer behaves like a standard ${\chi }^2$ but, instead, as a “half-${\chi }^2$”. This is essentially a ${\chi }^2$ plus a delta function at 0 (since, under the signal hypothesis, on average there will be an over-fluctuation half the time, for which ${\tilde{q}}_s=0$), as shown in Figure 8.12.

We can now revisit Example 8.2.1 to set an upper limit on $s$ rather than a confidence interval (Figure 8.13). $p({\tilde{q}}_s|s)$ is shifted to the left with respect to $p({\tilde{t}}_s|s)$; hence, the upper limit of $24$ is slightly lower than the upper bound of the 95% confidence interval we derived using ${\tilde{t}}_s$.

The CL${}_s$ criterion

We now introduce two conventions related to hypothesis testing and searches in particle physics. Firstly (the simple one), the POI $s$ is usually re-parametrized as $s\to \mu \cdot s$, where $\mu$ is now considered the POI, referred to as the “signal strength”, and $s$ is a fixed value representing the number of signal events we expect to see for the nominal signal strength $\mu$ of $1$. For the example in Figure 8.13, if we expect $s=10$ signal events, then we would quote the upper limit as $24∕s=2.4$ on $\mu$ at 95% CL.

The second, important, convention is that we use a slightly different criterion for confidence intervals, called “CL${}_s$”. This is motivated by situations where we have little sensitivity to the signal we’re searching for, as in the below example.

This is an extremely aggressive limit on $\mu$, where we’re excluding the nominal $\mu =1$ signal at a high confidence level. Given the complete lack of sensitivity to the signal, this is not a sensible result. The CL${}_s$ method solves this problem by considering both the $p$-value of the signal + background hypothesis $H_s$ (referred to as $p_{s+b}$ or just $p_{\mu }$ for short), and the $p$-value of the background-only hypothesis $H_0$ ($p_b$), to define a new criterion:

In cases where the signal region is completely background-dominated, the compatibility with the background-only hypothesis should be high, so $p_b\sim 1$ and, hence, $p_{\mu }^{\prime }$ will be increased. On the other hand, for more sensitive regions, compatibility should be lower $\Rightarrow p_b\sim 0$ and $p_{\mu }^{\prime }\sim p_{\mu }$.

To be explicit, here

where we should note that:

1.

We’re looking at the distribution of ${\tilde{t}}_s$ — not ${\tilde{t}}_0$ — under the background-only hypothesis, since the underlying test is of $H_s$, not $H_0$; and

2.

We’re integrating up to ${\tilde{t}}_{\mathrm{obs}}$, unlike for $p_s$, because lower $\tilde{t}$ means greater compatibility with the background-only hypothesis.

The effect of the CL${}_s$ criterion is demonstrated in Figure 8.14 for Examples 8.2.1 and 8.2.2. In the former, the background-only distribution is shifted to the right of the $s+b$ distribution. This indicates that the experiment is sensitive to $\mu$ and, indeed, we find $p_{\mu }^{\prime }\sim p_{\mu }$. In Example 8.2.2, however, the search is not sensitive and, hence, the background-only and $s+b$ distributions almost completely overlap, meaning $p_b\sim 1$ and $p_{\mu }^{\prime }>>p_{\mu }$.³

Finally, if we repeat the Neyman construction using the CL${}_s$ criterion $p_{\mu }^{\prime }$ instead of $p_{\mu }$ for Example 8.2.2, we can find an upper limit of $\mu \le 1.2$ at 95% CL, which is indeed a looser, more conservative, upper limit. The upper limit for Example 8.2.1 remains unchanged at $\mu \le 2.4$, as we would expect.

8.2.4 Expected significances and limits

Expected significance

The focus so far has been only on evaluating the results of experiments. However, it is equally important to characterize the expected sensitivity of the experiment before running it (or before looking at the data).

One way is to calculate the significance with which we expect to exclude the background-only hypothesis if the signal were, in fact, to exist. Practically, this means we are testing $H_0$ and, hence, need $p({\tilde{t}}_0|\mu =0)$ as before. However, now we also need the distribution of the test statistic ${\tilde{t}}_0$ under the background + signal hypothesis $p({\tilde{t}}_0|\mu =1)$. Then, by calculating the significance for each sampled ${\tilde{t}}_0$ under $H_{\mu =1}$, we can estimate the distribution of expected significances. This is illustrated for Example 8.2.3 in Figure 8.15.

Importantly, we usually quote the median of this distribution as the expected significance, since the median is “invariant” to monotonic transformations (i.e., the median $p$-value will always correspond to the median $Z$ as well, whereas the mean $p$-value will not correspond to the mean $Z$). Similarly, we quote the $16\%∕84\%$ and $2\%∕98\%$ quantiles as the $\pm 1\sigma$ and $\pm 2\sigma$, respectively, expected significances. These quantiles correspond to the cumulative probabilities for a standard Gaussian (Figure 8.16). For Example 8.2.3, we thus find the median expected significance to be 1.83.

Note that instead of converting each sampled $\tilde{t}$ under $H_{\mu =1}$ into a significance and finding the median of that distribution, as in Figure 8.15 (right), we can take advantage of the invariance of the median and directly use the significance of the median $\tilde{t}$ under $H_{\mu =1}$ (Figure 8.15, left). We will do this below for the expected limit.

Expected limits

The other figure of merit we care about in searches is the upper exclusion limit set on the signal strength. To derive the expected limit, we do the opposite of the above and ask, if the signal were not to exist, what value of $\mu$ would we expect to exclude at the 95% CL.⁴

This means we need:

1.

The distribution $p({\tilde{q}}_{\mu }|\mu )$ to solve for ${\mu }^{+}$ in Eq. 8.2.17 and be able to do the upper limit calculation (as in Section 8.2.3.);

2.

$p({\tilde{q}}_{\mu }|0)$ to get the median (and other quantiles’) expected ${\tilde{q}}_{\mu }^{\mathrm{obs}}$ for different signal strengths under the background-only hypothesis; and, furthermore,

3.

To scan over the different signal strengths to find the $\mu$ that results in a median $p$-value of $0.05$ — or, rather, $p_{\mu }^{\prime }$-value (Eq. 8.2.19), since we’re using the CL${}_s$ method for upper limits (Section 8.2.3.).

First, let’s look at the first two steps for just the $\mu =1$ signal strength in Example 8.2.3. These steps are similar to, and essentially an inversion of, the procedure for the expected significance: we’re now finding the $p_{\mu =1}^{\prime }$-value with respect to the signal + background hypothesis, for the median ${\tilde{q}}_{\mu }$ sampled under the background-only hypothesis. This is demonstrated in Figure 8.17.

The key difference with respect to calculating the expected significance is step 3, in which this procedure has to be repeated for a range of signal strengths to find the value that gives a median (and $\pm 1\sigma ,\pm 2\sigma$ quantile-) $p_{\mu }^{\prime }$ of $0.05$. This is thus the minimum value of $\mu$ that we expect to be able to exclude at 95% CL, as shown in Figure 8.18.

Thus, we have our expected limits. The right plot of Figure 8.18 is colloquially known as a “Brazil-band plot”, and is the standard way of representing limits. For example, Figure 8.19 is the corresponding plot by ATLAS for the Higgs discovery (scanning over the Higgs mass).