Statistical analysis of event classification
in experimental data

The analysis of experimental data – e.g. events created at a particle collider and recorded in a detector – often requires selecting those rare events with a signature of interest (the “signal”) and rejecting the bulk of “background”. This is performed by some “test”, modelled according to the signature, in order to make a decision (positive or negative) of whether the hypothesis “event is signal” is true or false.

In practice, such a test is never perfect: some signal events may wrongly be rejected (error type I), and some background may wrongly be selected (error type II). It is therefore essential to know the fraction of “wrong events” in the selected and rejected data samples, which are characteristic values of the test being used.

These values are not known a priori, but may be obtained by applying the very same test on two samples of Monte Carlo generated data, simulating as accurate as possible pure signal and pure background, respectively. The test’s “characteristic values” (error types I and II and their complements: sensitivity and specificity) are conditional probabilities, represented by the $2\times 2$ “contingency matrix”.

The test criteria depend on parameters which influence the selection and rejection characteristics. Optimization may be achieved by defining an appropriate “figure of merit”, to be maximized by varying the test parameters.

The conditional posterior probabilitiy of a selected event being genuine signal, provided the test result is positive, can be calculated by Bayesian inference. It depends on the prior probability (prevalence) which is in general unknown, however, can be estimated by using Monte Carlo information together with real data.

Part I presents a tutorial on the basics (Section 2), the signal and background simulations for obtaining the test charcteristics as encoded in the contingency matrix (Section 3), and their use in Bayesian inference (Section 4). The scenario of the base sample being created by two independent Poisson processes for signal and background is regarded in Section 5. Postulating these numbers to be fixed, estimates for the number of selected and rejected events are given (Subsection 5.1); test optimization and the choice of a “figure of merit” are discussed in Subsection 5.2. Section 6 presents a method for estimating the unknown numbers of signal and background events from the known numbers of selected and rejected events, using Monte Carlo information; estimates of the Bayesian prior probabilities follow straight-forward. Section 7 shows how to assess the quality of selected data by its purity or its fraction of background contamination.

Part II follows up with a comprehensive and rigorous statistical analysis. Section 8 recapitulates the basic assumptions and the notation. Bayesian estimation is the topic of Section 9. The posterior distribution of the fraction of events classified as signal is computed for various conjugate priors, both “uninformative” and informative. A simple affine transformation yields the posterior mean and the posterior variance of the unknown signal content of the sample. Frequentist estimation, on the other hand, presupposes a virtual parent population from which the observed sample is drawn (Section 10). The frequentist estimators and their covariance matrices are calculated for three choices of the parent population: fixed sample size (Subsection 10.1), Poisson-distributed sample size (Subsection 10.2), and sample size drawn from a mixture of Poisson distribution (Subsection 10.3). The pros and cons of the two approaches are summarized in the final Section 11.

Part I: The tutorial

It is assumed that some test is available which aims at determining, for each element $\mathfrak{e}_{i}\in\mathfrak{D}$, whether it does comply with hypothesis $\mathcal{H}$ and therefore will be selected and put into sub-sample $\mathfrak{S}^{\ast}$; otherwise it will be rejected and put into sub-sample $\mathfrak{B}^{\ast}=\mathfrak{D}\setminus\mathfrak{S}^{\ast}$. The “positive test result” $\mathcal{T}\equiv\mathfrak{e}_{i}\in\mathfrak{S}^{\ast}$ and the “negative test result” $\bar{\mathcal{T}}\equiv\mathfrak{e}_{i}\in\mathfrak{B}^{\ast}$ are defined accordingly.

However, any test is inherently not perfect but prone to errors of type I and II, the rates of which must be evaluated.

In order to achieve this, two separate samples are created by Monte Carlo simulation: $\mathfrak{S}^{\prime}$ containing $N^{\prime}_{S}$ pure signal events, and $\mathfrak{B}^{\prime}$ containing $N^{\prime}_{B}$ pure background events. Their ratio may be chosen arbitrarily; often $N^{\prime}_{S}/N^{\prime}_{B}\gg R_{S/B}$.

