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 through Wilks' 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 their distributions are better approximated by the asymptotic ones; these improvements go under the general 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 to find 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 Barndorff-Nielsen [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 $n_d$ degrees of freedom. A new statistic is defined as $t' = tn_d/E[t]$, which by construction has mean $E[t'] = n_d$, as for the chi-square distribution. Further details and applications of these methods are described in Refs. [30, 31].

One may also treat model uncertainties in a Bayesian manner but then use the resulting model in a frequentist test. Suppose the uncertainty in a set of nuisance parameters $\nu$ is characterized by a Bayesian prior p.d.f. $\pi(\nu)$. This can be used to construct the marginal (also called the prior predictive) model for the data $\boldsymbol{x}$ and parameters of interest $\theta$,

$$P_m(x|\theta) = \int P(x|\theta, \nu)\pi(\nu) d\nu. \quad (40.50)$$

The marginal model does not represent the probability of data that would be generated if one were really 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 the experiment is carried out with new values of $\nu$, randomly sampled from $\pi(\nu)$. It is in effect an average of models each with a given $\nu$, where the average is carried out with respect to the prior p.d.f. $\pi(\nu)$.

The marginal model for the data $\boldsymbol{x}$ can be used to determine the distribution of a test statistic $Q$, which can be written

$$P_m(Q|\theta) = \int P(Q|\theta, \nu)\pi(\nu) d\nu. \quad (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, $L_{s+b}$, to that of background only, $L_b$. Often the likelihoods are evaluated with the profiled values of the nuisance 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. 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 search region. 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 can require 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 \times 10^{-7}$, corresponding to a $5\sigma$ or larger effect. Determining this with Monte Carlo thus requires generating and fitting a very large number of experiments, perhaps several times $10^7$. 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 the need 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 the peak position that is only defined under the signal model (there is no peak in the background-only model). Furthermore, suppose a test statistic $q_0$ is defined using the profile likelihood ratio, so that by 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 value $u$, then an approximation for the global p-value is found to be

$$p_{\text{global}} \approx p_{\text{local}} + \langle N_u \rangle, \quad (40.52)$$

where $\langle N_u \rangle$, which is much smaller than one in cases of interest, is the mean number of “upcrossings” (as defined in Ref. [34]) of the statistic $q_0$ above the level $u$ in the range of the nuisance parameter considered (e.g., the mass range). The value of $\langle N_u \rangle$ can be estimated from the number of upcrossings $\langle N_{u_0} \rangle$ above some much lower value, $u_0$, by using a relation due to Davis [35],

$$\langle N_u \rangle \approx \langle N_{u_0} \rangle e^{-(u-u_0)/2}. \quad (40.53)$$

By choosing $u_0$ sufficiently low, the value of $\langle N_u \rangle$ can be estimated by simulating only a very small number of experiments, or even from the observed data, rather than of the order $10^7$ needed if one is dealing with a $5\sigma$ effect.

### 40.3.3 Goodness of fit

At times one wants to quantify the level of agreement between the data $\boldsymbol{x}$ and a hypothesis $H_0$ 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 $t(\boldsymbol{x})$ as done in Sec. 40.3.2 with a significance test. Here, however, the goodness-of-fit statistic is defined using general considerations of what constitutes greater or lesser compatibility between data and hypothesis and not with reference to alternative hypotheses or power. Nevertheless, for a given statistic $t(\boldsymbol{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) can be investigated.

The Neyman-Pearson lemma (see Sec. 40.3.1) states that a test of $H_0$ that has maximum power with respect to an alternative $H_1$ should be based on a statistic that is monotonic in the likelihood ratio $P(x|H_1)/P(x|H_0)$. Therefore a goodness-of-fit statistic derived without reference to a specific $H_1$ will not in general attain this maximum power. And even if a given $t(x)$ happens to be optimal or nearly so for some alternative $H_1$, there may be other relevant alternatives for which its power is lower. In practice, goodness-of-fit statistics are defined so as to provide sensitivity to a broad class of alternatives, e.g., those that differ in the location of the bulk of the distribution's probability or 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 the method 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 discussion of 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 = (y_1, \dots, y_N)$ are used to test the hypothesis that the values are Gaussian distributed with given mean values $E[y_i] = \mu_i$ and variances $V[y_i] = \sigma_i^2$ (or covariance $V_{ij} = \text{cov}[y_i, y_j]$ 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 $\mu_i$ may be given as functions of some other parameters $\theta = (\theta_1, \dots, \theta_M)$, with $M < N$. The (co)variances are treated 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 $\chi^2(\theta)$ with respect to $\theta$, i.e., $t = \chi^2_{\min} = \chi^2(\hat{\theta})$, e.g., in the case of uncorrelated measurements

$$t = \sum_{i=1}^{N} \frac{(y_i - \mu_i(\hat{\theta}))^2}{\sigma_i^2}. \qquad (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 $y_i$ and predicted $\mu_i$. The $p$-value of the Gaussian hypothesis is therefore taken to be

$$p = \int_{\chi_{\min}^2}^{\infty} f(t; n_d) dt, \qquad (40.55)$$

where $f(t; n_d)$ is the $\chi^2$ p.d.f. and $n_d = 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 $\chi^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 Carlo calculation.

Often the data represent numbers of events in $N$ bins of a histogram, i.e., $n = (n_1, \dots, n_N)$. An important case is when $n$ follows a multinomial distribution with $n_{tot} = \sum_{i=1}^{N} n_i$ total entries and mean values $\mu_i$, or equivalently bin probabilities $p_i = \mu_i/n_{tot}$. Alternatively one may model the $n_i$ as independent and Poisson distributed with means $\mu_i$. If there are no adjustable parameters in the hypotheses, then the goodness-of-fit can be quantified with Pearson's $\chi^2$ statistic, defined for multinomial data as

$$t_M = \sum_{i=1}^{N} \frac{(n_i - n_{tot}p_i)^2}{n_{tot}p_i}, \qquad (40.56)$$

or for the independent Poisson data as

$$t_P = \sum_{i=1}^{N} \frac{(n_i - \mu_i)^2}{\mu_i}. \qquad (40.57)$$

In the limit where the means $\mu_i$ are large, the statistics $t_M$ and $t_P$ 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 $n_{tot}$ is fixed and thus the $n_i$ are correlated.

Alternatively one can use the test statistic based on the likelihood ratio $t = -2 \ln \lambda(\hat{\theta})$ given previously in Eq. (40.16). Here

the likelihood ratio $\lambda$ 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 the likelihood-ratio statistic approaches the asymptotic limit faster than does Pearson's chi-square and thus when using the chi-square p.d.f. to compute a $p$-value one obtains in general a more accurate result from $-2 \ln \lambda$ (see [9]).

Since the mean of the chi-square distribution is equal to $n_d$, one expects in a “reasonable” experiment to obtain $\chi^2 \approx n_d$ (here $\chi^2$ refers to the minimized value). Hence the quantity $\chi^2/n_d$ is sometimes reported. Since the p.d.f. of $\chi^2/n_d$ depends on $n_d$, however, one must report $n_d$ as well if one wishes to determine the $p$-value. The $p$-values obtained for different values of $\chi^2/n_d$ are shown in Fig. 40.2.

If the minimized $\chi^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-of-fit, however, does not mean that one will have large statistical errors for parameter estimates. If, for example, the error bars (or covariance matrix) used in constructing the $\chi^2$ are underestimated, then this will lead to underestimated statistical errors for the fitted parameters and an increased value of the minimized $\chi^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 many times, assuming that the hypothesis and measurement errors used in the $\chi^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 $\chi^2$.

<figure><img src="image_258.png" /><figcaption>**Figure 40.1:** One minus the $\chi^2$ cumulative distribution, $1 - F(\chi^2; n)$, for $n$ degrees of freedom. This gives the $p$-value for the $\chi^2$ goodness-of-fit test as well as one minus the coverage probability for confidence regions (see Sec. 40.4.2.2).</figcaption></figure>

<figure><img src="image_259.png" /><figcaption>**Figure 40.2:** The ‘reduced’ $\chi^2$, equal to $\chi^2/n$, for $n$ degrees of freedom. The curves show as a function of $n$ the $\chi^2/n$ that corresponds to a given $p$-value.</figcaption></figure>
### 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 = (x_1, \dots, x_n)$. 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 of each of the $x_i$ (i.e., “unbinned”). The sample may be displayed graphically using the empirical distribution function (e.d.f.)

$$F_n(x) = \frac{\text{number of elements in sample with value } \le x}{n}. \qquad (40.58)$$

The e.d.f. $F_n(x)$ begins at zero for $x$ less than the lowest observed value, increases by $1/n$ at each observed $x_i$, and saturates at unity for $x$ equal to the highest value in the sample.

Given an i.i.d. data sample $x = (x_1, \dots, x_n)$, suppose one wants to test the hypothesis $H_0$ that $x$ follows the p.d.f. $f(x)$ using a goodness-of-fit test, i.e., without reference to a specific alternative. The hypothesis of the p.d.f. $f(x)$ is equivalent to saying the cumulative distribution function (c.d.f.) is $F(x) = \int_{-\infty}^{x} f(x') dx'$. As the e.d.f. $F_n(x)$ can be viewed as an estimator for the c.d.f. $F(x)$, we can derive 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

$$D_n = \sup_x |F_n(x) - F(x)|, \qquad (40.59)$$

where $sup$ denotes the supremum. That is, $D_n$ is the greatest vertical distance between $F_n(x)$ and $F(x)$ for any $x$. In the limit of large sample size $n$, one can show that the distribution of $D_n$, 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 value of $D_n$ corresponds to greater incompatibility between data and hypothesis, so the $p$-value is the probability $P(D_n \ge D_{n,\text{obs}}|H_0)$. 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 $F_n(x)$ and $F(x)$, giving a large maximum vertical distance $D_n$.

The Cramér-von Mises family of tests is based on the squared difference between $F_n(x)$ and $F(x)$ through an integral of the form

$$W^2 = \int_{-\infty}^{\infty} (F(x) - F_n(x))^2 w(x) f(x) dx, \qquad (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 tails of the distribution. For cases such as these where the weight function depends on $x$ only through the 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 the Bayes 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), $H_i$ and $H_j$, described by vectors of parameters $\theta_i$ and $\theta_j$, 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

$$\pi(H_i, \theta_i) = P(H_i)\pi(\theta_i|H_i). \qquad (40.61)$$

Here $P(H_i)$ is the overall prior probability for $H_i$, and $\pi(\theta_i|H_i)$ is the normalized p.d.f. of its parameters. For each model, the posterior probability is found using Bayes' theorem,

$$P(H_i|x) = \frac{\int P(x|\theta_i, H_i)P(H_i)\pi(\theta_i|H_i) d\theta_i}{P(x)}, \qquad (40.62)$$

where the integration is carried out over the internal parameters $\theta_i$ of the model. The ratio of posterior probabilities for the models is therefore

$$\frac{P(H_i|x)}{P(H_j|x)} = \frac{\int P(x|\theta_i, H_i)\pi(\theta_i|H_i) d\theta_i}{\int P(x|\theta_j, H_j)\pi(\theta_j|H_j) d\theta_j} \frac{P(H_i)}{P(H_j)}. \qquad (40.63)$$

The Bayes factor is defined as

$$B_{ij} = \frac{\int P(x|\theta_i, H_i)\pi(\theta_i|H_i) d\theta_i}{\int P(x|\theta_j, H_j)\pi(\theta_j|H_j) d\theta_j}. \qquad (40.64)$$

This gives what the ratio of posterior probabilities for models $i$ and $j$ would be if the overall prior probabilities for the two models were equal. If the models have no nuisance parameters, i.e., no internal 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 the data in favor of one hypothesis over the other.

Although the Bayes factor is by construction independent of the overall prior probabilities $P(H_i)$ and $P(H_j)$, it does require priors for all internal parameters of a model, i.e., one needs the functions $\pi(\theta_i|H_i)$ and $\pi(\theta_j|H_j)$. In a Bayesian analysis where one is only interested in the posterior 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 so for 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 meaningful probabilities over the entire range where they are defined.

An exception to this rule may be considered when the identical parameter appears in the models for both numerator and denominator of the Bayes factor. In this case one can argue that the arbitrary constants would cancel. One must exercise some caution, however, as parameters with 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

$$m = \int P(x|\theta)\pi(\theta) d\theta, \qquad (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 be done 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].

## 40.4 Intervals and limits

When the goal of an experiment is to determine a parameter $\theta$, the result is usually expressed by quoting, in addition to the point estimate, some sort of interval which reflects the statistical precision of the measurement. In the simplest case, this can be
given by the parameter's estimated value $\hat{\theta}$ plus or minus an estimate of the standard deviation of $\hat{\theta}$, $\hat{\sigma}_{\hat{\theta}}$. If, however, the p.d.f. of the 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 parameter with a specified probability;

* provide the information needed by the consumer of the result to draw conclusions about the parameter 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 the multiparameter case, the parameter estimates and covariance matrix) satisfy essentially all of these goals. 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 practical considerations 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 to which 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 regions that 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) $[\theta_{lo}, \theta_{up}]$ can be determined which contains a given fraction $1 - \alpha$ of the posterior probability, i.e.,

$$1 - \alpha = \int_{\theta_{lo}}^{\theta_{up}} p(\theta|x) \, d\theta. \qquad (40.66)$$

Sometimes an upper or lower limit is desired, i.e., $\theta_{lo}$ or $\theta_{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 $\theta$ values for which $p(\theta|x)$ is higher than for any $\theta$ not belonging to the set, which may constitute a single interval or a set of disjoint regions; these are called highest posterior density (HPD) intervals. Note that HPD intervals 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 zero for 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

$$\pi(s) = \begin{cases} 0 & s < 0 \\ 1 & s \ge 0 \end{cases}. \qquad (40.67)$$

This prior may be regarded as providing an interval whose frequentist properties can be studied, rather than as representing a degree of belief. For example, to obtain an upper limit on $s$, one may proceed as follows. The likelihood for $s$ is given by the Poisson distribution for $n$ with mean $s + b$,

$$P(n|s) = \frac{(s+b)^n}{n!} e^{-(s+b)}, \qquad (40.68)$$

along with the prior (40.67) in (40.37) gives the posterior density for $s$. An upper limit $s_{up}$ at confidence level (or here, rather,

credibility level) $1 - \alpha$ can be obtained by requiring

$$1 - \alpha = \int_{-\infty}^{s_{up}} p(s|n) ds = \frac{\int_{-\infty}^{s_{up}} P(n|s) \pi(s) ds}{\int_{-\infty}^{\infty} P(n|s) \pi(s) ds}, \qquad (40.69)$$

where the lower limit of integration is effectively zero because of the cut-off in $\pi(s)$. By relating the integrals in Eq. (40.69) to incomplete gamma functions, the solution for the upper limit is found to be

$$s_{up} = \frac{1}{2} F_{\chi^2}^{-1} [p, 2(n+1)] - b, \qquad (40.70)$$

where $F_{\chi^2}^{-1}$ is the quantile of the $\chi^2$ distribution (inverse of the cumulative distribution). Here the quantity $p$ is

$$p = 1 - \alpha (1 - F_{\chi^2} [2b, 2(n+1)]), \qquad (40.71)$$

where $F_{\chi^2}$ is the cumulative $\chi^2$ distribution. For both $F_{\chi^2}$ and $F_{\chi^2}^{-1}$ above, the argument $2(n+1)$ gives the number of degrees of freedom. For the special case of $b = 0$, the limit reduces to

$$s_{up} = \frac{1}{2} F_{\chi^2}^{-1} (1 - \alpha; 2(n+1)). \qquad (40.72)$$

It happens that for the case of $b = 0$, the upper limit from Eq. (40.72) coincides numerically with the frequentist upper limit discussed in Section 40.4.2.3. Values for $1 - \alpha = 0.9$ and $0.95$ are given by the values $\mu_{up}$ in Table 40.3. The frequentist properties of confidence intervals for the Poisson mean 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 assumption of 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 $\pi(s) \propto 1/\sqrt{s+b}$ for $s \ge 0$ and zero otherwise. 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. $\pi(b)$, then one would first use Bayes' theorem to find $p(s, b|n)$. This is then marginalized to find $p(s|n) = \int p(s, b|n)\pi(b) db$, from which one may determine an interval for $s$. One may not be certain whether to extend a model by including more nuisance parameters. In this case, a Bayes factor may be used to determine to what 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 a procedure 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 were to repeat the experiment many times. The coverage probability refers to the fraction of intervals in such an ensemble that contain the true parameter value. Confidence intervals are constructed so as to have a coverage probability greater than or equal to a given confidence level, regardless of the true 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; \theta)$ where $x$ represents the outcome of the experiment and $\theta$ is the unknown parameter for which we want to construct a confidence interval. The variable $x$ could (and often does) represent an estimator for $\theta$. Using $f(x; \theta)$, we can find using a pre-defined rule and probability $1 - \alpha$ for every value of
$\theta$, a set of values $x_1(\theta, \alpha)$ and $x_2(\theta, \alpha)$ such that

$$P(x_1 < x < x_2; \theta) = \int_{x_1}^{x_2} f(x; \theta) dx \geq 1 - \alpha. \quad (40.73)$$

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 - \alpha$, and one requires by convention $P(x_1 < x < x_2; \theta) \geq 1 - \alpha$.

This is illustrated for continuous $x$ in Fig. 40.3: a horizontal line segment $[x_1(\theta, \alpha), x_2(\theta, \alpha)]$ is drawn for representative values of $\theta$. The union of such intervals for all values of $\theta$, designated in the figure as $D(\alpha)$, is known as a *confidence belt*. Typically the curves $x_1(\theta, \alpha)$ and $x_2(\theta, \alpha)$ are monotonic functions of $\theta$, which we assume for this discussion.

<figure><img src="image_260.png" /><figcaption>Figure 40.3: Construction of the confidence belt (see text).</figcaption></figure>

Upon performing an experiment to measure $x$ and obtaining a value $x_0$, one draws a vertical line through $x_0$. The confidence interval for $\theta$ is the set of all values of $\theta$ for which the corresponding line segment $[x_1(\theta, \alpha), x_2(\theta, \alpha)]$ is intercepted by this vertical line. Such confidence intervals are said to have a *confidence level* (CL) equal to $1 - \alpha$.

Now suppose that the true value of $\theta$ is $\theta_0$, indicated in the figure. We see from the figure that $\theta_0$ lies between $\theta_1(x)$ and $\theta_2(x)$ if and only if $x$ lies between $x_1(\theta_0)$ and $x_2(\theta_0)$. The two events thus have the same probability, and since this is true for any value $\theta_0$, we can drop the subscript 0 and obtain

$$1 - \alpha = P(x_1(\theta) < x < x_2(\theta)) = P(\theta_2(x) < \theta < \theta_1(x)). \quad (40.74)$$

In this probability statement, $\theta_1(x)$ and $\theta_2(x)$, i.e., the endpoints of the interval, are the random variables and $\theta$ is an unknown constant. If the experiment were to be repeated a large number of times, the interval $[\theta_1, \theta_2]$ would vary, covering the fixed value $\theta$ in a fraction $1 - \alpha$ of the experiments.

The condition of coverage in Eq. (40.73) does not determine $x_1$ and $x_2$ uniquely, and additional criteria are needed. One possibility is to choose *central intervals* such that the probabilities to find $x$ below $x_1$ and above $x_2$ are each $\alpha/2$. In other cases, one may want to report only an upper or lower limit, in which case one of $P(x \le x_1)$ or $P(x \ge x_2)$ can be set to $\alpha$ 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 - \alpha$, regardless of the true value of the parameter. Because of the requirement $P(x_1 < x < x_2) \ge 1 - \alpha$ 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 - \alpha$.

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 $\theta$ (assume one constructs a test for all physical

values of $\theta$). One then excludes all values of $\theta$ where the hypothesis would be rejected in a test of size $\alpha$ or less. The remaining values constitute the confidence interval at confidence level $1 - \alpha$. If the critical region of the test is characterized by having a $p$-value $p_\theta \le \alpha$, then the endpoints of the confidence interval are found in practice by solving $p_\theta = \alpha$ for $\theta$.

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 the confidence belt. One possibility is to use a test statistic based on the *likelihood ratio*,

$$\lambda(\theta) = \frac{f(x; \theta)}{f(x; \hat{\theta})}, \quad (40.75)$$

where $\hat{\theta}$ is the value of the parameter which, out of all allowed values, maximizes $f(x; \theta)$. The confidence belt is taken to contain the values of $x$ that give the greatest values of $\lambda(\theta)$. This results in the intervals described in Ref. [42] by Feldman and Cousins.

If the model contains nuisance parameters $\nu$, then these can be incorporated into the test (or the $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 $\theta$ as excluded only if it is rejected for all possible values of $\nu$. The resulting interval for $\theta$ will then cover the true value with a probability greater than or equal to the nominal confidence level for all points in $\nu$-space.

If the $p$-value is based on the profiled values of the nuisance parameters, i.e., with $\nu = \hat{\nu}(\theta)$ used in Eq. (40.48), then the resulting interval for the parameter of interest will have the correct coverage if the true values of $\nu$ 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] (see also Ref. [26]).

#### 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 an estimator 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 $\sigma$, the probability that the measured value $x$ will fall within $\pm\delta$ of the true value $\mu$ is

$$\begin{aligned} 1 - \alpha &= \frac{1}{\sqrt{2\pi}\sigma} \int_{\mu - \delta}^{\mu + \delta} e^{-(x - \mu)^2 / 2\sigma^2} dx \\ &= \text{erf} \left( \frac{\delta}{\sqrt{2} \sigma} \right) = 2\Phi \left( \frac{\delta}{\sigma} \right) - 1, \end{aligned} \quad (40.76)$$

where erf is the Gaussian error function, which is rewritten in the final equality using $\Phi$, the Gaussian cumulative distribution. Fig. 40.4 shows a $\delta = 1.64\sigma$ confidence interval unshaded. The choice $\delta = \sigma$ gives an interval called the *standard error* which has $1 - \alpha = 68.27\%$ if $\sigma$ is known. Values of $\alpha$ for other frequently used choices of $\delta$ are given in Table 40.1.

<table><caption>**Table 40.1:** Area of the tails $\alpha$ outside $\pm\delta$ from the mean of a Gaussian distribution.</caption><thead><tr><th>α</th><th>δ</th><th>α</th><th>δ</th></tr></thead><tbody><tr><td>0.3173</td><td>1σ</td><td>0.2</td><td>1.28σ</td></tr><tr><td>4.55 × 10<sup>-2</sup></td><td>2σ</td><td>0.1</td><td>1.64σ</td></tr><tr><td>2.7 × 10<sup>-3</sup></td><td>3σ</td><td>0.05</td><td>1.96σ</td></tr><tr><td>6.3 × 10<sup>-5</sup></td><td>4σ</td><td>0.01</td><td>2.58σ</td></tr><tr><td>5.7 × 10<sup>-7</sup></td><td>5σ</td><td>0.001</td><td>3.29σ</td></tr><tr><td>2.0 × 10<sup>-9</sup></td><td>6σ</td><td>10<sup>-4</sup></td><td>3.89σ</td></tr></tbody></table>

We can set a one-sided (upper or lower) limit by excluding above $x + \delta$ (or below $x - \delta$). The values of $\alpha$ for such limits are half the values in Table 40.1.
<figure><img src="image_261.png" /><figcaption>**Figure 40.4:** Illustration of a symmetric 90% confidence interval (unshaded) for a Gaussian-distributed measurement of a single quantity. Integrated probabilities, defined by $\alpha = 0.1$, are as shown.</figcaption></figure>

The relation (40.76) can be re-expressed using the cumulative distribution function for the $\chi^2$ distribution as

$$ \alpha = 1 - F(\chi^2; n), \qquad (40.77) $$

for $\chi^2 = (\delta/\sigma)^2$ and $n = 1$ degree of freedom. This can be seen as the $n = 1$ curve in Fig. 40.1 or obtained using standard computer libraries. For multivariate measurements of, say, $M$ parameter estimates $\hat{\theta} = (\hat{\theta}_1, \dots, \hat{\theta}_M)$, construction of the confidence region requires the full covariance matrix $V_{ij} = \text{cov}[\hat{\theta}_i, \hat{\theta}_j]$, 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 the true (unknown) values $\theta$, and furthermore, the likelihood function itself will take on a Gaussian shape.

The standard error ellipse for the pair $(\hat{\theta}_i, \hat{\theta}_j)$ is shown in Fig. 40.5, corresponding to a contour $\chi^2 = \chi_{\min}^2 + 1$ or $\ln L = \ln L_{\max} - 1/2$. The ellipse is centered about the estimated values $\hat{\theta}$, and the tangents to the ellipse give the standard deviations of the estimators, $\sigma_i$ and $\sigma_j$. The angle of the major axis of the ellipse is given by

$$ \tan 2\phi = \frac{2\rho_{ij}\sigma_i\sigma_j}{\sigma_j^2 - \sigma_i^2}, \qquad (40.78) $$

where $\rho_{ij} = \text{cov}[\hat{\theta}_i, \hat{\theta}_j]/\sigma_i\sigma_j$ is the correlation coefficient.

The correlation coefficient can be visualized as the fraction of the distance $\sigma_i$ from the ellipse's horizontal center-line at which the ellipse becomes tangent to vertical, i.e., at the distance $\rho_{ij}\sigma_i$ below the center-line as shown. As $\rho_{ij}$ goes to $+1$ or $-1$, the ellipse thins to a diagonal line.

It could happen that one of the parameters, say, $\theta_j$, is known from previous measurements to a precision much better than $\sigma_j$, so that the current measurement contributes almost nothing to the knowledge of $\theta_j$. However, the current measurement of $\theta_i$ and its dependence on $\theta_j$ may still be important. In this case, instead of quoting both parameter estimates and their correlation, one sometimes reports the value of $\theta_i$, which minimizes $\chi^2$ at a fixed value of $\theta_j$, such as the PDG best value. This $\theta_i$ value lies along the dotted line between the points where the ellipse becomes tangent to vertical, and has statistical error $\sigma_{\text{inner}}$ as shown on the figure, where $\sigma_{\text{inner}} = (1 - \rho_{ij}^2)^{1/2}\sigma_i$. Instead of the correlation $\rho_{ij}$, one reports the dependency $d\hat{\theta}_i/d\theta_j$, which is the slope of the dotted line. This slope is related to the correlation coefficient by $d\hat{\theta}_i/d\theta_j = \rho_{ij} \times \frac{\sigma_i}{\sigma_j}$.

As in the single-variable case, because of the symmetry of the Gaussian function between $\theta$ and $\hat{\theta}$, one finds that contours of

<figure><img src="image_262.png" /><figcaption>**Figure 40.5:** Standard error ellipse for the estimators $\hat{\theta}_i$ and $\hat{\theta}_j$. In the case shown the correlation is negative.</figcaption></figure>

constant $\ln L$ or $\chi^2$ cover the true values with a certain, fixed probability. That is, the confidence region is determined by

$$ \ln L(\theta) \ge \ln L_{\max} - \Delta \ln L, \qquad (40.79) $$

or where a $\chi^2$ has been defined for use with the method of least-squares,

$$ \chi^2(\theta) \le \chi_{\min}^2 + \Delta\chi^2. \qquad (40.80) $$

Values of $\Delta\chi^2$ or $2\Delta\ln L$ are given in Table 40.2 for several values of the coverage probability $1 - \alpha$ and number of fitted parameters $M$. For Gaussian distributed data, these are related by $\Delta\chi^2 = 2\Delta\ln L = F_{\chi_M^2}^{-1}(1 - \alpha)$, where $F_{\chi_M^2}^{-1}$ is the chi-square quantile (inverse of the cumulative distribution) for $M$ degrees of freedom.

<table><caption>**Table 40.2:** Values of $\Delta\chi^2$ or $2\Delta\ln L$ corresponding to a coverage probability $1 - \alpha$ in the large data sample limit, for joint estimation of $M$ parameters.</caption><thead><tr><th>(1 - α) (%)</th><th>M = 1</th><th>M = 2</th><th>M = 3</th></tr></thead><tbody><tr><td>68.27</td><td>1.00</td><td>2.30</td><td>3.53</td></tr><tr><td>90.</td><td>2.71</td><td>4.61</td><td>6.25</td></tr><tr><td>95.</td><td>3.84</td><td>5.99</td><td>7.82</td></tr><tr><td>95.45</td><td>4.00</td><td>6.18</td><td>8.03</td></tr><tr><td>99.</td><td>6.63</td><td>9.21</td><td>11.34</td></tr><tr><td>99.73</td><td>9.00</td><td>11.83</td><td>14.16</td></tr></tbody></table>

For non-Gaussian data samples, the probability for the regions determined by Equations (40.79) or (40.80) to cover the true value of $\theta$ becomes independent of $\theta$ 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 the true 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 $\chi^2$ is that it is relatively simple to produce the interval for the combination of several experiments. If $N$ independent measurements result in log-likelihood functions $\ln L_i(\theta)$, then the combined log-likelihood function is simply the sum,

$$ \ln L(\theta) = \sum_{i=1}^{N} \ln L_i(\theta). \qquad (40.81) $$

This 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 $\sigma$ and selection efficiency $\epsilon$ in a fixed integrated luminosity $\mathcal{L}$, then it follows a Poisson distribution with mean $\mu = \sigma\epsilon\mathcal{L}$. If, on the other hand, one has selected a larger sample of $N$ events and found $n$ of them 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 $\mu$ 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 and lower limits are found to be

$$ \mu_{lo} = \frac{1}{2} F_{\chi^2}^{-1}(\alpha_{lo}; 2n), \quad (40.82a) $$

$$ \mu_{up} = \frac{1}{2} F_{\chi^2}^{-1}(1 - \alpha_{up}; 2(n + 1)), \quad (40.82b) $$

where confidence levels of $1 - \alpha_{lo}$ and $1 - \alpha_{up}$ refer separately to the corresponding intervals $\mu \ge \mu_{lo}$ and $\mu \le \mu_{up}$, and $F_{\chi^2}^{-1}$ is the quantile of the $\chi^2$ distribution (inverse of the cumulative distribution). For central confidence intervals at confidence level $1 - \alpha$, set $\alpha_{lo} = \alpha_{up} = \alpha/2$.

<table><caption>**Table 40.3:** Lower and upper (one-sided) limits for the mean $\mu$ of a Poisson variable given $n$ observed events in the absence of background, for confidence levels of 90% and 95%.</caption><thead><tr><th></th><th colspan="2">1 - α =90%</th><th colspan="2">1 - α =95%</th></tr><tr><th>n</th><th>μ<sub>lo</sub></th><th>μ<sub>up</sub></th><th>μ<sub>lo</sub></th><th>μ<sub>up</sub></th></tr></thead><tbody><tr><td>0</td><td>-</td><td>2.30</td><td>-</td><td>3.00</td></tr><tr><td>1</td><td>0.105</td><td>3.89</td><td>0.051</td><td>4.74</td></tr><tr><td>2</td><td>0.532</td><td>5.32</td><td>0.355</td><td>6.30</td></tr><tr><td>3</td><td>1.10</td><td>6.68</td><td>0.818</td><td>7.75</td></tr><tr><td>4</td><td>1.74</td><td>7.99</td><td>1.37</td><td>9.15</td></tr><tr><td>5</td><td>2.43</td><td>9.27</td><td>1.97</td><td>10.51</td></tr><tr><td>6</td><td>3.15</td><td>10.53</td><td>2.61</td><td>11.84</td></tr><tr><td>7</td><td>3.89</td><td>11.77</td><td>3.29</td><td>13.15</td></tr><tr><td>8</td><td>4.66</td><td>12.99</td><td>3.98</td><td>14.43</td></tr><tr><td>9</td><td>5.43</td><td>14.21</td><td>4.70</td><td>15.71</td></tr><tr><td>10</td><td>6.22</td><td>15.41</td><td>5.43</td><td>16.96</td></tr></tbody></table>

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 $\mu$. Values for confidence levels of 90% and 95% are shown in Table 40.3. For the case of binomially distributed $n$ successes out of $N$ trials with probability of success $p$, the upper and lower limits on $p$ are found to be

$$ p_{lo} = \frac{n F_F^{-1}[\alpha_{lo}; 2n, 2(N - n + 1)]}{N - n + 1 + n F_F^{-1}[\alpha_{lo}; 2n, 2(N - n + 1)]}, \quad (40.83a) $$

$$ p_{up} = \frac{(n+1)F_F^{-1}[1 - \alpha_{up}; 2(n+1), 2(N-n)]}{(N-n) + (n+1)F_F^{-1}[1 - \alpha_{up}; 2(n+1), 2(N-n)]}. \quad (40.83b) $$

Here $F_F^{-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 $\mu$ is defined such that it is proportional to the signal cross section. A statistical test may be carried out for hypothesized values of $\mu$, which may be done by computing a $p$-value, $p_{\mu}$, for all $\mu$. Those values not rejected in a test of size $\alpha$, i.e., for which one does not find $p_{\mu} \le \alpha$, constitute a confidence interval with confidence level $1 - \alpha$.

In general one will find that for some regions in the parameter space of the signal model, the predictions for data are almost indistinguishable from those of the background-only model. This corresponds to the case where $\mu$ 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) is very 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 $\alpha$, the type-I error rate. As one often takes $\alpha$ 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 the experimental 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 the resulting upper limit on $\mu$ 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 CL$_s$ method, where the measure used to test a parameter is increased for decreasing sensitivity [43, 44]. The procedure is based on a statistic called CL$_s$, which is defined as

$$ \mathrm{CL}_s = \frac{p_\mu}{1 - p_b}, \quad (40.84) $$

where $p_b$ is the $p$-value of the background-only hypothesis. In the usual formulation of the method, both $p_\mu$ and $p_b$ are defined using a single test statistic, and the definition of CL$_s$ 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 CL$_s \le \alpha$. As the denominator in Eq. (40.84) is always less than or equal to unity, the exclusion criterion based on CL$_s$ is more stringent than the usual requirement $p_\mu \le \alpha$. In this sense the CL$_s$ procedure is conservative, and the coverage probability of the corresponding intervals will exceed the nominal confidence level $1 - \alpha$. If the experimental sensitivity to a given value of $\mu$ is very low, then one finds that as $p_\mu$ decreases, so does the denominator $1 - p_b$, and thus the condition CL$_s \le \alpha$ is effectively prevented from being satisfied. In this way the exclusion of parameters in the case of low sensitivity is suppressed.

The CL$_s$ procedure has the attractive feature that the resulting intervals coincide with those obtained from the Bayesian method in two important cases: the mean value of a Poisson or Gaussian distributed measurement with a constant prior. The CL$_s$ intervals overcover for all values of the parameter $\mu$, however, by an amount that depends on $\mu$.

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 $\mu$ against the alternative of a lower rate, and therefore the critical region 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 physically meaningful signal parameter values, including those to which one has very little sensitivity, are rejected 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 $\mu$ 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 $\hat{\mu}$ is sufficiently close to the lower limit of zero, then only high values of $\mu$ are rejected, and the lower edge of the confidence interval is at zero. Note, however, that the coverage property of $1 - \alpha$ pertains to the entire interval, not to the probability for the upper edge $\mu_{up}$ to be greater than the true value $\mu$. For parameter estimates increasingly far away from the boundary, i.e., for increasing signal significance, the point $\mu = 0$ is excluded and the interval has nonzero upper and lower edges.
An additional difficulty arises when a parameter estimate is not significantly far away from the boundary, 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 $x_2 = \infty$ in Eq. (40.73). Then $x_1$ 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 with equal upper and lower tail areas). As pointed out by Feldman and Cousins [42], if the decision to report an upper limit or two-sided interval is made by looking at the data (“flip-flopping”), then in general there will be parameter values for which the resulting intervals have a coverage probability less than $1 - \alpha$. 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 (and are 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 $\alpha$ 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 - \alpha$.) The values of $1 - \alpha$ given here refer to the coverage of the true parameter by the whole interval $[\mu_1, \mu_2]$. In Table 40.3 for the one-sided upper limit, however, $1 - \alpha$ refers to the probability to have $\mu_{up} \ge \mu$ (or $\mu_{lo} \le \mu$ for lower limits).

<table><caption>**Table 40.4:** Unified confidence intervals $[\mu_1, \mu_2]$ for a the mean of a Poisson variable given $n$ observed events in the absence of background, for confidence levels of 90% and 95%.</caption><thead><tr><th rowspan="2">n</th><th colspan="2">1 - α =90%</th><th colspan="2">1 - α =95%</th></tr><tr><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th></tr></thead><tbody><tr><td>0</td><td>0.00</td><td>2.44</td><td>0.00</td><td>3.09</td></tr><tr><td>1</td><td>0.11</td><td>4.36</td><td>0.05</td><td>5.14</td></tr><tr><td>2</td><td>0.53</td><td>5.91</td><td>0.36</td><td>6.72</td></tr><tr><td>3</td><td>1.10</td><td>7.42</td><td>0.82</td><td>8.25</td></tr><tr><td>4</td><td>1.47</td><td>8.60</td><td>1.37</td><td>9.76</td></tr><tr><td>5</td><td>1.84</td><td>9.99</td><td>1.84</td><td>11.26</td></tr><tr><td>6</td><td>2.21</td><td>11.47</td><td>2.21</td><td>12.75</td></tr><tr><td>7</td><td>3.56</td><td>12.53</td><td>2.58</td><td>13.81</td></tr><tr><td>8</td><td>3.96</td><td>13.99</td><td>2.94</td><td>15.29</td></tr><tr><td>9</td><td>4.36</td><td>15.30</td><td>4.36</td><td>16.77</td></tr><tr><td>10</td><td>5.50</td><td>16.50</td><td>4.75</td><td>17.82</td></tr></tbody></table>

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 - \alpha = 0.9$. If the true signal parameter is zero, or in any case much less than the expected background, one will usually obtain a one-sided upper limit on $s$. In a certain fraction of the experiments, however, a two-sided interval for $s$ will result. Since, however, one typically chooses $1 - \alpha$ to be only 0.9 or 0.95 when setting limits, the value $s = 0$ may be found below the lower edge of the interval before the existence of the effect is well established. It must then be communicated carefully that in excluding $s = 0$ at, say, 90% or 95% confidence level from the interval, one is not necessarily claiming to have discovered the effect, for which one would usually require a higher level of significance (e.g., 5 $\sigma$).

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 the unphysical 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 $\mu$ of a Poisson distributed measurement has been the uniform distribution for $\mu \ge 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 $\mu$ 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 Poisson parameter, 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 the parameter, 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 new physics signal. As discussed in Sec. 40.3, an experimental analysis can often be formulated as a test of hypothetical model parameters. Therefore we may quantify the sensitivity by giving the results that we expect from such a test under specific assumptions about the signal process.

Here to be concrete we will consider a parameter $\mu$ 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 $\mu = 0$, and in addition one often wants to test nonzero values of $\mu$ to construct a confidence interval (e.g., limits) as described in Sec. 40.4. In the frequentist framework, the result of each tested value of $\mu$ is the p-value $p_\mu$ or equivalently the significance $Z_\mu = \Phi^{-1}(1 - p_\mu)$, where as usual $\Phi$ is the standard Gaussian cumulative distribution and its inverse $\Phi^{-1}$ is the standard Gaussian quantile.

Prior to carrying out the experiment, one generally wants to quantify what significance $Z_\mu$ is expected under given assumptions for the presence or absence of the signal process. Specifically, for the significance of a test of $\mu = 0$ (the discovery significance) one usually quotes the $Z_0$ one would expect if the signal is present at a given nominal rate, which we can define in general to correspond to $\mu = 1$. For limits, one often gives the expected limit under assumption of the background-only ($\mu = 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, and to say what values of $\mu$ 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_\mu$ 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 $\mu$, which is not necessarily the same as the hypothesized value of $\mu$ being tested. We may therefore consider the distribution $f(Z_\mu|\mu')$, i.e., the distribution of $Z_\mu$ that would be obtained by considering data samples generated under assumption of $\mu'$. In a similar way one can talk about the sampling distribution of an upper limit for $\mu$, $f(\mu_{up}|\mu')$.

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 - \Phi(Z)$, then gives $\text{med}[p_\mu|\mu'] = 1 - \Phi(\text{med}[Z_\mu|\mu'])$, whereas the corresponding relation does not hold in general for the mean.

In some cases one may be able to write down approximate formulae for the distributions of $Z_\mu$ and for limits, but more generally they must be determined from Monte Carlo calculations. In many cases of interest, the significance $Z_\mu$ and the limits on $\mu$ 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 $\mu s + b$. Here $s$ and $b$,
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 $\mu = 0$. To characterize the experimental sensitivity, we want to give the discovery significance expected under the assumption of $\mu = 1$.

In the limit where its mean value is large, the Poisson variable $n$ can be approximated as an almost continuous Gaussian variable with mean $\mu s + b$ and standard deviation $\sigma = \sqrt{\mu s + b}$. In the usual case where a physical signal model corresponds to $\mu > 0$, the p-value of $\mu = 0$ is the probability to find $n$ greater than or equal to the value observed,

$$p_0 = \Phi \left( \frac{n-b}{\sqrt{b}} \right), \qquad (40.85)$$

and the corresponding significance is $Z_0 = \Phi^{-1}(1 - p_0) = (n - b)/\sqrt{b}$. The median (here equal to the mean) of $n$ assuming $\mu = 1$ is $s + b$, and therefore the median discovery significance is

$$\text{med}[Z_0 | \mu = 1] = \frac{s}{\sqrt{b}}. \qquad (40.86)$$

The figure of merit "$s/\sqrt{b}$" has been widely used in particle physics as a measure of expected discovery significance. A better approximation for the Poisson counting experiment, however, may be obtained by testing $\mu = 0$ using the likelihood ratio (40.49) $\lambda(0) = L(0)/L(\hat{\mu})$, where

$$L(\mu) = \frac{(\mu s + b)^n}{n!} e^{-(\mu s + b)} \qquad (40.87)$$

is the likelihood function, and $\hat{\mu} = (n - b)/s$ is the estimator of $\mu$. In this example there are no nuisance parameters, as $s$ and $b$ are taken to be known. For the case where the relevant signal models correspond to positive $\mu$, one may test the $\mu = 0$ hypothesis with the statistic $q_0 = -2 \ln \lambda(0)$ when $\hat{\mu} > 0$, i.e., an excess is observed, and $q_0 = 0$ otherwise. One can show (see, e.g., Ref. [25]) that in the large-sample limit, the discovery significance is then $Z_0 = \sqrt{q_0}$, for which one finds

$$Z_0 = \sqrt{2 \left( n \ln \frac{n}{b} + b - n \right)} \qquad (40.88)$$

for $n > b$ and $Z_0 = 0$ otherwise. To approximate the expected discovery significance assuming $\mu = 1$, one may simply replace $n$ with the expected value $E[n|\mu = 1] = s+b$ (the so-called "Asimov data set"), giving

$$\text{med}[Z_0 | \mu = 1] = \sqrt{2 \left( (s+b) \ln \left(1 + \frac{s}{b}\right) - s \right)}. \qquad (40.89)$$

This has been shown in Ref. [25] to provide a good approximation to the median discovery significance for values of $s$ of several and for $b$ well below unity. The right-hand side of Eq. (40.89) reduces to $s/\sqrt{b}$ in the limit $s \ll b$.

Beyond the simple Poisson counting experiment, in general one may test values of a parameter $\mu$ with more complicated functions of the measured data to obtain a p-value $p_\mu$, and from this one can quote the equivalent significance $Z_\mu$ or find, e.g., an upper limit $\mu_{up}$. In this case as well one may quantify the experimental sensitivity by giving the significance $Z_\mu$ expected if the data are generated with a different value of the parameter $\mu'$. 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 method of 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 $Z_0$ obtained from the real data is subject to statistical fluctuations and will not in general be equal to its 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 is compared to expected statistical fluctuations. Therefore, in addition to the observed significance and limits it is useful to communicate not only their expected values but also a measure of the width 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 give the same probability content, i.e., $[0.1587, 0.8413]$ for $\pm 1\sigma$, $[0.02275, 0.97725]$ for $\pm 2\sigma$. 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 therefore find, usually with approximations or Monte Carlo studies, their sampling distribution and corresponding 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 with the likelihood from Eq. (40.87), the log of the Bayes factor that compares $\mu = 1$ to $\mu = 0$ is $\ln B_{10} = \ln(L(1)/L(0)) = n \ln(1 + s/b) - s$. That is, the expectation value, assuming $\mu = 1$, of $\ln B_{10}$ for this problem is

$$E[\ln B_{10} | \mu = 1] = (s+b) \ln \left(1 + \frac{s}{b}\right) - s. \qquad (40.90)$$

Comparing this to Eq. (40.89), one finds $\text{med}[Z_0|1] = \sqrt{2E[\ln B_{10}|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 $\hat{\mu} = (n-b)/s$ has a variance, assuming $\mu = 1$, of

$$V[\hat{\mu}] = V \left[ \frac{n-b}{s} \right] = \frac{1}{s^2} V[n] = \frac{s+b}{s^2}, \qquad (40.91)$$

so that the standard deviation of $\hat{\mu}$ is $\sigma_{\hat{\mu}} = \sqrt{s+b}/s$. One may therefore use $s/\sqrt{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/\sqrt{s+b}$ is also the expected significance with which one rejects $s$ assuming the signal is absent, 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: Classical Inference 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].

[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. Statist. Assoc.* **91**, 1343 (1996).

[19] J.M. Bernardo, J. R. Statist. Soc. **B41**, 113 (1979); J.M. Bernardo and J.O. Berger, *J. Am. Statist. 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*, CERN-2008-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 and Physical 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.* **A595**, 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).

[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 of x. Eq. 21.6 of Ref. [3] gives the rule: all such points are included in the acceptance region (the domain D(α)). Some authors have assumed the contrary, and shown that one can then obtain null 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].
# 41. Machine Learning

Written November 2021 by K. Cranmer (Wisconsin U.), U. Seljak (UC Berkeley; LBNL) and K. Terao (SLAC; Stanford U.).

41.1 Introduction 711

41.1.1 A gentle introduction with a representative example 712

41.2 Fundamental concepts 712

41.2.1 Loss, risk, empirical risk 712

41.2.2 Generalization 713

41.3 Common tasks and their associated loss functions 713

41.3.1 Supervised learning 713

41.3.1.1 Regression 713

41.3.1.2 A note on regularization 714

41.3.1.3 Classification 714

41.3.2 Unsupervised learning 715

41.3.2.1 Density estimation 715

41.3.2.2 Representation learning, compression, and auto-encoders 716

41.3.2.3 Clustering 716

41.3.3 Optimal control, reinforcement learning, and active learning 717

Reinforcement learning 717

Multi-arm bandits 717

Bayesian optimization 717

Connection to experimental design 718

Active learning 718

41.3.4 Anomaly detection and out-of-distribution detection 718

41.3.5 Simulation-based inference 719

41.3.5.1 Differentiable simulations 719

41.3.5.2 Unfolding as an inverse problem 719

41.4 Data representations, inductive bias, and example applications 720

41.5 Flavors of ML models 721

41.5.1 Support vector machines and kernel machines 721

Maximum-margin classifiers 721

Soft margins and slack variables 721

The dual problem 721

The kernel trick 721

Support vector regression 722

Kernel ridge regression 722

Gaussian Process Regression (krigging) 722

41.5.2 Decision trees 722

Tree-based models 722

Ensemble methods 723

Bagging 723

Random forests 723

AdaBoost 724

Gradient boosting 724

41.5.3 Neural networks 724

41.5.3.1 Feed-forward multi-layer perceptron 724

41.5.3.2 Activation functions 724

41.5.3.3 Softmax 725

41.5.3.4 The rise of deep learning 725

41.5.3.5 Convolutional neural networks 725

41.5.3.6 Pooling 725

41.5.3.7 CNN architectures for image analysis 726

Region Convolutional Neural Network 726

U-Net 726

41.5.3.8 Residual networks and skip connections 727

41.5.3.9 Recurrent neural networks 727

41.5.3.10 LSTM and GRU 729

41.5.3.11 Attention 729

41.5.3.12 Scaled dot-product attention 730

41.5.3.13 Transformer and multi-head attention 731

41.5.3.14 Graph networks and geometric deep learning 731

41.5.4 Deep generative models 732

41.5.4.1 Variational auto-encoders 733

41.5.4.2 Generative adversarial networks 733

41.5.4.3 Normalizing flows, autoregressive models, and score based models 734

41.6 Learning algorithms 735

41.6.1 Gradient-based optimization 735

41.6.2 Stochastic gradient descent 735

41.6.3 Optimization algorithms 735

41.6.4 Automatic differentiation and back propagation 735

41.6.5 The vanishing and exploding gradient problems 736

41.6.6 Early stopping 736

41.6.7 Initialization of model parameters 736

41.6.8 Input normalization 736

41.6.9 Batch normalization 736

41.6.10 Transfer learning: pre-training and fine-tuning 737

41.6.11 Zero, one, and a few shot learning 737

41.7 Incorporating uncertainty 737

41.7.1 Propagation of errors 738

41.7.2 Domain adaptation 738

41.7.3 Parameterized models 739

41.7.4 Data augmentation 739

41.7.5 Aleatoric and epistemic uncertainty 739

41.7.6 Model averaging and Bayesian machine learning 740

41.7.7 Connection to probabilistic machine learning 740

41.8 Infrastructure for deployment in experiments 741

## 41.1 Introduction

This chapter gives an overview of the core concepts of machine learning that are relevant to particle physics with some examples of applications to the energy, intensity, cosmic, and accelerator frontiers. Machine learning (ML) is an enormous field that has grown substantially in the last decade, propelled largely by the emergence of so-called deep learning (DL) [1, 2]. ML has a long history in particle physics going back to the late 1980s and early 1990s, see Refs. [3–5] for recent reviews.

Physicists are exploring and contributing to machine learning at an unprecedented rate, which poses a challenge for those that wish to have an up-to-date view of the field. This motivated an effort to create *A Living Review of Machine Learning for Particle and Nuclear Physics* [6], which can be downloaded here: https://github.com/iml-wg/HEPML-LivingReview. As of the time of this writing, the Living Review includes over 500 references organized hierarchically by topic. While we make references to some of these papers, this chapter focuses on the methodology and does not attempt to give a comprehensive review of the applications.

Despite the connotations of machine learning and artificial intelligence as a mysterious and radical departure from traditional approaches, we stress that machine learning has a mathematical formulation that is closely tied to statistics, the calculus of variations, approximation theory, and optimal control theory.

The topic can be organized along a few axes, which we use to organize this section. First, there are different learning paradigms, for example supervised learning, unsupervised learning, and reinforcement learning. We focus on the first two in this review, since reinforcement learning is less commonly used in particle physics. Within these paradigms there are various tasks; for example, classification and regression – which have been the primary use of ML in particle physics – are examples of supervised learning. In addition to the learning paradigm and tasks, there are various types of machine learning models that generically process some input and produce some output. The types of models vary based on what it is they are modelling (e.g. so-called discriminative vs. generative
models), as well as the way that they are implemented (e.g. neu-
ral networks, decision trees, and kernel machines). Next, there
are the issues around training or learning within the context of a
given task and model class, which connects to optimization and
regularization. We will briefly discuss the various considerations
that emerge in the application of machine learning methods to
physics, such as the treatment of systematic uncertainty, the in-
terpretability of the models, the incorporation of symmetry, etc..

**41.1.1 A gentle introduction with a representative example**

41.1.1 A gentle introduction with a representative example

We will use a specific, familiar example to introduce the vari-
ous ingredients in context before factorizing and abstracting them.
Consider the task of classifying energy deposits in a particle de-
tector as electrons or protons. For this example, let the detector
data consist of energy deposits in d sensors so that the data can
be represented as feature vector x ∈ R^d. Different components
of x may have different units (e.g. units of energy, momentum,
position etc.). Due to the complex interactions of particles in the
detector, we do not have an explicit probability model for the
high-dimensional data for the electron and proton scenarios, but
we do have a simulator that allows us to generate Monte Carlo
samples for each. This allows us to assemble a training dataset
{xᵢ, yᵢ}ᵢ=₁,...,n, where y is used as a label to identify how the
example was generated (e.g. y = 0 for electrons and y = 1 for
protons). We would like to find a function f: x → y that is able
to accurately predict the label on new data. Because we have
feature-label pairs, this is considered a supervised learning prob-
lem. We can use a neural network to provide a flexible family of
functions f<sub>φ</sub>: R<sup>d</sup> → R, where φ denotes the internal parameters
of the neural network (i.e. the weights and biases that we will dis-
cuss in Sec.41.5.3). The goal of the training procedure is to find
the value of the parameters φ that provide the 'best' predictions,
but since no model is perfect, we must be explicit about the trade-
offs. This is made concrete through a loss function L(f<sub>φ</sub>(x), y).
For this example, instead of the obvious zero-one loss (which is
0 if f<sub>φ</sub>(x) equals y and 1 if they are not), we choose to use the
squared-loss L<sub>sq</sub>(f<sub>φ</sub>(x), y) = (y - f<sub>φ</sub>(x))<sup>2</sup> (which may seem ad
hoc now, but will be motivated in Sec. 41.3.1.3). We can evaluate
the average of the loss on the training set of size n, which is re-
ferred to as the empirical risk R<sub>emp</sub>(f<sub>φ</sub>) = Σ<sub>i=1</sub><sup>n</sup> L(f<sub>φ</sub>(x), y)/n.
Training refers to numerically minimizing the empirical risk (often
referred to as the training loss through some abuse of terminol-
ogy). We can numerically optimize the model through gradient
descent, which iteratively adjusts the parameters of the network
according to φ<sup>t+1</sup> = φ<sup>t</sup> - λ∇<sub>φ</sub>R<sub>emp</sub>(f<sub>φ</sub>), where λ is referred to as
the learning rate. Once the optimization is complete and we ob-
tain the solution φ, it is natural to assess the quality of the trained
model f<sub>φ</sub> on an independent testing dataset. The empirical risk
evaluated on the testing set is often larger than on the training set,
and large differences indicate overfitting, which indicates that the
model does not generalize well to the unseen data. The ability to
accurately predict on unseen data is referred to as generalization
and the empirical risk on the test data provides a measure of the
generalization error. In order to reduce the generalization error
one might explore different model choices (e.g. neural network ar-
chitectures), additional regularization terms in the loss function,
different learning rates, optimization algorithms, or early stop-
ping criterion in the optimization. Once trained, the model can
be applied to data. In order to produce a binary electron vs. pro-
ton decision from the continuous output of the neural network,
we must threshold (i.e. classify as proton if f<sub>φ</sub>(x) > k). The
choice of the threshold k is often referred to as a working point
and it sets the tradeoff between electron and proton efficiency,
fake-rates, purity, etc. These familiar concepts in particle physics
are usually referred to in different terms in machine learning and
a receiver operating characteristic curve, or ROC curve, is used
to summarize the tradeoff between true positive rate (TPR) and
false positive rate (FPR). Importantly, the characterization of the
efficiency / rejection (or equivalently the ROC curve) requires la-
beled data. In a particle physics context, it is recognized that
the simulation is not perfect and the mismodelling is associated
to the presence of systematic uncertainty. In machine learning,
the discrepancy between the distribution of the training dataset

and the distribution of the data that the model will be applied
to in practice is referred to as *domain shift* or *distribution shift*.
While mismodelling in the training dataset might lead to a less-
than-optimimal classifier in practice, the real source of systematic
uncertainty comes from mismatch between the data used to char-
acterize the performance of the classifier and the unlabeled data
that the classifier is applied to. This motivates the use of data-
driven methods to calibrate the resulting model.

This example provides a vertical slice through the various as-
pects of supervised machine learning in particle physics. Now we
factorize and abstract the various ingredients in order to provide
a more general treatment with a broader scope.

41.2 Fundamental concepts

41.2.1 Loss, risk, empirical risk

The term learning in machine learning generally refers to op-
timization of some objective, which can be thought of as maxi-
mizing utility or minimizing risk. The risk brings together three
main ingredients. The first is the model family F (where f ∈ F
is the quantity that we vary during optimization), the second is
the loss function L, and the third is a data distribution p(u). The
risk for a model f ∈ F is defined as its expected loss

$$
\mathcal{R}[f] := \mathbb{E}_{p(u)}[\mathcal{L}(u, f(u))] \equiv \int \mathcal{L}(u, f(u)) p(u) \, du, \quad (41.1)
$$

where $\mathbb{E}_p[\cdot]$ refers to the expectation with respect to the distribution $p$. In the context of supervised learning, the distribution $p(u)$ describes a joint distribution over the features $x$ and the labels $y$ (i.e. $p(u) = p(x, y)$), the model only depends on the features $f(u) = f(x)$, and the loss function takes on the special form $\mathcal{L}(u, f(u)) = \mathcal{L}(y, f(x))$. In the context of unsupervised learning there are no labels, and $u = x$. Written this way, the risk is a functional, and the idealized goal for machine learning is to solve the optimization problem

$$
f^* = \arg \min_{f \in \mathcal{F}} \mathcal{R}[f], \quad (41.2)
$$

where $\mathcal{F}$ would include all possible functions.

One of the defining characteristics of machine learning in prac-
tice is that one does not know the data distribution p(u), but does
have access to samples from that distribution, i.e. {uᵢ}ᵢ=₁,...,n
with uᵢ i.i.d. p(u). This leads to the corresponding empirical risk

$$
\mathcal{R}_{\text{emp}}[f] := \mathbb{E}_{\hat{p}(u)}[\mathcal{L}(u, f(u))] \equiv \frac{1}{n} \sum_{i=1}^{n} \mathcal{L}(u_i, f(u_i)), \quad (41.3)
$$

where $\hat{p}(u) = \frac{1}{n} \sum_{i=1}^{n} \delta(u - u_i)$ is referred to as the empirical distribution of the dataset $\{u_i\}_{i=1,...,n}$. The empirical risk minimization principle is a core idea in statistical learning theory [7], which approximates $f^*$ with its empirical analogue

$$
\hat{f} = \arg \min_{f \in \mathcal{F}} \mathcal{R}_{\text{emp}}[f], \quad (41.4)
$$

where $\hat{\mathcal{F}}$ are all possible functions parametrized by the model
parameters $\phi$. In an idealized infinite parameter limit machine
learning functions, such as neural networks, are universal approx-
imators, such that they cover all functions and $\mathcal{F} = \hat{\mathcal{F}}$. For finite
size networks this may or may not be a valid assumption. Expres-
sivity of the network characterizes this universality property and
is a function of the network architecture and its parameters such
as width and depth of neural network layers. If the expressivity is
too small it leads to underfitting. However, an equally important
consideration is the risk of overfitting if we optimize equation 41.4
for too long.

While the loss function may quantify some well-motivated no-
tion of (negative) utility, it is also common to design loss functions
so that $f^*$ has some desired property. While in practice one does
not know the data distribution $p(u)$, it is constructive to imagine
that one does and analyze Eq. 41.2 with the calculus of variations.
In Secs. 41.3 we will consider several such loss functions where one
can show that the corresponding $f^*$ has the desired property even
if the form of the loss is not obvious from the point of view of utility. Furthermore, there are often multiple loss functions that can lead to the same $f^*$. Then one can think of machine learning as applied calculus of variations where one solves Eq. 41.4 with a sufficiently flexible model, powerful optimization algorithms, and practical considerations to break the degeneracy between different loss functions that lead to the same $f^*$.

### 41.2.2 Generalization

With a sufficiently flexible model, it is possible to fit the training dataset very well, though the model might not generalize well to unseen data due to overfitting. More concretely, for a non-negative loss function one might have $\mathcal{R}_{\text{emp}}[\hat{f}] \to 0$, while the true risk might be large ($\mathcal{R}[\hat{f}] \gg 0$). The gap between the $\mathcal{R}[\hat{f}] - \mathcal{R}_{\text{emp}}[\hat{f}]$ is typically referred to as the excess risk<fnref n="1" />. While it is generally not possible to evaluate $\mathcal{R}[\hat{f}]$ exactly because we do not know $p(u)$, we can use an independent testing dataset (also called validation dataset) to obtain an unbiased estimate of it. This cross-validation method motivates the test-train split of the data.

Intuitively, a model with many parameters has more flexibility and is more prone to overfitting. A common and intuitive heuristic is that one should not fit a model with more parameters than there are data points. However, a more careful treatment reveals that this heuristic can be both pessimistic and optimistic. For example, the single-parameter model $f_{\phi}(x) = \text{sign}(\sin(\phi x))$ can perfectly classify any assignment of labels on data $(x_i, y_i)$ with $x \in \mathbb{R}$ and $y \in \{-1, 1\}$ and generalize poorly. Conversely, sometimes highly over-parameterized models (that have large subspaces of their parameters where $\mathcal{R}_{\text{emp}}[\hat{f}_{\phi}] \to 0$) might generalize well [8, 9]. Often this is achieved through regularization, both explicit and implicit (section 41.3.1.2).

Structural risk minimization is a modification to the empirical risk minimization principle that was introduced by Vapnik and Chervonenkis to account for the potential for overfitting [7]. However, the bounds are based on a worst-case type analysis and are often very weak. Recall that while one cannot calculate the true risk, one can obtain an unbiased estimate of it with a held-out, independent testing sample. Thus one can empirically compare the generalization error of two models and find that one generalizes better than the other even if the bounds might suggest the opposite. One of the major conceptual shifts that happened with the rise of deep learning was to more fully appreciate that these bounds and structural risk minimization were not a good learning principle in practice and that more theoretical work is needed to close the gap between formal bounds and empirical estimates of generalization error.

## 41.3 Common tasks and their associated loss functions

We now move to common tasks encountered in machine learning and their associated loss functions.

### 41.3.1 Supervised learning

Supervised learning generally refers to the class of problems where the training dataset are presented as input-output pairs $\{x_i, y_i\}_{i=1,...,n}$, where $x_i \in \mathcal{X}$ are the input features and $y_i \in \mathcal{Y}$ are the corresponding target labels. Furthermore, it is typically assumed that $(x_i, y_i) \stackrel{i.i.d.}{\sim} p(x, y)$, though $p(x, y)$ is usually not known explicitly. Finally, the loss function in supervised learning takes on the special form $\mathcal{L}(y, f(x))$. The resulting trained model is then used to predict the labels for dataset where labels are not available.

In what follows we will make use of the equalities $p(x, y) = p(x|y)p(y) = p(y|x)p(x)$ and Bayes theorem, which is discussed in Sec. 39.1. We will also use the notation $x \sim p(x)$ to indicate that the random variable $x$ is being drawn, sampled, or generated from the distribution $p(x)$.

#### 41.3.1.1 Regression

The goal of regression is to predict a label $y \in \mathcal{Y}$ given an input feature vector $x \in \mathcal{X}$. Typically, the label is a real-valued scalar, but $\mathcal{X}$ can be $\mathbb{R}^d$ or some more structured target (e.g. an image, sequence, graph, quantile, or distribution). When $\mathcal{Y}$ is discrete, the task is usually referred to as classification; however, the two are closely related and logistic regression is an example where the model predicts a continuous probability associated to the possible label values. In elementary statistical language, the target label $y$ is often called a dependent variable, while the feature $x$ is called the independent variable. In classical statistics, one often assumes a model for the data such as

$$y_i = f_{\phi}(x_i) + e_i, \quad (41.5)$$

where $e_i$ is an additive error term that is often assumed to be independent of $x$ and often assumed to be normally distributed. This leads to classic approaches like least-squares (see Sec. 40.2.3), and when the model $f_{\phi}$ is linear in $\phi$ (not in $x$!) one has linear regression that has a closed-form solution. However, we can relax these assumptions and consider the general case of an arbitrary joint distribution $p(x, y)$, which can be written as $p(y|x)p(x)$ without loss of generality. Consider the squared error as a loss function, which leads to the mean-squared error (MSE) for the empirical risk:

$$\mathcal{L}_{\text{MSE}}(y, f(x)) = (y - f(x))^2. \quad (41.6)$$

One might expect that the squared error would only be appropriate in the case that the conditional distribution $p(y|x)$ is normally distributed, but one can use the calculus of variations to show that in general

$$f_{\text{MSE}}^{*}(x) = \mathbb{E}_{p(y|x)}[y], \quad (41.7)$$

that is the optimal regressor for the MSE is the conditional expectation of $y$ given $x$.

One issue with the squared-error as a loss function is that it is sensitive to outliers. Alternatively, one can use the absolute error $|y - f(x)|$ as a loss function<fnref n="2" />. However, the discontinuous derivative of the absolute (L1) error leads to challenges in optimization. As a result there are various other loss functions, such as the Huber loss, that aim to be both robust and more amenable to optimization that we do not discuss here.

Note that these framing of regression yields a function $f(x)$ that only provides a point estimate for $y$. An alternative approach to regression is to model the full conditional distribution $p(y|x)$. One such example is Gaussian Process regression, which is discussed in Sec. 41.5.1. In that probabilistic approach, one can still obtain a point estimator, such as the conditional expectation or the maximum a posteriori (MAP) estimator

$$f^*(x) = \arg\max_y p(y|x), \quad (41.8)$$

and one can also derive uncertainty estimates on the predicted value $y$. In this setting, the prior distribution on the model family is closely related to the concept of regularization, which we touch on in Sec. in Sec. 41.3.1.2 and in Sec. 41.5.1.

When one directly models $p(y|x)$, or goes further to model the joint distribution $p(x, y) = p(y|x)p(x)$, then one can use maximum likelihood for the loss function. In that approach, the problem is really one of density estimation, which is a type of unsupervised learning that we discuss in Sec. 41.3.2.1. These two approaches are a classic examples of two different approaches to modelling. Regression with $f_{\text{MSE}}^*(x)$ is the prototypical example of descriptive modelling, while modelling the joint distribution is a prototypical example of generative modelling. Generally, discriminative approaches with supervised learning out perform generative approaches when there is sufficient data, but generative approaches can be beneficial in data-starved settings [10].

<footnote n="1">¹Similarly, the gap between the true risk of the learned model and the true risk of the optimal model (i.e. $\mathcal{R}[\hat{f}] - \mathcal{R}[f^*]$) is referred to as the regret. This quantity is mainly of interest for theoretical analysis of machine learning algorithms, and not of practical concern since usually neither term is tractable.</footnote>

<footnote n="2">²The absolute error and squared error are often denoted as L1 and L2 errors, respectively, in reference to the corresponding norms.</footnote>
### 41.3.1.2 A note on regularization

The trained model $\hat{f}$, or equivalently, the parameters of the trained model $\phi$ can be thought of as point estimates of $f^*$, and there is a correspondence to the issues of bias and variance discussed in Sec. 40.2 on parameter estimation. Generically, there is a bias-variance tradeoff, and when the number of parameters is large and the number of data points is not much larger, introducing a small bias can often lead to a significant reduction in variance. This motivates the explicit addition of a regularization term to the loss function, which will introduce some bias $f_{\text{reg}}^* \neq f^*$. A common form for of regularization is to penalize by the L2 norm of the parameters (i.e. $\|\phi\|^2$), which is referred to as Tikhonov regularization. This appears in the form of penalized maximum likelihood, and it is also commonly used in unfolding [11]. One can also interpret the regularization term as an explicit prior on the parameters, and the resulting model as the Bayesian maximum a posteriori (MAP) estimator. When paired with linear regression this is known as ridge regression, and when paired with kernel machines (see Sec. 41.5.1) this gives rise to kernel ridge regression or Gaussian process regression.

Another form of regularization is to restrict the model class $\hat{f}$. For example, a neural network and a sequence of narrow step functions (delta functions) can both be shown to be universal approximators in infinite parameter size limit, but on real world examples the former generalizes much better than the latter. Within the class of neural network models, convolutional neural networks are a subset of generic feedforward neural networks that enforce translational symmetry (see Sec. 41.5.3.5 for more discussion). Similarly, one might restrict to Lipschitz continuous functions. These types of choices are often encoded in the architecture of a neural network and are broadly referred to as inductive bias in the model.

In addition to explicit regularization terms in the loss function or through restrictions to the model class, it is also possible to regularize implicitly. One implicit regularization is through early stopping [11, 12], where we monitor the loss on training dataset and the loss on test dataset. While the training dataset loss continues to decrease with more gradient descent cycles, the test loss may not, and early stopping stops the training when test loss flattens out or begins to increase. Another powerful form of regularization used in deep learning models is known as dropout [13], which randomly removes some parts of the model during training and can be thought of as implementing a type of model averaging [14]. What is more surprising is that in the case of highly over-parameterized models where there is a large degenerate parameter space that achieves zero loss, $\Phi_0 = \{\phi | \mathcal{R}_{\text{emp}}[f_\phi]\} = 0$, that the dynamics of the optimization algorithm that is used will break the degeneracy and favor some particular $\hat{\phi} \in \Phi_0$ as if an additional regularization term was secretly included. Despite zero loss and over-parametrization, the corresponding generalization error may be small, a phenomenon called benign overfitting [15]. Importantly, the dynamics of different optimization algorithms will have different implicit regularization effects, and thus favor different parameter points in $\Phi_0$ that will have different generalization error [16]. Understanding this interaction is a topic of contemporary research in machine learning [17].

### 41.3.1.3 Classification

The goal of classification is to predict one of a finite number of class labels $y \in \mathcal{Y}$ given an input feature vector $x \in \mathcal{X}$. It is similar to regression in this way, but the focus is on discrete target space $\mathcal{Y}$. An important special case is when the label can only take on one of two values (e.g. “signal” or “background”), which is referred to as binary classification and is equivalent to simple hypothesis testing in statistics. It is common for a classifier to be the composition of a model $g: \mathcal{X} \to \mathbb{R}^{|\mathcal{Y}|}$ that predicts continuous probabilities for each class (i.e. $g(x) \approx p(y|x)$) followed by an operation that then chooses the discrete label $y \in \mathcal{Y}$, such as a fixed threshold or $f(x) = \arg\max g(x) \approx \arg\max_y p(y|x)$. This is the case for both classical methods like logistic regression and modern, deep learning approaches to classification; therefore, we will use the term probabilistic classifier for $g(x)$ or just classifier

when it is clear in context.

An intuitive loss function for classification is the zero-one loss, which simply counts the number of mis-classifications:

$$ \mathcal{L}_{0/1}(y, f(x)) = \begin{cases} 0, & \text{if } f(x) = y \\ 1, & \text{otherwise.} \end{cases} \qquad (41.9) $$

The zero-one loss can also be written as $\mathcal{L}_{0/1}(y, f(x)) = \mathbf{1}(y \neq f(x))$, where $\mathbf{1}(\cdot)$ is the indicator function. The zero-one loss is non-differentiable, so it does not pair well with gradient-based optimization.

For binary classification, one can use $y = \{0, 1\}$ as numerical values for the class labels and the mean-squared error $\mathcal{L}_{\text{MSE}}(y, f(x))$ in Eq. 41.6 for the loss function. The resulting model will approximate $f_{\text{MSE}}^*$, the conditional expectation of Eq. 41.7 takes on the form

$$ \begin{aligned} f_{\text{MSE}}^{*}(x) &= \mathbb{E}_{p(y|x)}[y] \to p(y=1|x) \\ &= \frac{p(x|y=1)p(y=1)}{p(x|y=0)p(y=0) + p(x|y=1)p(y=1)}. \end{aligned} \qquad (41.10) $$

That is the MSE loss for binary classification leads to the Bayesian posterior probability that the label $y=1$ given the feature vector $x$.

Equation 41.10 highlights an important general feature of supervised learning relevant for particle physics, which is that the joint distribution $p(x, y)$ of the training dataset implies a prior distribution $p(y)$ on the labels or classes. This prior distribution need not reflect a degree of belief or the frequency in real data, it represents the frequency in the training dataset. However, it is important to keep in mind that when applying the resulting model to a different dataset with the same conditional distribution (data likelihood) $p(x|y)$ for the features and a different prior $p'(y)$ for the labels that the probabilistic interpretation of the result will not be properly calibrated. A common choice for binary classification is to use a balanced training dataset with $p(y=0) = p(y=1) = \frac{1}{2}$, while in many cases the true $p'(y=1)$ in the experimental data might be very small (i.e. low signal-to-background), unknown, or zero (i.e. a hypothetical particle that does not exist).

If $p'(y)$ and $p(y)$ are known then the Bayes theorem can be used to re-calibrate the posterior $p(y|x)$ from one prior to another. One example of such re-calibration is the correspondence of binary classification to simple hypothesis tests in frequentist statistics discussed in Sec. 40.3.1 of the Statistics chapter. In that setting, the Neyman-Pearson lemma states that the optimal classifier is given by the likelihood-ratio. Adapting to the notation of this chapter, we have

$$ f_{\text{N.P.}}^{*}(x) = \frac{p(x|y=1)}{p(x|y=0)}, \qquad (41.11) $$

which does not depend on the prior probabilities $p'(y=0)$ or $p'(y=1)$ as in Eq. 41.10, or, equivalently, assumes equal priors $p'(y=0) = p'(y=1)$. From Bayes theorem, we have the identity $p(y=1|x)/p(y=0|x) = [p(x|y=1)p(y=1)]/[p(x|y=0)p(y=0)]$, which can be used to show that the two functions are related by a one-to-one, monotonic transformation

$$ f_{\text{N.P.}}^{*}(x) = \frac{p(y=0)}{p(y=1)} \frac{f_{\text{MSE}}^{*}(x)}{1 - f_{\text{MSE}}^{*}(x)}, \qquad (41.12) $$

which is referred to as the likelihood-ratio trick, which plays an important role in simulation-based inference (see Sec. 41.3.5). Importantly, the monotonic transformation does not impact the tradeoff of type-I and type-II error (or, equivalently, the FPR and TPR), therefore the ROC curve for $f_{\text{N.P.}}(x)$ and $f_{\text{MSE}}(x)$ are identical and do not depend on the prior probabilities $p(y)$. This property has been leveraged in the context of weakly supervised approaches [18] and enables one to train a classifier on data without access to labels as long as one has two datasets with different ratios $p(y=1)/p(y=0)$ and the same conditional distribution $p(x|y)$ of the features given the labels.
A generalization of the binary loss function for classification of Eq. 41.10, which applies to multiple classes, is the cross-entropy loss

$$ \mathcal{L}_{\text{xe}}(y, f(x)) = - \sum_{c \in |\mathcal{Y}|} \mathbf{1}(y=c) \log(f_c(x)), \quad (41.13) $$

where $f: \mathcal{X} \to \mathbb{R}^{|\mathcal{Y}|}$ and the indicator function picks out the term in the sum for the corresponding class label $y$. This loss can be derived from maximizing the posterior of Eq. 41.57 using a discrete set of class labels $y$, which identifies $f_c(x) = \tilde{f}(y=c|x) = p(y=c|x)$ and thus enforces the constraint $\sum_c f_c(x) = 1$ and $f_c(x) \ge 0$ (for all $x \in \mathcal{X}$, e.g. by using the softmax function). The notation is aligned with the interpretation of $\tilde{f}(y|x)$ as a conditional distribution, i.e. an approximation to the true posterior $p(y|x)$. The risk associated to the cross entropy loss function is

$$ \begin{aligned} \mathcal{R}_{\text{xe}}[f] &= \mathbb{E}_{p(x,y)} \left[ - \sum_{c \in |\mathcal{Y}|} \mathbf{1}(y=c) \log f_c(x) \right] \\ &= - \sum_{c \in |\mathcal{Y}|} p(y=c) \mathbb{E}_{p(x|y)} [\log \tilde{f}(y=c|x)]. \end{aligned} \quad (41.14) $$

This is equivalent to $\mathcal{R}_{\text{xe}}[f] = \mathbb{E}_{p(x)}[H[p(y|x), \tilde{f}(y|x)]]$, where $H[p, f] = \mathbb{E}_{p}[-\log f]$ is the cross-entropy between the two distributions. One can use a Lagrange multiplier $\lambda$ to enforce the normalization constraint (e.g. equation 41.10) and the calculus of variations to show that

$$ f_{x.e.,c}^{*}(x) = \lambda p(x, y=c) = \lambda p(y=c|x)p(x) = p(y=c|x). \quad (41.15) $$

This approach is closely related to the loss functions that are used for density estimation, the forward Kullback-Leibler (KL) divergence, and the maximum likelihood estimation. Minimizing cross entropy $H[p, f_{\phi}]$ with respect to $\phi$ is equivalent to minimizing the forward KL divergence

$$ KL(p||f_{\phi}) := \mathbb{E}_{p}[\log p(x)] - \log f_{\phi} = H[p, f_{\phi}] - H[p], \quad (41.16) $$

where $H[p] := \int p(x) \log p(x) dx$ is the entropy and independent of $f_{\phi}$. The KL divergence $KL[p||f] \ge 0$, and equal if and only if $p = f$.

One can also consider the reverse KL divergence $KL[f_{\phi}||p]$, which is also minimized by $f_{\phi} = p$; however, this requires one to be able to generate samples $x_i \sim f_{\phi}(x)$ and be able to evaluate the probability density $p(x_i)$. Often this is not the case for real world data, but this approach is useful in the context of variational inference, where $f_{\phi}(x)$ is an approximation to the posterior of $x$, and $p(x)$ is the likelihood times the prior, which can be evaluated on samples $x_i$ drawn from $f_{\phi}$. Because likelihood times prior is not normalized this KL divergence optimizes the lower bound to the normalization (Evidence Lower Bound Objective or ELBO).

The forward KL is also closely related to the variational free energy principle in statistical mechanics where $H[f_{\phi}]$ represents the entropy of the variational distribution, $p(x) \propto \exp(-E(x)/kT)$ is the Boltzman factor for the state $x$ with energy $E(x)$, and $H[f_{\phi}, p] = \mathbb{E}_{f_{\phi}}[E(x)]$ represents the expected energy for the variational distribution. As is well known to physicists, minimizing the free energy involves a balance between minimizing the energy and maximizing the entropy.

In some cases it is possible to augment the training dataset with an unbiased, stochastic estimate of $p(y=c|x)$ that we denote $s_c(x, z)$. For example, when the simulation involves latent variables $z$ (i.e. Monte Carlo truth quantities), then the simulation encodes $p(x, z|y)$, the joint distribution over the observed features $x$ and the latent variables $z$ conditioned on the class $y$. In many cases the simulation evolves through a Markov process (e.g. the detector response only depends on the momenta of the incoming particles $z$, not the details of the hard scattering $y$). In that case, it is often possible to calculate $s_c(x, z) = p(y=c|x, z)$ for each training sample [19]. Using the identity, $\mathbb{E}_{p(z|x)}[p(y|x, z)] = p(y|x)$, we see that $s_c(x, z)$ is an unbiased estimator of $p(y=c|x)$.

In this case, one can use $s_c(x, z)$ in place of the indicator function in Eq. 41.13 to construct an improved (lower-variance) loss function that reproduces the same cross-entropy risk function$^{3,4}$

$$ \begin{aligned} \mathcal{R}'_{\text{xe}}[f] &= \mathbb{E}_{p(x,y,z)} \left[ - \sum_{c \in |\mathcal{Y}|} s_c(x, z) \log f_c(x) \right] \\ &= - \sum_{c \in |\mathcal{Y}|} \mathbb{E}_{p(x)}[p(y=c|x) \log f_c(x)], \end{aligned} \quad (41.17) $$

which yields the same optimal classifier $f_{x.e.,c}^{*}(x) = p(y=c|x)$ [20, 21].

We note that unlike in the binary classification case, the multi-class classifier is sensitive to the priors $p(y)$ used in training. This leads to complications as often the class proportions are unknown. For example, one might be interested in classifying a signal when multiple backgrounds are present and the relative proportion of those different background components is uncertain. Ideally one would like the class proportions for the background components used in training to match those in the data, which presents an additional training challenge if those proportions are heavily unbalanced.

### 41.3.2 Unsupervised learning

Unsupervised learning generally refers to the class of problems that use unlabeled training dataset $\{x_i\}_{i=1,...,n}$, where $x_i \in \mathcal{X}$ are the input features. Furthermore, it is typically assumed that $(x_i) \stackrel{i.i.d.}{\sim} p(x)$, though $p(x)$ is usually not known explicitly. Finally, the loss function in supervised learning takes on the special form $\mathcal{L}(x, f(x))$.

A related concept is that of self-supervised learning, which also aims to distill useful features in the data without requiring supervision labels for every sample in the input data. Self-supervised methods can make use of large unsupervised datasets and build meaningful representations by performing data augmentation, and learning the latent space mapping that is insensitive to it. For example, in the case of galaxy images one may augment the data by performing image rotations, adding noise, size scaling, adding point spread function smoothing etc., all of which are realistic transformations expected in a real galaxy image survey [22]. Self-supervised learning then learns the latent space representation where all of these augmentations of the same training sample result in the same latent space position. This training is augmented with contrastive learning, which ensures that different training samples do not all collapse to the same latent space position. When this latent space representation is used for downstream tasks such as classification it outperforms other forms of supervised learning [23].

#### 41.3.2.1 Density estimation

The goal of density estimation is to estimate a distribution $p(x)$ based on samples $\{x_i\}_{i=1,...,n}$ with $x_i \stackrel{i.i.d.}{\sim} p(x)$. Conceptually, this is the same goal as when fitting a parameterized distribution $f(x; \theta)$ to data using the method of maximum likelihood as described in Sec. 40.2.2 of the chapter on statistics. In practice, the difference in the machine learning context has to do with the flexibility of the model and the dimensionality of the data. A highly-flexible model, which can effectively approximate any distribution, is referred to as a non-parametric model (though, ironically, usually this means the model has many parameters). In contrast, typical maximum likelihood fits in particle physics are based on restricted families of distributions with relatively few parameters and the data is typically one- or two-dimensional, though occasionally five- or six-dimensional.

Maximizing the likelihood function in Eq. 40.10 $\mathcal{L}(\theta) =$

<footnote>$^3$A similar approach can also be used for the squared-error, see Ref. [20].</footnote>

<footnote>$^4$The right hand side of Eqs. 41.14 and 41.17 are written in a different form, but are equivalent.</footnote>
$\prod_{i=1}^{n} f(x_i; \theta)$ is equivalent to minimizing the empirical risk:

$$ \mathcal{R}_{\text{emp,xe}}[f_{\phi}] = -\frac{1}{n} \sum_{i=1}^{n} \log f_{\phi}(x), \quad (41.18) $$

where we adopt the notion used in this chapter. The loss is simply $\mathcal{L}(x, f_{\phi}(x)) = -\log f_{\phi}(x)$, and the corresponding risk is

$$ \mathcal{R}_{\text{xe}}[f_{\phi}] = \mathbb{E}_{p(x)}[-\log f_{\phi}(x)], \quad (41.19) $$

which is the cross entropy $H[p, f_{\phi}]$. For density estimation, the model is usually constructed to enforce $\int f_{\phi}(x)dx = 1$ and $f_{\phi}(x) \ge 0$ so that it can be interpreted as a distribution. With this constraint, one can show that $f_{\phi}^*(x) = p(x)$.

The concepts of generalization and overfitting are particularly acute in unsupervised learning, where the likelihood maximization of equation 41.18, combined with universal approximator assumption, must converge onto $\hat{p}(x) = \frac{1}{n} \sum_{i=1}^{n} \delta(x - x_i)$, the empirical distribution of the dataset $\{x_i\}_{i=1,...,n}$. This distribution has the highest likelihood on the training dataset and the lowest likelihood on the test data where it gives $\hat{p}(x) = 0$. as long as the test dataset are not identical to the training dataset. So the empirical distribution of the training dataset has the worst possible generalization property, yet it is the solution we converge to for sufficiently expressive architectures in the absence of any regularization. In contrast, in supervised learning we often observe the phenomenon of benign overfitting, where even zero loss can generalize well.

In addition to approaches to density estimation that involve learning in the sense of minimizing a loss or risk function, we note that there are also classical density estimation techniques such as histogramming and kernel density estimation [24–26].

### 41.3.2.2 Representation learning, compression, and auto-encoders

A recurring topic in machine learning and statistics is how to represent the data. Much of classical statistics involves constructing a low-dimensional summary statistic that extracts the relevant information from the data for a particular task (a sufficient statistic in language of classical statistics). There is a spectrum of representations with trade-offs. At one end of this spectrum is lossless compression that allows you to encode the data into a smaller, intermediate representation that carries all the information since it can be decoded back into the original data. At the other end of the spectrum is something like the likelihood-ratio, which is a single scalar that carries the relevant information needed for hypothesis testing for a single hypothesis, but it discards all the other information that might be needed for other tasks (such as testing other hypotheses). An intermediate point in this spectrum is the process of feature engineering, which refers to the creation of new features $\mathcal{X}'$ from the original features $\mathcal{X}$ in hopes that the down-stream task will be easier with the new features. For example, instead of working directly with the energy and momentum of particles, one might compute invariant masses, angles between particles, etc. This type of feature engineering generally improved performance for shallow neural networks, decision trees, etc.; however, with the rise of deep learning this is often no longer necessary and often limits performance compared to working with the original features. One can think of the intermediate layers of a neural network between the input and the output a representation of the data that is good for the task at hand, and by training all the layers of the network simultaneously (or “end-to-end”) one can see the intermediate layers as a learned representations. For a review see Ref. [27].

An example of a linear dimensionality reduction representation and data compression is principal component analysis (PCA) of data $x \in \mathbb{R}^d$ at fixed latent space dimensionality $k$ ($k<d$), which finds the orthogonal linear transformation, $O$,

$$ O: \mathbb{R}^k \to \mathbb{R}^d, z \mapsto Oz, OO^T = I_d \quad (41.20) $$

that maximizes the data variance in the latent space. Maximizing the variance of the transformed data is equivalent to minimizing

the average reconstruction error (the residual variance in data space),

$$ \mathcal{L}_{\text{a.e.}}(x, f(x)) = \|x - f(x)\|^2. \quad (41.21) $$

A PCA can thus be interpreted as a linear, orthogonal model that is trained to minimize the $L_2$-distance between the input data and the reconstructed data given the fixed dimensionality $k$. In practice, the PCA problem can be solved analytically without the use of optimization algorithms or the loss function: the principal components are given by the eigenvectors of the data covariance matrix.

A suitable latent space dimensionality, $k$, is chosen by ordering the eigenvalues, $\lambda_i$, of the data covariance in descending order, and keeping only the first few eigenvectors that correspond to the largest eigenvalues. The cut is often made at dimensionality that capture around 90% of the data variance. For many data sets this results in $k \ll d$. The average reconstruction error that originates from the discarded eigenvalues is $\sigma_{\text{recon}}^2 = \sum_{i=k+1}^{d} \lambda_i$.

Another common type of representation learning and nonlinear dimensionality reduction is based on the auto-encoder $f = g \circ e: \mathcal{X} \to \mathcal{X}$, where $e: \mathcal{X} \to \mathcal{Z}$ is referred to as the encoder and $g: \mathcal{Z} \to \mathcal{X}$ is referred to as the generator or decoder. Typically the dimensionality of $\mathcal{Z}$ is much less than $\mathcal{X}$, and $z = e(x)$ can be thought of as a compressed representation of the input. The intermediate space $\mathcal{Z}$ is sometimes referred to as the bottleneck or the latent space of the auto-encoder. If the bottleneck is sufficiently large and the encoder and decoder are sufficiently flexible, then the function $f$ could just be the identity (i.e. lossless compression). However, if the encoder and decoder are not sufficiently flexible or the dimensionality of the latent space is not large enough there will be some reconstruction error. Therefore the reconstruction error of Eq. 41.21 serves as a natural loss function of an auto-encoder.

Once trained, the encoder $e(x)$ can be used independently of the decoder to provide a generic low-dimensional representation of the data. The flexibility of this approach is attractive; however, there are no guarantees that this representation will be optimal for the other task. Indeed, the transition from pre-trained auto-encoders to end-to-end learning is one of the important trends that characterized the onset of the deep learning era.

While achieving zero reconstruction error may seem good as it would imply lossless compression, it often performs poorly in practice. First, the encoder may be overfit to the training dataset and not generalize well to held out data. Secondly, it may not be robust to domain shift (see Sec. 41.7.2). One approach to address these issues is the denoising auto-encoder, which uses a form of regularization that corrupts the input with noise $x' \sim q(x'|x)$ while still targeting reconstruction of the uncorrupted input $x$.

$$ \mathcal{L}_{\text{d.a.e.}}(x, f(x)) = \|x - f(x')\|^2 \quad \text{with} \quad x' \sim q(x'|x), \quad (41.22) $$

where $q(x'|x)$ is some probability distribution such as a multivariate normal.

### 41.3.2.3 Clustering

The goal of clustering is to group the data $\{x_i\}_{i=1,...,n}$ into $k$ groups, or clusters, usually with $k \ll n$. Intuitively, if two data points belong to the same cluster, then they should be similar in some sense. Conversely, if two data points are very different, then they should be assigned to different clusters. The notion of similarity usually is based on some heuristic, and there are a variety of algorithmic and probabilistic clustering algorithms. In some cases $k$ is specified, while in others it is determined by the clustering algorithm. There is also a distinction between flat clustering that directly partitions the data into $k$ clusters and hierarchical clustering where clusters are nested hierarchically as the name suggests. In many cases clustering uses some notion of distance $d(x_i, x_j)$, which may be the $L_p$ norm $\|x_i - x_j\|_p$.

One of the most common clustering algorithms is known as k-means, where $k$ is specified by the user and results in sets $S = \{S_1, ..., S_k\}$ that minimize the variance of each cluster. Thus,
the objective is

$$
\begin{align}
\arg\min_{\mathbf{S}} \sum_{i=1}^{k} \sum_{\mathbf{x} \in S_i} \|\mathbf{x} - \boldsymbol{\mu}_i\|^2 &= \arg\min_{\mathbf{S}} \sum_{i=1}^{k} |S_i| \mathrm{Var} S_i \tag{41.23} \\
&= \arg\min_{\mathbf{S}} \sum_{i=1}^{k} \frac{1}{2|S_i|} \sum_{\mathbf{x}, \mathbf{y} \in S_i} \|\mathbf{x} - \mathbf{y}\|^2 \nonumber
\end{align}
$$

where $\mu_i$ is the mean of points in $S_i$.

While k-means and many other clustering algorithms are de-
fined in terms of an optimization problem, the optimization ob-
jective is often not representative of the actual notions of per-
formance in a given context. As a result, there are a number of
quantities used for the evaluation and assessment of the result-
ing clustering. These include Davies-Bouldin index, Dunn index,
Rand index, Jaccard index, purity, F-measure, Hopkins statistic,
etc. [28].

41.3.3 Optimal control, reinforcement learning, and active learning

Many problems in science and engineering can be cast as a con-
trol problem, which comprises a cost functional that is a function
of state and some control variables that specify some underlying
dynamical system. This is relevant for the control of accelera-
tors where the dynamical system is physical. This formalism can
also be used to describe the design of experiments, planning of
an observational survey, and other decision making processes rel-
evant to the scientific method. There is a tremendous amount of
literature on the subject, and it is closely connected to planning,
dynamic programming, and reinforcement learning. Optimal con-
trol can be seen as an extension of the calculus of variations, and
thus generalizes the framing of learning presented in Sec. 41.2.

**Optimal control** theory deals with finding a control for a dynamical system over a period of time such that the objective function is optimized. The underlying system can be discrete or continuous and may be deterministic or stochastic. The commonalities and differences between optimal control and reinforcement learning can be best understood through the formalism of a Markov decision process (MDP), which is a discrete-time stochastic control process. We follow common notation, such as that found in Wikipedia.

A Markov decision process comprises four components often organized as a 4-tuple (S, A, P_a, R_a), where: S is a set of states called the state space, A is a set of actions called the action space, P_a(s, s') = Pr(s_{t+1} = s' | s_t = s, a_t = a) is the probability that action a in state s at time t will lead to state s' at time t + 1, R_a(s, s') is the immediate reward (or expected immediate reward) received after transitioning from state s to state s', due to action a.

The policy function $\pi$ is a mapping from state space to action space that can be either deterministic or probabilistic. For examples, the policy that drives a computer chess playing system, decides which move to make given the current state of the board. Similarly, policies dictate which experiment should be built next, which field of the sky should be observed, or how to adjust the operational parameters of an accelerator. The dynamics of the resulting system are then fixed by combining the policy with the underlying MDP. The evolution of the resulting dynamical system behaves like a Markov chain since the action chosen in state s is completely determined by $\pi(s)$ and $\Pr(s_{t+1} = s' | s_t = s, a_t = a)$ implies the Markov transition matrix $\Pr(s_{t+1} = s' | s_t = s)$.

The objective optimal control is to choose a policy $\pi$ that will maximize a cumulative function of the instantaneous rewards $R_a$. A common choice is the expected discounted sum:

$$
\mathbb{E} \left[ \sum_{t=0}^{\infty} \gamma^t R_{a_t}(s_t, s_{t+1}) \right], \quad (41.24)
$$

where $a_t \sim \pi(s_t)$ are the actions given by the policy,
the expectation computed with respect to the distribution
$s_{t+1} \sim P_{a_t}(s_t, s_{t+1})$, and $\gamma$ is the discount factor satisfying

$0 \le \gamma \le 1$. The discount factor is usually close to 1 and some-
times reparameterized as $\gamma = 1/(1+r)$, where r is called the dis-
count rate. A lower discount factor motivates the decision maker
to favor taking actions early, rather than postpone them indefi-
nitely.

A policy that maximizes the objective function is called an op-
timal policy and denoted $\pi^*$, though the optimal policy need not
be unique. Importantly, the Markov property implies that the op-
timal policy is only a function of the current state. Dynamic pro-
gramming can be used to find the optimal policy for MDPs with
finite state and action spaces. For instance, in value iteration
(a.k.a. backward induction) can be used to solve the “Bellman
equation” [29]. For continuous-time systems, the optimal policy
is defined by the Hamilton–Jacobi–Bellman equation [30].

In many settings, it is assumed that the state s is fully known when action is to be taken and there are no latent variables. When this assumption is not true, the problem is called a partially observable MDP. These problems are generally more difficult and the dynamic programming algorithms do not directly apply [31].

**Reinforcement learning** The main difference between the classical dynamic programming methods and reinforcement learning (RL) algorithms is that the latter do not assume knowledge of an exact mathematical model of the MDP and they target large MDPs where exact methods become infeasible. For example, RL was used in the context of jet physics to search for the most likely jet clustering when the number of constituents was too large for the exact dynamic programming algorithm to be used [32]. In addition, RL can be used when the probabilities or rewards are unknown. Instead, the transition probabilities are often accessed indirectly through interaction with a real or simulated environment. For example, in Q-learning one uses experience to estimate the array $Q(s,a)$ defined as

$$
Q(s, a) = \sum_{s'} P_a(s, s')(R_a(s, s') + \gamma V(s')). \quad (41.25)
$$

While this function is also unknown, it can be estimated during
learning based on (s, a) pairs together with the outcome s'.

Numerous variations to RL exist, which include so-called model-
based and model-free approaches (referring to models of the in-
stantaneous rewards and the state transitions) and on-policy and
off-policy (which describes how the actions taken during learning
are related to the current policy). See Ref. [33] for an introduction
and Ref. [34] for a recent review.

**Multi-arm bandits** Multi-arm bandit problems are a classic reinforcement learning problem where one tries to maximize the expected gain by allocating a limited set of resources to various alternatives. The name is a reference to a gambler with a fixed amount of money that must choose between multiple slot machines (or “one-armed” bandits) when the payoff for the individual machines is unknown. A hallmark of multi-arm bandit problems is that they involve a trade-off between exploration (playing machines to estimate their payoff) and exploitation (playing machines with the highest estimated payoff). Multi-armed bandits are used to manage large projects, organizations, and scheduling problems. The theory has a long history going back to Robbins in 1952 that used it to study the sequential design of experiments [35] and Gittins who derived an optimal policy under some conditions [36].

**Bayesian optimization** A closely related set of techniques involve optimizing some expensive black box function $f(x)$. For instance, the function may be computationally expensive to evaluate or low-latency, e.g. it may involve manually re-configuring a system. This is particularly relevant for analysis optimization in particle physics where evaluating $f(x)$ involves processing large numbers of simulated collisions. Another common use case involves optimizing the hyperparameters of a learning algorithm.

Without any assumptions about the function $f(x)$ this is hope-
less; however, if one assumes something about the functions (e.g.
some notion of smoothness) then one can leverage function evaluations. evaluations $\{f(x_t)\}_{t=1,...,T}$ to say something about what value the function might take on at other values of $x$. This is usually cast in Bayesian terms, and Gaussian processes (Section 41.5.1) are often used to model the distribution over $f(x)$. The optimization techniques that use this framing are generically referred to as Bayesian Optimization [37].

Optimization in this context is usually characterized by an exploration-exploitation trade-off, similar to what is found in multi-arm bandits. Here, exploration refers to function evaluations that characterize the function in regions that haven't been evaluated, while exploitation refers to evaluations near what is predicted to be its maximum. This setting is similar to reinforcement learning in that it involves sequential decisions (i.e. where to evaluate the function next), but usually the target function $f(x)$ is assumed to be static. In that sense, the state referred to in the language of an MDP is the state of knowledge about the function after sequential evaluations $\{f(x_t)\}_{t=1,...,T}$. The reward at time $t$ is not the value of the function $f(x_t)$, but some quantity that characterizes what was learned about the function's maximum. In this literature, one often refers to the acquisition function, which plays a similar role as the expected value of the reward in RL. Common acquisition functions include the probability of improvement, the expected improvement, and an upper-confidence bound [38].

### Connection to experimental design
In the context of physics experiments we often want to build the next generation experiment which can reach certain target precision on parameters of interest that the experiment can measure. To achieve this precision we may deploy the concept of experimental design, where we the objective is to optimize some objective that quantifies the expected performance of an experiment. Moreover, we would like to achieve this within some constraints such as a minimal cost. Such problems are often solved using the Fisher matrix optimization, where Fisher matrix can be viewed as the expectation of the inverse covariance matrix. We can define

$$t(x|\theta_0) := \nabla_{\theta} \log p(x|\theta)|_{\theta_0}, \qquad (41.26)$$

where $t(x|\theta_0)$ is referred to by statisticians as the score [20]. The score plays an important role in the classical statistics as it is a sufficient statistic when $p(x|\theta)$ is in the exponential family, and through the Rao-Cramér-Fréchet bound on the variance of an estimator for $\theta$, and is used to define the Fisher information matrix $I_{ij}(\theta) := E_p(x|\theta)[t_i(x|\theta)t_j(x|\theta)]$. The Fisher information, in turn is an important object in experimental design. In particle physics, the score is closely related to the concept of optimal observables. The corresponding diagonal element of the inverse of the Fisher matrix thus provides expected uncertainty estimation of a given parameter, and is a lower bound to the error of an unbiased estimator (Cramér-Rao bound). In experimental design we vary the experiment parameters within the constraints such as the total cost to minimize this uncertainty. This framework has been widely used in cosmology, where Fisher analysis is the foundation of any experiment proposal.

### Active learning
Active learning is closely related to Bayesian optimization, described above. In Bayesian optimization one estimates the function $f(x)$ from some set of evaluations $\{y_t = f(x_t)\}_{t=1,...,T}$; however, the goal is to find the maximum $x^* = \arg\max_x f(x)$. In active learning, the goal is not to find the maximum of $f(x)$, but to approximate the function as one does in supervised learning. The main difference compared to vanilla supervised learning is that the labeled training dataset isn't provided a priori in a passive way, but the learning algorithm actively decides where to generate ($x_t, y_t = f(x_t)$) pairs. The function $f(x)$ is sometimes referred to as an oracle. Active learning is particularly attractive when obtaining labelled data is a costly process.

More broadly, a challenge of many machine learning applications is obtaining labeled data, which can be a costly process. If a system could learn from small amounts of data, and choose by itself what data it would like the user to label via an external process called oracle, it would make machine learning more

powerful. Such frameworks are also called Experiment Design or Active Learning. In Active Learning, a model is trained on a small amount of data (the initial training dataset), and an acquisition function (often based on the model's uncertainty) decides on which data points to ask for a label. The acquisition function selects one or more points from a pool of unlabelled data points, with the pool points lying outside of the training dataset. Once we label the selected data points, these are added to the training dataset, and a new model is trained on the updated training dataset. This process is then repeated, with the training dataset increasing in size over time. The advantage of such systems is that they often result in dramatic reductions in the amount of labelling required to train an ML system (and therefore cost and time).

### 41.3.4 Anomaly detection and out-of-distribution detection

Unsupervised anomaly detection techniques detect anomalies in an unlabeled test data set under the assumption that the majority of the in-distribution data are normal under some measure, while out-of-distribution (OOD) data are not. In the context of auto-encoders a popular technique is to use the reconstruction error of Eq. 41.21 to identify an outlier as one with a large reconstruction error [39, 40]. One issue with this method is that for higher dimensional latent space and flexible neural network architectures the encoder-decoder map become identity for any input data, $f(x) = x$, regardless of whether input $x$ is from the in-distribution training dataset or from the out-of-distribution data. The choice of auto-encoder latent space dimensionality is thus an important hyperparameter that must be tuned.

Another set of anomaly detection techniques construct a model representing normal behavior from a given In Distribution training dataset, and then test the likelihood of a test instance to be generated by the utilized model. For instance, one can use density estimation methods such as normalizing flows (section 41.5.4.3) to learn the density (likelihood) of the In Distribution training dataset $p(x)$, and apply it to the test data. The expectation is that out-of-distribution data will have a lower density (likelihood) under the in-distribution density model. This expectation is however not always met in high dimensions and the method suffers because likelihood based training is sensitive to the smallest variance directions [41]: for example, a zero variance pixels leads to an infinite $p(x)$, and noise must be added to regularize it. But low variance directions may contain little or no information on the global structure of the image, so there is a mismatch between the training objective and outlier detection objective. A related issue is that of typicality: an in-distribution likelihood will typically be lower than the maximum value. For a Gaussian likelihood the typical set is on average $n/2$ in log likelihood below the peak value, where $n$ is the dimensionality of the data. So an out-of-distribution dataset that is closer to the peak would be preferred in terms of likelihood even though its distribution does not match the in-distribution data. If this happens by chance on low variance directions which dominate the likelihood, normalizing flows can assign higher likelihoods to out-of-distribution data than to in-distribution training dataset [42]. A number of techniques have been proposed to circumvent these limitations, such as likelihood regret [43], likelihood-ratio [41], likelihood in auto-encoder latent space [44], and Wasserstein distance training of the likelihood $p(x)$ [45]. These methods can achieve better anomaly detection performance than the auto-encoder reconstruction error [44, 46].

Supervised anomaly detection techniques require a data set that has been labeled as in-distribution and out-of-distribution and involves training a classifier (the key difference to many other statistical classification problems is the inherent unbalanced nature of outlier detection). These methods assume some form for what out-of-distribution data may look like, and their success relies on whether the assumed form is a realistic representation of actual out-of-distribution data. When this assumption is valid these methods can be more powerful than unsupervised methods, but the reverse is also true. A hybrid between the two approaches is to train a classifier without labels [47]. All these approaches are largely complementary to each other [48]. An example of different anomaly detection methods applied to HEP is LHC Olympics
2020 and Dark Machines challenges [49, 50].

### 41.3.5 Simulation-based inference

The goal of simulation-based inference (related to, but distinct from, likelihood-free inference) is to extend the statistical procedures described in the Chapter on Statistics (e.g. parameter estimation, hypothesis tests, confidence intervals, and Bayesian posterior distributions) to the situation where one does not know the explicit likelihood $p(x|\theta)$, the probability of the data given the parameters $\theta$, but has access to a simulator that defines the likelihood $p(x|\theta)$ implicitly [51, 52]. In a typical setup we would like to solve the so called inverse problem of getting the posterior of the parameters given the data, $p(\theta|x)$, but we cannot use Bayes theorem directly because we do not have explicit $p(x|\theta)$.

In particle physics, the simulators usually use Monte Carlo event generators (see Sec. 43) to sample unobserved latent variables $z$, such as the $z_p$ phase space of the hard scattering (see Sec. 49.4), $z_s$ associated to showering and hadronization, and $z_d$ associated to the interaction of particles with the detector (see Sec. 34). As such, the full simulation chain can be expressed symbolically as

$$ 
\begin{aligned}
p(x|\theta) &= \int dz p(x, z|\theta) \\
&= \int dz_d \int dz_s \int dz_p p(x|z_d)p(z_d|z_s)p(z_s|z_p)p(z_p|\theta),
\end{aligned} 
\quad (41.27) 
$$

where $\theta$ are the Lagrangian parameters that dictate the hard scattering. Evaluating the likelihood is intractable as it would require evaluating the integral above for each event.

While the likelihood is intractable, simulators provide the ability to generate synthetic data $x_i$ i.i.d. $p(x|\theta)$ for any value of the parameters $\theta$. One can use a suitable proposal distribution $\tilde{p}(\theta)$, sample $\theta_i$ i.i.d. $\tilde{p}(\theta)$, generate synthetic data $x_i \sim p(x|\theta_i)$, and then assemble a training dataset $\{x_i, \theta_i\}_{i=1,...,n}$ that can be used to train various machine learning models.

There is thus a close analogy between Simulation-based Inference and data driven machine learning tasks discussed so far, replacing $\theta$ with $y$. One difference is that in simulation-based inference we can always generate new samples by running additional simulations, while we typically view training dataset in machine learning as fixed. This property of Simulation-based Inference enables Active Learning, where the additional simulations are chosen such as to minimize the error on the desired statistical inference task. Another difference is that we often have access to the joint likelihood $p(x, z|\theta)$, where $z$ are unobserved latent variables$^5$.

Typically in particle physics, one uses histograms or kernel density estimation to model the distribution of observables (low-dimensional summary statistics such as the invariant mass) of simulated data [53]. Alternatively, one can use an explicit parametric family (such as a falling exponential, a Gaussian distribution, etc.) to model $\hat{f}(x|\theta) \approx p(x|\theta)$. That model is then used as a surrogate for the unknown density implicitly defined by the simulator. A related approach is known as Approximate Bayesian Computation (ABC), which approximates the likelihood through an acceptance probability that synthetic data is sufficiently close to the observed data [54, 55]. In practice, these techniques are limited to low-dimensional representations of the data. Thus the potential of recent machine learning approaches to simulation-based inference is to extend this approach to higher-dimensional data, while maintaining the already well-established statistical procedures.

For instance, one can use normalizing flows (see Sec. 41.5.4.3) and the loss functions for density estimation (see Sec. 41.3.2.1) to learn a surrogate model for the likelihood $\hat{f}(x|\theta) \approx p(x|\theta)$ [56]. Similarly, one can use conditional density estimation to learn a surrogate model for the posterior $\hat{f}(\theta|x) \approx p(\theta|x)$, which may involve including the prior-to-proposal ratio $\tilde{p}(\theta)/p(\theta)$ [57]. In addition to the unsupervised learning

techniques, one can also use supervised learning to learn the likelihood-ratio $r(x|\theta_0, \theta_1) = p(x|\theta_0)/p(x|\theta_1)$ by leveraging the likelihood-ratio trick of Eq. 41.12 [20, 58].

In some cases one can also augment the training dataset to include the joint likelihood-ratio

$$ r(x_i, z_i|\theta_0, \theta_1) := p(x_i, z_i|\theta_0)/p(x_i, z_i|\theta_1), \quad (41.28) $$

which can be used to reduce the variance for the squared-error loss or the improved cross-entropy estimator of Eq. 41.17 [20, 21]. While the marginal likelihood $p(x|\theta)$ is intractable due to the high-dimensional integral over the latent space, the joint likelihood is often tractable since no integration is necessary. In addition, one can often augment the training dataset with the joint score

$$ t(x_i, z_i|\theta_0) := \nabla_{\theta} \log p(x_i, z_i|\theta) \Big|_{\theta_0}. \quad (41.29) $$

Regressing on the joint score with the squared loss function $\mathcal{L}(t(x, z|\theta_0), f(x)) = (t(x, z|\theta_0) - f(x))^2$ corresponds to risk functional

$$ \mathcal{R}_{\text{score}}[f] := -\mathbb{E}_{p(x,z|\theta_0)}[(\nabla_{\theta} \log p(x_i, z_i|\theta)|_{\theta_0} - f(x))^2]. \quad (41.30) $$

One can show with the calculus of variations that the risk is minimized by the score of equation 41.26.

#### 41.3.5.1 Differentiable simulations

One of the approaches to make inference feasible in high dimensional latent space is to make simulations differentiable with respect to all of its parameters, global variables $\theta$ and latent variables $z$. Latent variable models are known in Machine Learning as generative models and are discussed below, but latent variable models are also common in scientific applications. As an example, in cosmology one can view the initial conditions of the large scale structure or cosmic microwave background as Gaussian distributed latent variables.

While differentiable simulations have traditionally not been developed for scientific applications, the success of Machine Learning, where backpropagation algorithm combined with gradient descent optimization (see Sec. 41.6.1) is the basis of its recent advances, has inspired a new look at this. One recent example is FlowPM cosmological N-body simulation, which takes advantage of Mesh-Tensorflow to achieve a GPU-accelerated, distributed, and differentiable simulation [59].

One broad class of inference problems where gradients make the problem simpler are the so called inverse problems. As a simple example let's assume the data $x$ are observed with some Gaussian noise with known variance $\sigma^2$, equal for all data points. The likelihood of the data can be written as

$$ \log p(x|z, \theta) = -\frac{N}{2} \log[2\pi\sigma^2] - \frac{\|x-g(z, \theta)\|^2}{2\sigma^2}. \quad (41.31) $$

Here $g(z, \theta)$ is the forward model (generator or decoder), i.e. the simulation output in the absence of noise. The joint distribution is $p(x, z|\theta) = p(x|z, \theta)p(z|\theta)$, where $p(z|\theta)$ is the prior of latent space variables, assumed to be known, with a known gradient with respect to $\theta$ and $z$. In this case the gradients $\nabla_{\theta} \log p(x, z|\theta)$ and $\nabla_z \log p(x, z|\theta)$ are available when the simulation forward model gradients $\nabla_z g(z, \theta)$ and $\nabla_z g(z, \theta)$ are available.

Availability of simulation gradients in turn enables gradient based optimization or sampling. In sampling approach one can use Monte Carlo Markov Chains to obtain the posterior samples of latent space $z$ and parameter space $\theta$ [60]. In optimization approach one can find the best-fit point $\hat{z}$, which is the MAP estimate of the latent space variables $z$ given $x$ and fixed $\theta$ [61]. Another optimization approach is Variational Inference, which attempts to model the posterior probability distribution using optimization. This is discussed further in Section 41.5.4.

#### 41.3.5.2 Unfolding as an inverse problem

While much of the work on simulation-based inference described above is aimed at inferring the parameters $\theta$ of the simulator, there is also work that aims to to infer the latent variables $z$. Here it is useful to think of the parameters $\theta$ as parameters of a theory, such as masses, coupling constants, or Lagrangian parameters, while $z$

<footnote n="5">⁵For this reason we prefer to use simulation-based inference instead of likelihood-free inference: joint likelihood $p(x, z|\theta)$ is often available, it is the marginal integral over latent space $z$ that is assumed to be intractable.</footnote>
might describe the kinematics of a collision before the detector response.

Inferring the distribution $p(z|\{x_1, \dots, x_n\})$ from a dataset of multiple observations is commonly referred to as unfolding in particle physics, and deconvolution in other contexts. Unfolding is a classic inverse problem, and the collection of ideas being used for machine-learning based simulation-based inference are also being applied in this setting. Examples of recent work exploring these approaches are Refs. [62–73]. In addition, there has been work to infer the posterior distribution for an individual event $p(z|x_i)$ using probabilistic programming techniques [74–76].

## 41.4 Data representations, inductive bias, and example applications

In Sec. 41.3 we describe the input data as living in an abstract space $x_i \in \mathcal{X}$. In this section, we briefly describe some of the common types of structured data that are encountered in physics and refer to the corresponding models classes that have been developed to work with them. We will describe the model classes in more detail in the following section.

The most basic and common type of data structure is when $\mathcal{X} = \mathbb{R}^d$. This is often referred to as tabular data since the entire data set $\{x_i\}_{i=1,\dots,n}$ can be thought of as a table with $n$ rows and $d$ columns. It is common to think of an individual entry $x_i$ as a vector in $d$-dimensional Euclidean space, where the coordinates correspond to the columns of this table. In some cases individual components of $x_i$ might be integers or take on only discrete values, in which case describing the space of the data as real-valued is a slight abuse of notation and representation. For many years this was the dominant type of data in high energy physics as it is a natural input type for shallow neural networks, multi-layer perceptrons, support vector machines, and tree-based methods found in popular tools such as TMVA [77].

For categorical data, one typically uses a numerical representation such as integer encoding where different categories are mapped to integers with a corresponding dictionary. Another common representation of categorical data is based on the so-called one-hot encoding (aka ‘one-of-K’ or ‘dummy’), in which case the category is mapped to a $k$-dimensional binary vector where $k$ is the number of categories and each component of this vector corresponds to a particular category. In the one-hot encoding, only one of the components is non-zero. Finally, there are approaches in one learns an embedding that maps discrete categories into $\mathbb{R}^d$; an example of this is Word2Vec [78]. Interestingly, such embeddings can preserve various types of semantics; for instance, the vector walking-walk is similar to the vectors swimming-swam as are the vectors connecting countries and their capital cities. This allows for a loose sense of algebra on the word embeddings such as walking-swimming+swam = walk. Similar types of embeddings have also been used in a number of scientific use-cases including biological sequences (e.g. DNA, RNA, and Proteins) for bioinformatics applications [79].

Particle physics data often is represented with an extension of the simple tabular data structure where the number of columns is not fixed. For instance, if the rows correspond to data for individual collisions, the number of electrons (and positrons) reconstructed in the event is variable. Thus the number of columns needed to represent the energy, momentum, and charge of these particles is also variable. A common solution to this problem is to fix a maximum number of particles and then truncate and zeropad to fit the data into a fixed tabular representation, though this is not the natural representation of the data and it leads to a loss of information.

Sequential data is also commonly encountered in physics (e.g. in time series). Here an individual entry $x_i = (x_i^1, \dots, x_i^t, \dots, x_i^{T_i})$ where $t$ is index for the ordered sequence, $T_i$ is the length of the sequence (which might be variable), and the data associated to each “time” $x_i^t \in \mathbb{R}^d$. This is similar to the previous example where the energy, momentum, and charge of the $t$-th electron in the $i$-the event would be $x_i^t$ and the electrons might be sorted according to their energy or transverse momentum. Sequential data is also encountered in natural language processing, where $x_i^t$ correspond to individual words in a sentence. Recurrent neural networks (see

Sec. 41.5.3.9) are particularly well suited to sequential data. Examples applications from the Living Review include [80–85].

Image-like data is one of the most dominant forms of data in industrial applications of deep learning, is very relevant for astronomy and cosmology, and also appears in particle physics in various forms. Image-like data typically involves $d$-dimensional features associated to a regular grid or lattice that does not vary across the individual instances $x_i$. The canonical example is a standard image from a camera with $W \times H$ pixels where the $p$-th pixel has data $x_i^p \in \mathbb{R}^3$ corresponding to the three channels in the RGB color model. It is important to recognize that the data corresponding to the 2-dimensional image is not 2-dimensional; instead, it is $(W \times H \times c)$-dimensional, where $c$ is the number of channels. In astronomy, an image may be grey scale ($c = 1$) or there may be more channels ($c > 3$) corresponding to different color filters. In other applications, the grid or lattice might be 3- or 4-dimensional. For example, the data associated to a regularly segmented calorimeter can be thought of as a 3-dimensional image and the data associated to a lattice simulation of a classical or quantum system can be thought of as a 4-dimensional image. Convolutional neural networks, described in Sec. 41.5.3.5, are particularly well suited to image-like data. Example applications from the Living Review include [81, 86–108].

It is also possible that the data (or features) associated to one “pixel” or lattice site may itself be structured. For example, the single read-out plane of a Liquid Argon time projection chamber (LArTPC) may involve a 1-dimensional or 2-dimensional grid, but the data associated to each “pixel” is itself a sequence or waveform. Example applications in Neutrino Physics from the Living Review include [109–139]. Similarly, in lattice quantum chromodynamics, the data associate to a particular site (or link) would be group valued (e.g. $x_i^p \in SU(3)$ as in Refs. [140,141]).

Both sequential and image-like data have a notion of temporal or spatial structure. While it is possible to unroll an image into a $(W \times H \times c)$-dimensional vector, that would erase the spatial structure and obfuscate the fact that nearby pixels are highly correlated. Similarly, one could permute the time index for sequential data, but that would destroy the temporal structure of the data. The complementary point of view is that the model class should also be aware of the structure of the data. Concurrent neural networks and convolutional neural networks are good examples of what is called inductive bias as the models do incorporate something about the structure of the data that they are expected to be used with. In some cases this can be formalized in terms of symmetry. For example, if we train model to classify images of cats and dogs, we would like it’s prediction to be invariant to where in the image the cat is. This type of translational invariance can be enforced in the design of the model class.

While permuting the elements of a sequence destroys the temporal structure of a time series, attaching a temporal index $t$ to a set of objects with features $x_i^t$ can also be problematic. If the data corresponding to $x_i$ really are a set $\{x_i^1, \dots, x_i^{T_i}\}$ (e.g. a point cloud), then we would like the output of the model to be permutation invariance. A standard sequential or convolutional model will not generally be permutation invariant, but recently models such as Deep Sets, various types of graph neural networks, and transformers can be made to enforce permutation symmetry. Example applications from the Living Review include [100, 142–152].

The temporal and spatial structure of sequences and image like data can also be generalized. For instance, a 1-dimensional sequence can be generalized to a tree structured data like one finds in the hierarchical clustering of jets or as in a directed-acyclic graph (DAG). Generalizations of recurrent neural networks have been constructed that can operate over these more complex data structures [153, 154]. More generally, one can considered graph-structured data composed of nodes and edges or multi-graphs that group together three nodes into faces or $k$ nodes into $k$-edges. Graph neural networks are a class of models that work with this type of data. The emerging subfield of geometric deep learning aims to unify the notation, terminology, and theory that connect these considerations of the structure of the data and the corresponding model architecture. Example applications in the Living Review include [121, 128, 131, 132, 134, 155–181].
If the data are expected to have a symmetry associated to them but one is working with a model class that does not enforce this symmetry, then data augmentation is a common procedure used to improve generalization performance. Here one starts with an initial dataset {$x_i$}$_{i=1,...,n}$ and produces an augmented dataset {$x'_i$}$_{i'=1,...,n'}$ through some data augmentation strategy. For example, one might apply a random rotation $R_{i'}$ to an image to produce $x'_i = R_{i'}(x_i)$ if one assumes rotational invariance in the underlying problem.

In some cases some of the individual features (components) of $x$ are functions of other features. For instance, one may include components of a four-vector ($E, p_x, p_y, p_z$) as well as redundant information such as transverse momentum, azimuthal angles, rapidity, etc. In this case, the data is restricted to a lower-dimensional surface embedded in $\mathcal{X}$. Even if the features aren't redundant, statistically the data are often effectively restricted to a small subspace of statistically likely samples and those that are exceedingly unlikely. For instance, the space of natural images is a small and highly structured subspace of all possible images, which are dominated by what we would perceive visually as noise. The term data manifold is used to describe this restricted subspace where the data are to be found, even though it does not necessarily satisfy the formal requirements of a manifold in the mathematical sense.

These considerations on the structure of the data not only apply not to the input data $x_i \in \mathcal{X}$, but also to the output data $y_i \in \mathcal{Y}$. For instance, one might want a sequence-to-sequence model as in machine translation of written text [182] or to learn a function that takes sets as input and produces graphs as output as in the Set2Graph mode [183]. One might also want the input and output of the model to be different in representations of an underlying symmetry group and for the model to enforce group-equivariance [140, 141]. The development of the necessary modelling components to enable practitioners to compose and train these types of models is a significant development for the field of physics.

## 41.5 Flavors of ML models

### 41.5.1 Support vector machines and kernel machines

Support Vector Machines (SVM) are a class of supervised learning models used for classification and regression. The learning algorithm involves a convex optimization problem that has a unique solution and can be solved with quadratic programming techniques. In this sense, they are robust and easier to characterize than neural networks that involve non-convex optimization. While their history goes back to Vapnik and Chervonenkis in the 1960s, they gained popularity after Bernhard Boser, Isabelle Guyon and Vladimir Vapnik suggested a way to create nonlinear classifiers by applying the kernel trick in 1992 [184]. Originally they were developed for binary classification and only supported the restricted case where the training dataset can be separated without errors, but in 1995 Cortes and Vapnik extended the technique to non-separable training dataset [185]. A variant targeting regression, known as support vector regression (SVR), was developed in 1997 [186]. The theory around support vector machines was well developed, and they dominated machine learning research for roughly a decade before the rise of deep learning (see Sec. 41.5.3.4).

**Maximum-margin classifiers** Linear support vector machines are used for binary classification, where $\mathcal{X} = \mathbb{R}^d$ and the target labels are conventionally defined as $\mathcal{Y} = \{-1, 1\}$. The classification is simply based on which side of a hyperplane the data lie. Any hyperplane can be written as the set of points $x$ satisfying $w^T x - b = 0$, where $w, b \in \mathbb{R}^d$ are the parameters of the model. The vector $w$ is normal to the hyperplane, but not necessarily normalized. The quantity $\frac{b}{\|w\|}$ quantifies the offset of the hyperplane from the origin along the normal vector $w$.

If the training dataset is linearly separable, then there is a region bounded by two parallel hyperplanes, called the margin, that separate the two classes of data. The maximum margin classifier is uniquely defined by making the distance between these two hyperplanes as large as possible. The boundaries of the margin

can be defined by $w^T x_i - b = \pm 1$, and the width of the margin is given by $\frac{2}{\|w\|}$. Figure 41.1 illustrates this for $x \in \mathbb{R}^2$.

Since the width of the margin is maximized when $\|w\|$ is minimized, we can state the goal of the (hard) maximum-margin classifier in the linear separable case as the following constrained optimization problem: Minimize $\|w\|^2$ subject to the constraint $y_i(w^T x_i - b) \ge 1$ for $i = 1, ..., n$. The $w$ and $b$ that solve this problem uniquely determine the resulting classifier, $\hat{y}(x) = \text{sgn}(w^T x - b)$.

This geometric description makes it clear that the maximum-margin hyperplane is completely determined by those $x_i$ that lie nearest to it: the eponymous support vectors.

**Soft margins and slack variables** Often the data will not be linearly separable and thus the hard-margin optimization problem described above has no solution. Such a non-separable case can be accommodated with the use of so-called slack variables. One slack variable $\xi_i > 0$ is introduced for each data point. As shown in Figure 41.1 (right), if a data point $x_i$ is outside the margin and correctly classified as in the separable case, then $\xi_i = 0$. If it within the margin but still correctly classified, then $0 < \xi_i < 1$. If it is misclassified, then $\xi_i > 1$. The new constrained optimization problem is to minimize

$$ \sum_i \xi_i + \lambda \|w\|^2 \quad (41.32) $$

subject to the constraints $y_i(w^T x_i - b) \ge 1 - \xi_i$ and $\xi_i > 0$ for $i = 1, ..., n$. Effectively, this corresponds to minimizing $\sum_i \max(0, 1 - y_i(w^T x_i - b)) + \lambda \|w\|^2$, where the first term is called the hinge loss.

**The dual problem** In the language of convex optimization, linear programming, and quadratic programming the optimization problems stated above for the hard and soft maximum-margin classifiers define the primal objective to be minimized. There is a Lagrange dual problem where the solution corresponds to maximization with $\min_{\phi} L(\phi) = \sup_{\alpha > 0} \bar{L}(\phi, \alpha)$, where $\alpha$ are Lagrange multipliers. While it is not apparent in the primal formulation of the constrained optimization problems, the solution to the dual problem takes on the form

$$ \hat{y}(x) = \text{sgn}\left(\sum_i y_i \alpha_i x_i^T x - b\right), \quad (41.33) $$

where the $\alpha_i > 0$ for the support vectors and otherwise zero. This reformulation is plausible since the orientation of the maximum-margin hyperplanes defined by $w$ is completely specified by the support vectors. The crucial observation for the kernel trick described below is that the final model is defined in terms of coefficients and the inner product $\langle x_i, x \rangle_{\mathcal{X}} = x_i^T x$.

**The kernel trick** One approach to obtain non-linear decision boundaries is to introduce a non-linear mapping $\phi: \mathcal{X} \to \mathcal{V}$ that embeds the data in some higher-dimensional space $\mathcal{V}$. One can then use the linear SVM described above in the space $\mathcal{V}$. This is not unlike the many neural network models where the first layers use a non-linear activation function while the last layer is linear.

The kernel trick avoids the explicit mapping $\phi$ that is needed to get linear learning algorithms to learn a nonlinear function or decision boundary. This is possible because the solution to the dual problem in Eq. 41.33 is represented in terms of inner products $\langle \phi(x_i), \phi(x) \rangle_{\mathcal{V}} = \phi(x_i)^T \phi(x)$ corresponding to the candidate point $x$ and the support vectors $x_i$. Thus all that is needed is a way to evaluate the inner products. A kernel or a kernel function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ is a function that can be expressed as an inner product in another space $\mathcal{V}$ via $k(x, x') = \langle \phi(x), \phi(x') \rangle_{\mathcal{V}}$. Critically, one does not need to ever evaluate the map $\phi(x')$, which is implicitly defined by the kernel. In some cases the kernel is straight forward to evaluate, even though the target space $\mathcal{V}$ implicitly defined by the kernel would be infinite dimensional.

The choice of the kernel is roughly the analog of the architecture of a neural network. It dictates the types of non-linearities the
<figure><img src="image_263.png" /><figcaption>**Figure 41.1:** Illustration of a maximum margin classifier for a linear support vector machine in the separable case (left) and with the soft-margin and slack variables $\xi$ (right).</figcaption></figure>

SVM can implement in the input space. It also provides an opportunity to encode prior knowledge from physics such as symmetries, length scales, etc. There are also connections between very wide neural networks with random weights and biases and the kernel of a Gaussian process [187–189], which, when applied to deep networks, lead to the concept of neural tangent kernel [190]. While it is beyond the scope of this chapter, the interested reader may also be curious to understand the connections as formalized in the language of Reproducing kernel Hilbert spaces (RKHS) [191].

**Support vector regression** There are a few flavors of SVR. In Vapnik's $\epsilon$-SVR formulation is similar to the picture for classification, but the boundary hyperplane serves as the prediction and the width of the margin for which there is no penalty is $\epsilon$. Again, one can introduce a non-linear map $\phi(x)$ explicitly or implicitly via the corresponding kernel $k(x, x')$. One can view the predictive model as $\hat{y}(x) = w^T\phi(x) - b$ and the goal of the optimization problem is to minimize

$$ \|w\|^2 + C \sum_i (\xi_i + \xi_i^*) \qquad (41.34) $$

where $\xi_i$ and $\xi_i^*$ are slack variables associated to the model over or under-estimating $y_i \in \mathbb{R}$. The constrained optimization is subject to $w^T x_i - b - y_i \le \epsilon + \xi_i$, $y_i - w^T x_i + b \le \epsilon + \xi_i^*$, $\xi_i, \xi_i^* \ge 0$ for $i = 1, \dots, n$. The solution to the dual problem leads to a kernelized representation of the solution

$$ \hat{y}(x) = \sum_i (\alpha^* - \alpha) K(x, x_i) - b, \qquad (41.35) $$

where the $\alpha_i, \alpha_i^*$ are the parameters of the dual problem. In libraries like LibSVM the output of the algorithm is just $\alpha_i^* - \alpha_i$ for $i = 1, \dots, n$.

**Kernel ridge regression** Kernel Ridge Regression (KRR) uses the same model as in SVR; however, it uses a different loss function. KRR uses the squared error loss while support vector regression uses the $\epsilon$-insensitive loss function (i.e. no penalty is enforced if the prediction is within $\epsilon$). Both are combined with l2 regularization (the $\|w\|^2$ term). While SVR leads to a sparse solution that only depends on the support vectors, it can be slow to solve the constrained convex optimization problem. In contrast, KRR can be done in closed-form and is typically faster to train for medium-sized datasets. The consequence is that the learned model is non-sparse and thus slower than SVR to evaluate at prediction-time.

**Gaussian Process Regression (kriging)** Another variation on kernel-based regression is Gaussian process regression (GPR), also referred to as krigging. Again the model is based on the kernel trick. While KRR aims to minimize the mean-squared error loss $(y_i - \hat{y}(x_i))^2$ along with an l2 regularization term proportional to $\|w\|^2$, GPR elevates these terms to a probabilistic interpretation where the first is proportional to a log-likelihood for Gaussian $N(y_i|\hat{y}(x_i)), \sigma)$ likelihood and the second corresponds to the log of a Gaussian prior on $w_i$. This is related to the discussion on Tikhonov regularization found in Sec. 41.3.1.2. GPR goes further and uses a mean function $m(x)$ and kernel $k(x, x')$ to define a probabilistic model for $y(x)$ and $y(x')$ that includes correlations. For any finite number of points $\{x_j\}_{j=1,\dots,m}$, the Gaussian process defines the joint distribution for the values that $f(x_i)$ might take. Typically, one specifies the prior mean $m(x)$ and prior covariance function $k(x, x')$ and then conditions on the training dataset $\{x_i, y_i\}_{i=1,\dots,n}$. The posterior can then be evaluated for an independent set of points $\{x_j\}_{j=1,\dots,m}$. The posterior mean is typically used as the prediction equivalent to SVR and KRR, but GPR also gives a meaningful notion of uncertainty in the predictions, and can be viewed as a realization of kernel Bayesian regression [192].

One advantage of GPR is that one can work with a family of kernel functions parameterized by some hyperparameters $\eta$. One can then optimize the hyperparameters via gradient-ascent on the marginal likelihood function. In contrast, hyperparameter tuning in SVR and KRR typically requires a grid search or some other black-box optimization procedure evaluated on held out data or some form of cross-validation.

Rasmussen and Williams provides an excellent review of Gaussian processes [193]. Numerically, Gaussian Process libraries are confronted with computing the log determinant of the covariance kernel, which naively scales like $O(n^3)$ in computational complexity. However, recently there have been a number of numerical advances that make these methods fast and scalable to large datasets [194,195]. While Gaussian Processes are used more widely in cosmology, they are not widely used in particle physics or astro-particle physics. Recent works explore the design of physics-inspired kernels and use Gaussian Processes to model the intensity for a Poisson point process like those found in experimental particle physics and $\gamma$-ray and X-ray astronomy [196–198]. Gaussian Processes are also extensively used in Bayesian Optimization (Section 41.3.3).

### 41.5.2 Decision trees

**Tree-based models** Classification and Regression Trees (CART) typically partition the input space into $J$ disjoint regions $\mathcal{X} = \mathcal{X}^1 \cup \dots \cup \mathcal{X}^J$ through a sequence of $J-1$ binary splits based
<figure><img src="image_264.png" /><figcaption>**Figure 41.2:** Figure depicting the nonlinear SVM decision boundary in $\mathcal{X}$ (left) and the corresponding linear maximum margin classifier in $\mathcal{V}$ (right). Figure adapted from Kernel_Macine.svg Wikimedia Commons (CC BY-SA 4.0) by Cranmer.</figcaption></figure>

on an individual components of $x \in \mathcal{X}$ (e.g. $x_4 < 0.7$) [199]. The model is piecewise constant and assigns the value $b_j \in \mathcal{Y}$ to the j-th terminal region $\mathcal{X}^j$. The model can be written

$$ \hat{y}(x) = f_{\phi}(x) = \sum_{j} b_j \mathbf{1}(x \in \mathcal{X}^j). \quad (41.36) $$

The parameters $\phi$ of the model comprise the components index and thresholds for the successive splittings and the coefficients $b_j$.

Tree learning refers to the algorithm used to choosing the tree structure and determining the predictions at leaf nodes. Optimization of the tree structure involves a difficult discrete optimization since the change in the loss with respect to the tree structure is non-differentiable and it is intractable to explore the combinatorially large space of possible trees with brute force. Therefore, the discrete optimization component of tree learning typically involves some approximate algorithm based on heuristics. In contrast, optimization of the $b_j$ for a given tree structure can exploit gradient-based optimization algorithms.

Common approaches to building the decision tree start with a root node and grow with splits based on individual attributes (components of $x$). These are referred to as top-down induction strategies. There are various impurity heuristics used for choosing the best attribute to split on such as the Gini index, cross-entropy and mis-classification error. Generally they aim to find a split that will refine the terminal nodes such that they have higher purity than the parent node.

Because most tree learning algorithms consider splits aligned with individual feature components, there are some failure modes for tree-based models. However, tree-based models work well with tabular data that is composed of a mix of continuous and discrete features. Tools such as XGBoost [200] are widely used in data science competitions such as Kaggle and the boosted decision trees (BDTs) implemented in StatPatternRecognition [201] and TMVA [77] have been one of the most used techniques in particle physics [3].

Individual trees are often referred to as weak-learners and they can be combined in various ways described below. Regularization is also an important consideration with tree-based models as one can always learn a tree that assigns exactly one training dataset point per terminal node and memorize the training dataset exactly. One approach to this is called pre-pruning, which simply terminates the growing of the trees if the number of training samples reaching the terminal node drops below some threshold, the purity of a terminal nodes is below some threshold, or if the improvement in purity due to a proposed split is not above a

threshold. Another regularization approach is called post-pruning, which uses a validation data set that is disjoint from the training dataset to probe generalization performance. In this approach, after initially growing a tree with the training dataset, one considers a sequences of pruned trees where splits are removed based on some heuristic. One finds the tree in this sequence of pruned trees that minimizes the generalization error on the validation set. Alternatively, in tools such as XGBoost there is an explicit regularization term included in the loss function (see Eq. 41.43).

**Ensemble methods** The idea of ensemble methods is to combine multiple models into a more performant one by exploiting the bias-variance tradeoff [202]. This is most commonly achieved through averaging (e.g. bagging and random forests), which primarily reduces variance, or boosting (e.g. AdaBoost and gradient boosting), which primarily reduces bias. Here bias refers to the difference between the Bayes optimal model and the average model produced by the learning procedure with different training sets and variance quantifies how much the learned model varies from one training set to another.

The motivation of boosting is to combine the outputs of many “weak” models to produce a more expressive model. Compared to averaging techniques like bagging and random forests, the model is built sequentially on modified versions of the data and the final predictions are combined through a weighted sum

$$ \hat{y}(x) = \sum_{t=1}^{T} \beta_t \hat{y}_t(x), \quad (41.37) $$

where $\beta_t$ expand the parameters of the model $\phi$.

**Bagging** The idea behind bagging (bootstrap aggregating) is to create $T$ bootstrap training datasets $B_1, \dots, B_T$ drawn from the training dataset $\{x_i, y_i\}_{i=1,\dots,n}$, then learn a model $\hat{y}_t$ for each, and finally construct an average model $\hat{y}(x) = (1/T) \sum_t \hat{y}_t(x)$. If one had $T$ independent training datasets each of size $n$, then the bias of the average model would be the same as the original model, but the variance would be reduced by a factor of $T$. By using bootstrap resampling, the bias may increase but the reduction in variance often dominates, which leads to improved performance.

**Random forests** Random forests refers to a type of “perturb and combine algorithm” that combines bagging and random attribute subset selection. Again one builds trees $\hat{y}_t(x)$ from bootstrap training datasets $B_t$, but instead of choosing the best split
among all attributes, one select the best split among a random subset of k attributes. If k includes all attributes, then it is equivalent to bagging.

**AdaBoost** In AdaBoost (adaptive boost) the sequence of trees $\hat{y}_1, \dots, \hat{y}_T$ are trained with reweighted versions of the original training dataset such that the weight of individual training sample is based on the prediction error in the previous iteration [203]. This requires working with a loss function that and learning procedure for the individual iterations that is amenable to weighted training dataset $\{x_i, y_i, w_i\}_{i=1,\dots,n}$. Incorporating the weights $w_i$ is straight forward when the risk is expressed as an expectation, since the empirical risk of Eq. 41.3 is just replaced with the weighted average. Similarly, the heuristic for many of the tree-based learning algorithms (e.g. the Gini index) also have natural generalizations with weighted events.

In the context of classification, the weighted error of the model $\hat{y}_t(x)$ is

$$ \text{err}_t = \frac{\sum_i w_i^{(t)} [1[y_i \neq \hat{y}_t(x_i)]]}{\sum_i w_i^{(t)}}. \quad (41.38) $$

Based on this weighted error, the coefficient $\beta_t$ of the component $\hat{y}_t(x)$ in Eq. 41.37 is given by

$$ \beta_t = \log \left( \frac{1 - \text{err}_t}{\text{err}_t} \right). \quad (41.39) $$

Then for the next iteration the weights of the misclassified events are updated as $w^{(t+1)} = w^{(t)} \exp(\beta_t)$ and then renormalized so that the sum of all weights is 1. This reweighted dataset is then used to train the next model $\hat{y}_{t+1}(x)$ and the entire procedure is initialized with uniform weights $w_i^{t=0} = 1/n$.

There is an analogous procedure for regression with the squared loss function based on the residuals $r_i = y_i - \hat{y}_t(x_i)$ (see for example Ref. [77] for details).

**Gradient boosting** One of the most powerful forms of tree based models, which is implemented in the tool XGBoost is referred to as *gradient boosting* [204]. In this setup, the model is purely additive as in the case of random forests, so the model is Eq. 41.37 with all $\beta_t = 1$. Note this is without loss of generality since the $\beta_t$ can be absorbed into the $b_j$ of Eq. 41.36. As with AdaBoost, the model is built sequentially through the sequence $\hat{y}_1, \dots, \hat{y}_T$.

At each iteration, a new term $f_t$ will be added to the sum in Eq. 41.37. For a given decision tree defined by splits on attributes, one can approximate the objective function (the loss function $\mathcal{L}$ plus a regularization term $\Omega$) as a function of $b_j$ in a second order Taylor series:

$$ \text{obj}^{(t)} = \sum_{i=1}^{n} [\mathcal{L}(y_i, \hat{y}_i^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i)] + \Omega(f_t) + \text{constant}, \quad (41.40) $$

where

$$ g_i = \partial_{\hat{y}_i^{(t-1)}} \mathcal{L}(y_i, \hat{y}_i^{(t-1)}) \quad (41.41) $$

and

$$ h_i = \partial_{\hat{y}_i^{(t-1)}}^2 \mathcal{L}(y_i, \hat{y}_i^{(t-1)}) \quad (41.42) $$

In XGBoost, the regularization term is taken to be

$$ \Omega(f) = \gamma J + \frac{1}{2} \lambda \sum_{j=1}^{J} b_j^2, \quad (41.43) $$

where J is the number of terminal nodes in the tree. With the second-order approximation of the objective, one can directly solve for the optimal $b_j$ for the next tree and the corresponding value of the optimized objective function. The improvement in the objective function can then be used as a heuristic for choosing the best split. Specifically, define $G_j = \sum_{i \in I_j} g_i$ and

$H_j = \sum_{i \in I_j} h_i$, where $I_j$ is the set of indices of data points assigned to the j-th leaf. The heuristic used in XGBoost for splitting a node is

$$ \text{Gain} = \frac{1}{2} \left[ \frac{G_L^2}{H_L + \lambda} + \frac{G_R^2}{H_R + \lambda} - \frac{(G_L + G_R)^2}{H_L + H_R + \lambda} \right] - \gamma. \quad (41.44) $$

This formula can be interpreted as the score on the new left leaf plus the score on the new right leaf minus the score on the original leaf minus a regularization penalty on the additional leaf. If the gain from splitting a leaf is smaller than $\gamma$, then the total Gain is negative and the split will not be added, which can be seen as implementing a form of pruning.

### 41.5.3 Neural networks

In this section we focus on the different types of components used in modern neural network architectures. Gradient-based optimization techniques are most commonly used for training neural networks, and they are described in Sec. 41.6.1. Similarly, other important aspects to effectively training neural network models such as parameter initialization, early stopping, etc. are discussed in Sec. 41.6.

The vanishing and exploding gradient problem is a common challenge for gradient-based optimization of neural networks and is described in Sec. 41.6.5. That problem is referred to repeatedly in this section because it has motivated the development of numerous architectural components described below.

#### 41.5.3.1 Feed-forward multi-layer perceptron

One of the core components in neural networks is the fully-connected, feedforward network or multi-layer perceptron (MLP), which is composed of L layers: $f = f^{(L)} \circ \dots \circ f^{(1)}$. The l-th layer defines a function that maps a $d_{l-1}$-dimensional input vector, called features, to an $d_l$-dimensional output $f^{(l)}: \mathbb{R}^{d_{l-1}} \to \mathbb{R}^{d_l}$. A unit responsible for producing an individual output in $d_l$-dimensional output is called a neuron or a filter interchangeably. For $l < L$, the functions $f_l$ are called hidden layers, and the number of neurons ($d_l$) is referred to as the width of the hidden layers. The layers in an MLP take on the form:

$$ f^{(l)}(u) = \sigma^{(l)}(W^{(l)}u + b^{(l)}), \quad (41.45) $$

where $W^{(l)} \in \mathbb{R}^{d_l \times d_{l-1}}$ is called the weight matrix, the components of the vector $b^{(l)} \in \mathbb{R}^{d_l}$ are referred to as the biases, $u \in \mathbb{R}^{d_{l-1}}$ is the input from the previous layer, $W^{(l)}u$ denotes a matrix-vector product, and $\sigma^{(l)}$ is a non-linear activation function that is usually applied element-wise. The parameters of the network comprise the full collection of weights and biases, $\phi = (W^{(1)}, \dots, W^{(L)}, b^{(1)}, \dots, b^{(L)})$.

#### 41.5.3.2 Activation functions

The activation functions $\sigma$ in neural networks are nonlinear functions and key to the expressiveness of the resulting family of functions. Two traditionally used functions are logistic and hyperbolic tangent functions. These functions are bounded to be (0, 1) and (-1, 1) respectively, and are symmetric about the input value of zero. On the other hand, away from the zero input value, a gradient of both functions quickly vanishes and this poses a challenge in using gradient-based optimization method (see Sec. 41.6.1). This can be avoided, to some extent, by normalizing the input values and carefully initializing the values of $W^{(l)}$ and $b^{(l)}$. These are discussed in Sec. 41.6.7, 41.6.8 and 41.6.9. Yet, it becomes difficult to maintain a null input value for a deep neural network, a model with many layers. Instead, a popular choice for a deep neural network is Rectified Linear Unit (ReLU):

$$ \sigma(x) = \begin{cases} x & \text{if } x > 0 \\ 0 & \text{otherwise} \end{cases} \quad (41.46) $$

which computational cost is small and provides the gradient does not vanish for $x \in (0, +\inf)$ [205,206]. An alternative to preserve a non-zero gradient in negative input values are called *Leaky ReLU* and modifies the output to $0.01x$ for $x \in (-\inf, 0]$ [207]. Another variant, called *Parametric ReLU* (PReLU), turns the coefficient
0.01 into a variable that is optimized as a part of the model during optimization [208].

The choice of activation functions depends on the model architecture and applications. As described, while the use of ReLU types are a typical choice for a deep neural network, a logistic function is a popular choice at the final layer for classification tasks. Recently, in the area of neural scene representation, sinusoidal activation functions have been found to be surprisingly effective [209].

### 41.5.3.3 Softmax

A softmax function is often used to normalize elements of a discrete vector $u$, or to interpret the output as a probability over a set of $n$ discrete categories. Given a real-valued input vector $u \in \mathbb{R}^n$, the softmax function computes the output vector $v \in \mathbb{R}^n$ the $i$-th component is given by:

$$v_i = \frac{\exp(u_i)}{\sum_{j=1}^{n} \exp(u_j)} \quad (41.47)$$

The result has the property that $v_i \in (0, 1)$ and $\sum v_i = 1$. The components of the input vector $u$ are often referred to as *logits* in reference to their connection to the logistic function used in logistic regression. The softmax function is commonly used as the last layer in multi-class classifier. The softmax is also used in the context of attention (see Sec. 41.5.3.11).

### 41.5.3.4 The rise of deep learning

There are a number of universal approximation theorems in the theory of neural networks. One of the first was that even with one hidden layer ($L = 2$), an MLP can approximate any continuous function if the non-linear activation function $\sigma$ not a polynomial and the width $d_1$ is large enough [210]. However, it is often more efficient (in terms of the number of parameters) to increase the depth of the network $L$ [211].

Training a deep network (i.e. $L > 2$) that generalizes well can be more difficult, requiring large training datasets, many gradient updates, and suitable regularization. The introduction of large labeled training sets, advances in computing (e.g. graphic processing units or GPUs which enabled orders of magnitude acceleration in parallel computation including matrix multiplies [212]), development of ReLU, research progress in initialization and optimization algorithms for model parameters, and regularization techniques like *dropout* [13] all played an important role in the rise of *deep learning* [2,213]. Though the name deep learning was originally a reference to the depth $L$ of such networks, modern deep learning is characterized more by the composition of various types of modules that are trained through gradient-based optimization. Below we introduce some other common network architectures.

### 41.5.3.5 Convolutional neural networks

Convolutional Neural Networks (CNNs) are widely used for image-like data. They implement the convolution of the input image $u$ and a *filter* $W$ (also referred to as a kernel). The parameters of the filter are learnable and the convolution involves traversing over input and calculating the inner product of the filter $W$ with the part of the input in the *receptive field*, which has the same spatial shape as the filter and is centered at the target pixel. At each location - indexed by $i$ and $j$ below - there is a pixel that may have a vector of features associated with it. In the context of CNNs, these components of these features - indexed by $c$ and $c'$ below - are often referred to as *channels* in reference to the red, green, and blue color channels in a traditional image. The convolution operation is often denoted with a $*$, and the result can be expressed as

$$v_c(j) = (W * u_c)(j) = \sum_{c'} \sum_{i} W_{c,c'}(i) u_{c'}(j - i), \quad (41.48)$$

where "j - i" is shorthand for the pixel index corresponding to the translation from pixel j to i. By repeating the operation over all pixels, the result of a kernel convolution is also an image as

illustrated in Figure 41.3. Note that the number of channels in the output $v$ does not need to be the same as in the input, and the collection of filters $W_{c,c'}$ is often referred to as a filter bank. The entire image for a fixed channel index is often referred to as a *feature map*.

A key feature of the CNN architecture is that it is *equivariant* to translations, meaning that if the input image is shifted (e.g. $u(i) \to u'(i) = u(i - k)$), then the output is also shifted by the same amount (e.g. $v(j) \to v'(j) = v(j - k)$). This equivariance property is a natural consequence of using convolutions. A fully connected MLP would not generally have this symmetry; however, it is enlightening to imagine transferring the computation performed by a CNN to the weights and biases of a fully connected MLP, which would result in duplicating the weights of the filters multiple times. In this view, the CNN can be interpreted as a fully connected MLP with *shared weights*, which would maintain the equivariance property. This view is helpful for gaining intuition about the inductive bias of models and makes clear that a CNN is a subset of the fully connected MLPs that satisfy the translation equivariance property.

The discussion above makes it clear that CNNs will have fewer parameters than the corresponding fully connected MLP, which can alleviate the optimization challenge. In addition, the convolutional structure allows for data to be reused effectively as patterns in one part of an image in the training dataset effectively contribute to learning that pattern anywhere in the image. Larger kernel sizes allow for the filters to learn more complicated patterns, but at the expense of having more parameters and needing more data to train (the extreme case being a kernel size that is the size of the entire image, which would be equivalent to a fully connected MLP). In practice, kernel sizes of 3 are popular, 5 is sometimes used, and larger kernels are rarely used.

A kernel convolution involves three hyperparameters, the size of kernel (typically an odd number so that the filter has an unambiguous center), *stride*, and *padding*. The stride is the number of pixels between each target pixel (i.e. the center of the filter). The stride size of 1 implies the target pixels are adjacent, and a stride of 2 implies skipping 1 pixel in between (along each spatial axis). Padding is an operation to expand the input image by a specified number of pixels (i.e. "padding size") for when the target pixel is near the edge and the filter would extend beyond the input image. In Figure 41.3, the input image is 4x4 pixels, the kernel is of size 3, and no (zero) padding is applied. The result is a smaller image. Alternatively, With the padding of size 1, the output would retain the 4x4 shape of the input image.

A kernel size of 1 is also frequently used and is referred to as a *1x1 convolution*. While it cannot capture a geometrical features, it can perform linear operations on the input features, including increasing or decreasing the number of features. When combined with a non-linearity and stacked repeatedly this can be seen as a small MLP attached to an individual pixel, and for this reason the 1x1 convolution is sometimes referred to as a *network-in-network* [214]. This is often used to extract more powerful features to be used by the latter layers, and also to compress features when the next layers may be computationally demanding such as a block of many convolution layers with a large kernel size [215, 216].

One may wonder how CNNs identify features with a spatial size larger than a typical kernel size. One mechanism for this is by stacking multiple convolutional layers - e.g. the composition of two 3x3 kernels will lead to an effective 4x4 kernel. In addition, a typical CNN architecture uses pooling (described below), which effectively down-samples the image so that it can be processed at different resolutions. The effective receptive field in the input image may be much larger than the kernel size in this case. An alternative approach is to use an *inception module* which is designed to extract features simultaneously using kernels of different size [215].

### 41.5.3.6 Pooling

Pooling plays an important role in convolutional neural networks both practically and in terms of their mathematical properties. A pooling operation is a type of aggregation or downsampling that takes many pixels as input and produces one pixel for output. Typically, the pooling operation is applied indepen-