Machine learning and multivariate goodness of fit

Multivariate goodness-of-fit and two-sample tests are important components of many nuclear and particle physics analyses. While a variety of powerful methods are available if the dimensionality of the feature space is small, such tests rapidly lose power as the dimensionality increases and the data inevitably become sparse. Machine learning classifiers are powerful tools capable of reducing highly multivariate problems into univariate ones, on which commonly used tests such as $\chi^2$ or Kolmogorov-Smirnov may be applied. We explore applying both traditional and machine-learning-based tests to several example problems, and study how the power of each approach depends on the dimensionality. A pedagogical discussion is provided on which types of problems are best suited to using traditional versus machine-learning-based tests, and on the how to properly employ the machine-learning-based approach.


Introduction
Multivariate goodness-of-fit and two-sample tests are important components of many nuclear and particle physics analyses. For example, an unbinned maximum likelihood fit of a probability density function (PDF) to a data sample is commonly performed in such analyses. The maximum likelihood value cannot be used to determine the level of agreement between the fit PDF and the data; therefore, some other goodness-of-fit test must be employed for this important task. A related problem is determining whether two data samples -of which, one may be simulated -share the same parent PDF, where the PDF itself may be unknown, i.e. no parametric expression of the PDF is required. The most widely used statistic in nuclear and particle physics analyses for both goodness-of-fit and two-sample testing is the χ 2 [1], though many other tests are also applicable to both univariate and multivariate problems (see, e.g., Refs. [2,3]). While a variety of powerful tests exist when the dimensionality of the feature space is small, all such tests rapidly lose power as the dimensionality increases and the data inevitably become sparse; this is known as the curse of dimensionality [4].
Machine learning (ML) classifiers are powerful tools capable of reducing highly multivariate problems into univariate ones, on which commonly used tests such as χ 2 or Kolmogorov-Smirnov (KS) [5,6] may be applied. To the best of our knowledge, this idea was first proposed in the literature in Ref. [7]. It has since become a popular method in the data-science toolkit (see, e.g., Ref. [8] and references therein, especially Refs. [9,10]) for the following reasons: ML classifiers are well suited to problems where the dimensionality is large; powerful ML classifier algorithms are now available in easy-to-use open-source software packages, and in several programming languages; the previous two points have made the usage of ML methods commonplace in classification problems, which in turn has led to analysts developing both skill and comfort using them; and, finally, once the ML classifier has reduced the problem into a univariate one, quantifying the level of agreement with a p-value is simple, e.g., using the well known χ 2 or KS test.
Despite the fact that ML algorithms are now commonly used in nuclear and particle physics analyses, e.g., to classify events as signal or background, including in online data-selection trigger algorithms [11], most of the community appears to be unaware of the possibility of employing ML classifiers for goodness-of-fit or two-sample testing. In this article, we explore applying both traditional and ML-based methods to various example problems, and study how the power of each depends on dimensionality. A pedagogical discussion is provided on which types of problems are best suited to the use of traditional versus ML-based approaches to goodness-of-fit and two-sample testing, and on how to properly employ the machine-learning-based approach. The primary goals of this article are: • to present the ML-based approach to the nuclear and particle physics community in a clear and concise manner, using language familiar to members of our community; • to provide a step-by-step guide to proper application of the ML-based approach, including the inclusion of systematic uncertainties; • and to impart in the reader some intuition of when to expect ML-based goodness-of-fit tests to outperform traditional methods.
This article is structured as follows: Sec. 2 discusses traditional testing methods; an overview of applying ML algorithms for these tasks is presented in Sec. 3; example problems are presented in Sec. 4; and a summary is provided in Sec. 5.

Traditional Methods
This section provides a brief overview of traditional, i.e. non-ML-based, methods for performing goodness-of-fit and two-sample testing. For a more complete discussion, we refer the interested reader to Refs. [2,3]. Throughout this article, we denote the d-dimensional vector of event features as x and PDFs as f ( x), each with subscripts where necessary for clarity.