The simulated to real background ratio is

and its inverse $1/r_{B}$ may serve as the normalization factor for scaling simulated to real background, as used in Section 7.

Similarly, Eq. (18) asks: when applying the test, which is the conditional posterior probability $p_{S-}$ of $\mathfrak{e}_{i}$ being signal ($\mathcal{H}$ = true) if the test result is negative ($\mathcal{T}$ = false)? Applying Bayes’ theorem to this case yields:

Reverting to the variables of Eqs. (2, 9, 11, 13), Eq. (24) becomes

Similarly to above, the posterior probability $p_{S-}$ depends only on the prior probability $p_{S}$ and on the negative likelihood ratio $\ell_{-}=\alpha/\eta$.

Calculation of the remaining conditional posterior probabilities $p_{B+}$ and $p_{B-}$ of $\mathfrak{e}_{i}$ being background ($\mathcal{H}$ = false) if the test result is positive ($\mathcal{T}$ = true) or negative ($\mathcal{T}$ = false), respectively, is trivial by using Eqs. (19–20).

A priori information are the likelihood ratios $\ell_{+}$ and $\ell_{-}$ which are obtained from the Monte Carlo simulations, see Section 3. The prior probabilities $p_{S}$ and $p_{B}$ are the signal and background fractions, Eq. (2). They are determined by Eq. (4) if the signal to background ratio $R_{S/B}$ is known; if not, estimates of them can be calculated from the real data sample by a method presented in Section 6.

In many experiments, notably those performed by particle collisions, the real data sample $\mathfrak{D}$ has been stochastically accumulated by two independent Poisson processes for signal and background events, respectively.⁶⁶6   In real experiments, the raw data are further affected by triggers of data acquisition, the detector acceptance, reconstruction efficiencies, pre-selection for the analysis proper, etc. These complex aspects, which are beyond the scope of this paper, are discussed e.g. in [5, 6]. Thus, their numbers $N_{S}$ and $N_{B}$ behave as random variables with expectation values

where $p_{S}$ and $p_{B}$ are fixed values, determined only by the physics laws underlying the production of signal and background events; they are a-priori unknown, but their values can be estimated as shown in section 6.

The base sample $\mathfrak{D}=\mathfrak{S}\cup\mathfrak{B}$ thus obtained may be regarded, from a frequentist’s point-of-view, as one stochastic draw from a virtual population, consisting of many repetitions of the experiment at same conditions.

The marginal standard deviations and the correlation ⁷⁷7   The correlation $\rho\,(N_{S},N_{B})$ is a characteristic of these distributions, hence it is a constant. It must not be confused with the empirical correlation $r\,(N_{S},N_{B})$ which is a random variable, dependent on the actual draw, however with expectation value $\left<\,r\,\right>=\rho\,(N_{S},N_{B})=0$. As can be shown by simulation with an increasing number of draws, plots of the $r$-distribution approach the Dirac function $\delta\,(r)$, thus illustrating the fact that $r\,(N_{S},N_{B})$ is a consistent estimator of $\rho\,(N_{S},N_{B})$. of the bivariate distribution of $N_{S}$ and $N_{B}$ are as follows:

Hence, for the total number $N=N_{S}+N_{B}$

However, for dealing with the following selection, $\mathfrak{D}$ is postulated to be a base sample of fixed albeit unknown sizes $N_{S}$ of $\mathfrak{S}$ (signal) and $N_{B}$ of $\mathfrak{B}$ (background), and their sum $N=N_{S}+N_{B}$ being the fixed and known size of $\mathfrak{D}=\mathfrak{S}\cup\mathfrak{B}$.

A rigorous analysis is deferred to part II of this study.

The fractions $p_{S}$ and $p_{B}=1-p_{S}$ must be obtained from independent prior knowledge. If the signal to background ratio (SBR) is known, i.e. $R_{S/B}=\left<R_{S/B}\right>\pm\sigma(R_{S/B})$, their values are determined by Eq. (4), and error propagation yields

Estimates of $N_{S}$ and $N_{B}$ are derived straight-forward from Eq. (2).

In the following we assume that SBR is unknown. In this non-trivial case the fractions $p_{S}$ and $p_{B}$ are unknown parameters, however, their estimation is possible by using as prior knowledge the test’s characteristic values $\left[\,\varepsilon,\alpha,\beta,\eta\,\right]$, applied to the real data sample $\mathfrak{D}$ at hand. The method is described in detail below.

The number of real data events selected and rejected are given by $E_{S}=\left<E_{S}\right>\pm\sigma(E_{S})$ and $E_{B}=\left<E_{B}\right>\pm\sigma(E_{B})$, respectively, with $E_{S}+E_{B}=N$. These numbers can be used for calculating estimates of the true number of signal and background events, $\tilde{N}_{S}$ and $\tilde{N}_{B}$, respectively, by re-writing Eqs. (30–31) as

or in matrix notation (see Eq. (15)) as

Using Eqs. (14, 16) yields the solutions ⁹⁹9   $\det\,(\bm{A})=\delta=\varepsilon-\beta>0$ (see Eq. (16)) ensures non-singularity. The simpler formulae $\tilde{N}_{S}\approx E_{S}/\varepsilon$, $\tilde{N}_{B}\approx N-E_{S}/\varepsilon$, as often used instead, are valid only for a negligible contamination rate $\beta\ll 1$.

obeying $\tilde{N}_{S}+\tilde{N}_{B}=N$. Error propagation, observing $\rho\,(E_{S},E_{B})=-1$ (see Eq. (33)), yields the corresponding standard deviations as

with $\sigma(E_{S})=\sigma(E_{B})$ given by Eq. (32) and approximately setting $N_{S}\approx\tilde{N}_{S}$ and $N_{B}\approx\tilde{N}_{B}$. $\tilde{N}_{S}+\tilde{N}_{B}=N$ implies maximal anti-correlation:

and error propagation yields $\sigma(\tilde{N}_{S}+\tilde{N}_{B})=|\,\sigma(\tilde{N}_{S})-\sigma(\tilde{N}_{B})\,|=0$.

Note that the values of Eqs. (41–43) are derived empirically from the observed numbers $E_{S}$ and $E_{B}$. They must not be confused with Eq. (28) which represents characteristics of the distributions determined by the data acquisition.

Accurate estimates of the true number of signal and background events in the real data sample $\mathfrak{D}$, see Eqs. (39–42), are essential for data analyses based on event counting, e.g. measurement of particle decay branching fractions.

Investigation of the sample $\mathfrak{S}^{\ast}$, containing $E_{S}$ selected real events, requires an accurate knowledge of the expected numbers $\left<E_{S/S}\right>$ of pure signal events and $\left<E_{B/S}\right>$ of contaminating background events within this sample.¹⁰¹⁰10   Such an investigation would typically be summarized like “analysis performed on $E_{S}$ events, $\left<E_{B/S}\right>$ of which are expected to be background”.

Their fractions w.r.t. all selected events, purity ¹¹¹¹11   The term “purity” is defined ambigously in the literature. and fraction of contamination, are relevant quantities for judging the data quality:

These fractions $x_{S}$ and $x_{B}=1-x_{S}$ are the probabilities of a selected event to be signal or background, respectively. The number $n_{S}$ ($n_{B}$) of signal (background) events within the sample $\mathfrak{S}^{\ast}$ of size $E_{S}$ is distributed according to the binomial pdf:

Let $\mathfrak{D}$ be a data set consisting of $N$ recorded events $\mathfrak{e}_{i},\;{}i=1,\ldots,N$. $\mathfrak{D}$ can be partitioned into the subsets $\mathfrak{S}$ containing the “signal” events and $\mathfrak{B}$ containing the “background” events:

The respective fractions of “signal”and “background” events are denoted by:

A binary indicator $I_{i}\in\{0,1\}$ is attached to each event _i such that

The actual values of the indicators are unknown a priori.