Goodness-of-Fit & Two-Sample Testing
The goal of a goodness-of-fit test is to obtain a p-value for the hypothesis f ( g., f 0 may have been determined by performing an unbinned maximum likelihood fit to the data. The p-value is the probability of observing a test statistic that corresponds to lesser agreement than the actual observed value when the null hypothesis f = f 0 is true. Therefore, if f = f 0 , the p-values are uniformly distributed on the interval (0, 1). A related problem is determining whether two data samples, { x i } n a i=1 and { x j } n b j=1 , share the same parent PDF, i.e., the goal is to determine the p-value for the hypothesis f a = f b , where parametric expressions of f a and f b may not be available, e.g., searching for CP violation in Dalitz decays [12]. In each case, one ideally wants to use a test that has power against all alternatives in the asymptotic limit n (a,b) → ∞. N.b., one can recast a goodness-of-fit test into a two-sample test by generating a simulated (Monte Carlo) data sample from the PDF f 0 , then performing a two-sample test on the simulated and experimental data samples. This is a common approach in cases where an analytic expression for the detector response does not exist, and is a valid way of assessing the goodness of fit in all cases.

Univariate Methods
For the univariate case, each observation involves a single feature x i = x i . The most commonly used one-dimensional tests in nuclear and particle physics are the χ 2 and KS tests; however, a variety of other powerful tests are available for univariate problems, e.g., the Anderson-Darling (AD) [13] and Cramer-von Mises (CvM) [14,15] tests. We assume that the reader is familiar with both the χ 2 and KS tests and omit any detailed discussion on them; however, we do feel it is important to note that many algorithms for obtaining p-values from a two-sample KS test assume that the underlying distributions are continuous. The response of some ML classifiers is discrete, so some care needs to be applied here. If the fraction of repeated x i values is small, then using an algorithm that assumes x is continuous will give a good estimate of the p-value; however, if the number of ties is large, the p-value must be obtained using another method, e.g., the permutation test [16].

Multivariate Methods
As with the univariate case, there are many methods available for performing multivariate goodnessof-fit and two-sample tests; however, most are (in some way) based on local-density, and any density-based method loses power rapidly as d increases (even power techniques like the energy test [17]). For a detailed review of traditional multivariate tests, see Refs. [2,3]. The χ 2 test is applicable in any number of dimensions, but it is important to define a binning scheme that avoids producing too many low-occupancy bins when d is large. When comparing the χ 2 and ML-based methods below, we employ adaptive (probability) binning [18], where the binning scheme is defined using the pooled { x i } n a i=1 + { x j } n b j=1 data sample such that all bins in both the a and b samples are expected to have equal population if f a = f b when n a,b → ∞. The number of bins is optimized separately for each example problem studied in Sec. 4, and we note that this number is in all cases less than one might naively expect.

Machine Learning Methods
ML classifiers are powerful tools capable of reducing highly multivariate problems into univariate ones. This section describes the procedure for using a ML algorithm to reduce the problem down to a single dimension, then applying a well-known univariate two-sample test to obtain the pvalue. Inclusion of systematic uncertainties and some details about specific ML algorithms is also provided. N.b., a related topic is the connection between ML classifiers and likelihood ratio testing, see Ref. [19] for detailed discussion and Ref. [20] for implementation in software.

Strategy
The basic strategy employed is to use a ML classifier to reduce the d-dimensional problem down to a one-dimensional one, then to apply a well-known univariate two-sample test to obtain the p-value. Specifically, the steps are: • each sample is labeled and split into training, validation, and testing subsamples; • from the training samples, the ML algorithm learns a model m( x) designed to classify which sample an event belongs to based on the measured values of its x features (m( x) is often referred to as the score or response); • from the validation samples, the ML algorithm hyperparameters are optimized; • from the testing samples, two univariate distributions {m( x i )} are constructed, one for each label value; • and, finally, a univariate two-sample test, e.g., the χ 2 or KS test, is performed on the two testing-sample {m( x i )} distributions to obtain a p-value.
To instead perform a goodness-of-fit test, the same procedure is followed but where one of the data samples is simulated by sampling from the test PDF f 0 . This procedure is simple to implement as it involves sequentially performing two tasks that physicists already know how to do: (1) training a ML classifier, and (2) obtaining a p-value from a comparison of two one-dimensional distributions. N.b., if the sample sizes are limited, k-fold cross validation can be used.

Systematic Uncertainties
A simple method for including systematic uncertainties involves first choosing to bin the onedimensional {m( x i )} distributions and selecting the χ 2 two-sample test to obtain the p-value. This facilitates adding bin-by-bin systematic uncertainties in the standard way one includes them for any χ 2 test (e.g., see Eq. 4.1 below). The best approach for determining the systematic uncertainty in each {m( x i )} bin will depend on the specific problem, but will be similar to assigning systematics for any binned measurement. Some useful techniques include: weighting the events in the testing samples to account for systematic effects, then producing weighted {m( x i )} distributions where the shifts in the bin contents relative to the unweighted distributions are used to determine the uncertainty; and generating alternative data samples that include some systematic effect(s), then repeating the entire procedure, including the training, where again the differences in the {m( x i )} distributions are used to determine the uncertainties.

Algorithms Considered
Since our focus is on goodness-of-fit and two-sample tests, we restrict our study to ML algorithms that are easy for physicists to use and do not require significant CPU resources. 1 Even with this restriction, a large number of ML classifiers are now available in easy-to-use open-source software packages. In our study, we consider: BDT boosted decision trees using either AdaBoost [22], as implemented in SCIKIT-LEARN [23], or XGBoost [24], as implemented in REP [25], and a tunable number of trees (other hyperparameters are also tuned that are specific to the boosting algorithms); ANN artificial neural networks that are multilayered, as implemented in KERAS [26] and SCIKIT-LEARN, with zero or one hidden layer and a tunable number of neurons per layer; SVM support vector machines using the radial basis function kernel as implemented in SCIKIT-LEARN, with tunable C and γ parameters [23].
All ML algorithms are trained on the training subsamples, the hyperparameters are tuned on the validation subsamples using the Bayesian optimization framework as implemented in the SPEARMINT package [27], and the performance is evaluated on the testing subsamples (these are all independent data samples). For the sake of readability, only the best-performing ML algorithm from each category is shown for each example problem in Sec. 4.

Example Problems
In this section, we consider three example problems. The first two are chosen for their pedagogical value, and because each is easy to implement for any dimensionality. For these, we study how the performance of the χ 2 and ML-based tests depends on d. The third problem is more realistic, and involves searching for a two-body particle decay. For each problem, ensembles of data sets are generated each with a sample size of 10k events. 2 The notation used in this section is that in all cases the two data samples being compared have parent PDFs labeled by f a and f b . The a data sample will always originate from f 1 , i.e. f a = f 1 . When the null hypothesis is true, f b = f 1 ; however, when it is not true, then f b = f 2 . The power of each method is assessed by studying how frequently it is able to reject the null hypothesis ( f a = f b ) at 95% confidence level. 1 Our motivation here is pragmatic. If performing a test requires purchasing commercial software, dedicating significant effort to learning how to implement the ML algorithm in software, or if substantial CPU power is needed to adequately train the ML algorithm, we do not believe that such a method would be adopted by a substantial fraction of the physics community for use as a goodness-of-fit or two-sample test. Of course, the reader is free -even encouraged -to disagree, and could instead apply the ML-based procedure using a deep-learning neural network library, e.g. TENSORFLOW [21], or some other cutting-edge ML algorithm. The procedure presented in this paper is independent of which ML classifier is used. 2 We use 5-fold validation, where each ML algorithm is trained on 8k events and tested on the remaining 2k. This process is repeated 5 times such that all 10k events are utilized in both training and testing, but always in an unbiased way. To reduce the CPU time required for training, the hyperparameters are optimized for each algorithm only once per ensemble, then used for data sets in the ensemble. This should not need to be done in a real-world analysis, where only a single data sample is analyzed.
We universally assume a systematic uncertainty of 1% in all regions of feature space, largely just to demonstrate how to incorporate systematic uncertainties. Therefore, the two-sample χ 2 , where each data sample has the same number of events, is defined as where a k (b k ) denotes the yield of sample a(b) in the k th bin. The systematic uncertainties are implemented by randomly sampling a k (b k ) from the Gaussian functions G (µ = n k a,b , σ = 0.01 · n k a,b ), where n k a,b is the generated event count in bin k. This is done both for the multivariate χ 2 tests, and for the univariate χ 2 tests performed on the ML-response distributions (of course, implementing systematic uncertainties by smearing the bin contents is not required in a real-world analysis).

Gaussian
First we consider a d-dimensional Gaussian distribution, where the values {x k } d k=1 are the observed features for each event. Events are sampled from the d-dimensional PDF where σ k = 1 for f 1 and σ k = 0.95 for f 2 . Figure 1 shows example data sets sampled from f 1 and f 2 for the case d = 2, as this is the only dimensionality that is easily displayed. Figure 2 shows the optimal adaptive binning for an example two-dimensional data set, while Fig. 3 shows example optimal-binned distributions for d = 2, 6, and 10; the binned spectra are shown flattened into a single dimension by plotting versus bin number to aid in visually comparing the a and b samples. Example ML-response distributions are shown in Fig. 4. First, we verify that each test produces valid p-values under the null hypothesis by generating an ensemble of 2 × 100 data sets where f a = f b = f 1 . Figure 5 shows that both the multivariate χ 2 and ML-based methods reject the expected ≈ 5% of data samples when the null hypothesis is true. Next, we generate an alternative set of 100 data sets with f b = f 2 . Figure 5 shows that the MLbased tests outperforms the χ 2 test for d > 2. In this example, where the PDFs are hyperspherical, the ANN and SVM outperform the BDT. This is not surprising, since decision trees are constructed through a series of one-dimensional (rectangular) splittings of the data. The main conclusion drawn from this example problem is that, as expected, the power of the χ 2 test decreases rapidly as d increases, while the power of the ML-based tests is much less sensitive to dimensionality. If the underlying model is known, we can design a powerful feature for distinguishing between f 1 and f 2 . Figure 6 shows the d-dimensional Euclidean distance from the origin of each measurement for d = 10. This feature optimally reduces this d-dimensional problem into a one-dimensional one. Applying the one-dimensional χ 2 test to the reduced problem provides greater power than applying the ML-based approach to the d-dimensional one. Indeed, this is a common approach in physics analyses, e.g., when studying particle decays the invariant mass of the decay products is used to reduce all of the kinematic information into a single feature. 3 Our d-dimensional Gaussian  example problem is effectively a one-dimensional one, and so the adaptive-binned χ 2 is a powerful option if we treat it as such; however, if we are unable to design a feature to fully capture all of the relevant information, which is likely the case when the underlying model is unknown, then the ML-based approach is more powerful than traditional ones.

Oscillatory
Next we consider a d-dimensional oscillatory problem where the observed features are x = θ = {θ k } d k=1 and the PDF is defined according to where n k = 5 for all k in f 1 , but for f 2 the value of n 1 = 6 is used (all other n k>1 = 5, see Fig. 7). The results of applying χ 2 and ML-based two-sample tests are shown in Figs. 8 and 9. Here, the underlying rectangular nature of the BDT is advantageous, whereas our chosen set of ANN and SVM hyperparameters do not lead to good performance of these algorithms on this problem. This example demonstrates that in some cases, one ML algorithm may vastly outperform some others. Studying the performance of different ML-based tests on simulated data samples can reveal such situations. We note that this problem can also be reduced to a single dimension by projecting onto the θ 1 axis, if one knows that the two models are the same in all other dimensions. The fact that our first two example problems possess this feature should not be taken to mean that most realworld problems are inherently one-dimensional. Rather, the simplest problems to implement for any dimensionality d are inherently one-dimensional.

Particle Decays
For the final example problem, we consider a simple real-world physics use case of searching for a particle of mass m(X) decaying into two equal-mass lighter particles X → αβ . The energy and  We have a background-only data sample (in a real search, this could be from simulation, same-sign candidates in data, etc.) and another data sample in which we want to look for evidence of an X → αβ signal. First, we apply the ML-based and adaptive-binned χ 2 tests to the 8-dimensional problem, where the fraction of the search data sample that arises from X → αβ decays is varied from zero to 10%. The results are shown in Fig. 10. As expected, the ML-based methods outperform the χ 2 test; however, when studying particle decays, we know that calculating the invariant mass is an effective way to reduce the kinematic information into a single feature. Figure 10 also shows that, for this simple problem, the one-dimensional χ 2 test on the m(αβ ) distribution outperforms the ML-based methods performed on the 8-dimensional {E α , p α , E β , p β } feature space -at least for the ML algorithms employed in this study. More advanced deep-learning algorithms may be competitive with the human-led feature design approach of using the invariant mass here [28].

Summary
Multivariate goodness-of-fit and two-sample tests are important components of many nuclear and particle physics analyses; however, the tests commonly used in such analyses are only powerful when the dimensionality is small. Machine learning classifiers are powerful tools capable of reducing highly multivariate problems into univariate ones, on which commonly used tests such as χ 2 or Kolmogorov-Smirnov may be applied. We explored applying both traditional and machinelearning-based tests to several example problems, and studied how the power depends on dimensionality. For up to about three dimensions, the power of the adaptive-binned χ 2 test is comparable to that of ML-based tests; however, for higher dimensionality the ML-based approach is far superior to the χ 2 test. A caveat to the previous statement is that if some information about the underlying model is known, it may be possible to construct features that reduce the effective dimensionality of the problem. Indeed, such human-led feature design is common in physics analyses, e.g., when studying particle decays the invariant mass is used to reduce all of the kinematic information into a single feature. In such cases, the problem may effectively be a low-dimensional one and the adaptivebinned χ 2 is a powerful option. When one or two features cannot be designed to fully capture all of the relevant information -which is likely the case when detector-response features are included in the problem -ML-based tests are expected to be more powerful than traditional ones. We encourage the community to investigate adopting this approach to goodness-of-fit and two-sample testing for problems where the dimensionality is large and cannot be reduced without significant loss of information.