A binary classifier is a function $C$ that computes a function value $c_{i}=C(\mathfrak{e}_{i})\in\{0,1\}$ . If $c_{i}=1$, event $\mathfrak{e}_{i}$ is classified as “signal”; if $c_{i}=0$, event $\mathfrak{e}_{i}$ is classified as “background”. Such a classifier cannot be expected to perform perfectly. The probabilities of correct and wrong decisions can be summarized as follows:

where $\alpha\geq 0$ is the probability of classifying a “signal” event as being “background” (type I error), and $\beta\geq 0$ is the probability of classifying a “background” event as being “signal” (type II error). It is required that $\alpha+\beta<1$, see Eq. (16). The difference $\varepsilon-\beta=1-\alpha-\beta>0$ will be denoted by $\delta$.

The probabilities $\alpha$ and $\beta$ are unknown a priori, but can be estimated to arbitrary precision from an independent, sufficiently large simulated sample $\mathfrak{D}^{\prime}$, for which the classification of events into “signal”and “background” is known, see Section 3. In the following, $\alpha$ and $\beta$ will be considered as known constants.

The respective probabilities $q_{S},q_{B}$ of classifying an arbitrary event as “signal” or “background” are given by:

The respective numbers $E_{S}$ and $E_{B}$ of events classified by $C$ as “signal” and “background” are given by:

Conditional on $q_{S}$, $E_{S}$ is distributed according to the binomial distribution Bino($N,q_{S}$) with the probability function:

The posterior pdf of $q_{S}$, given an observation $E_{S}$, can be written down in closed form if the prior pdf $\pi(q_{S})$ is taken from the conjugate family of priors, which is the Beta family in this case [9].

If the prior distribution of $q_{S}$ is Beta($a,b$), the posterior pdf is given by:

The Bayes estimator, i.e., the posterior mean and its variance are equal to:

where $\langle X\rangle$ is the expectation of the random variable $X$ and $\sigma^{2}(X)$ is its variance.

If no prior information on $q_{S}$ is available, there are three popular “uninformative” priors ¹²¹²12   The maximum-entropy prior and Jeffrey’s prior are not really uninformative, as they assign a specific probability to any interval $q_{1}\leq q_{S}\leq q_{2}$ that is contained in $(0,1)$, see [10]. Neither is Haldane’s improper prior, as it can be considered as the limit of Beta($\varepsilon,\varepsilon$) for $\varepsilon\longrightarrow 0$, thus in the limit assigning the probability of 0.5 to both $q_{S}=0$ and $q_{S}=1$. in the unit interval:

1. 1.

   the maximum-entropy prior Beta($1,1$)=Unif($0,1$): $\pi(q_{S})=1,\ 0\leq q_{S}\leq 1$;

2. 2.

   Jeffrey’s prior Beta($\textstyle\frac{1}{2},\textstyle\frac{1}{2}$): $\pi(q_{S})\propto q_{S}^{-1/2}(1-q_{S})^{-1/2},\ 0<q_{S}<1$;

3. 3.

   Haldane’s improper prior Beta($0,0$): $\pi(q_{S})\propto q_{S}^{-1}(1-q_{S})^{-1},\ 0<q_{S}<1$.

Prior knowledge about $q_{S}$ can be incorporated via an informative prior. For example, assume that it is known that $q_{S}$ is about 0.4, with a standard uncertainty of about 0.1. This information can be encoded by a Beta prior Beta($a,b$) with mean 0.4 and variance 0.01:

Credible intervals for $q_{S}$ can be computed via quantiles of the posterior distribution. A symmetric credible interval $C$ with probability content of $1-\alpha$ is given by $C=\left[c_{1},c_{2}\right]$, where $c_{1}$ is the $\alpha/2$-quantile of the posterior pdf in Eq. 68 and $c_{2}$ is the $(1-\alpha/2)$-quantile.

An HPD interval $[h_{1},h_{2}]$ can be computed by solving the following system of non-linear equations:

where $f(h)=p(h|E_{S}=k)$ is the posterior pdf in Eq. 68 and $F(h)$ is its cumulative distribution function. A MATLAB code snippet that solves the system is shown in Fig. 2.