{'fraction_non_alphanumeric': 0.03259564481798832, 'fraction_numerical': 0.012866342847786934, 'mean_word_length': 4.685153583617748, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Empirical Bayes unfolding of elementary particle spectra at the Large Hadron Collider

We consider the so-called unfolding problem in experimental high energy physics, where the goal is to estimate the true spectrum of elementary particles given observations distorted by measurement error due to the limited resolution of a particle detector. This an important statistical inverse problem arising in the analysis of data at the Large Hadron Collider at CERN. Mathematically, the problem is formalized as one of estimating the intensity function of an indirectly observed Poisson point process. Particle physicists are particularly keen on unfolding methods that feature a principled way of choosing the regularization strength and allow for the quantification of the uncertainty inherent in the solution. Though there are many approaches that have been considered by experimental physicists, it can be argued that few -- if any -- of these deal with these two key issues in a satisfactory manner. In this paper, we propose to attack the unfolding problem within the framework of empirical Bayes estimation: we consider Bayes estimators of the coefficients of a basis expansion of the unknown intensity, using a regularizing prior; and employ a Monte Carlo expectation-maximization algorithm to find the marginal maximum likelihood estimate of the hyperparameter controlling the strength of the regularization. Due to the data-driven choice of the hyperparameter, credible intervals derived using the empirical Bayes posterior lose their subjective Bayesian interpretation. Since the properties and meaning of such intervals are poorly understood, we explore instead the use of bootstrap resampling for constructing purely frequentist confidence bands for the true intensity. The performance of the proposed methodology is demonstrated using both simulations and real data from the Large Hadron Collider.


Introduction
This paper studies a generalized linear inverse problem (Bochkina, 2013), called the unfolding problem (Prosper and Lyons, 2011;Cowan, 1998;, arising in the analysis of the data produced at the Large Hadron Collider (LHC) at CERN, the European Organization for Nuclear Research. The LHC is the world's largest and most powerful particle accelerator. It collides two beams of protons in order to study the properties and interactions of elementary particles produced in such collisions. The trajectories and energies of these particles are recorded using four gigantic underground particle detectors and the vast amounts of data produced by these experiments are analyzed in order to draw conclusions about fundamental laws of physics. Due to their complex structure and huge quantity, the analysis of these data poses significant statistical and computational challenges.
Physicists use the term "unfolding" to refer to correcting the distributions measured at the LHC for the limited resolution of the particle detectors. Let X be some physical quantity of interest measured in the detector. This could, e.g., be the energy, mass or production angle of a particle. Due to the noise induced by the detector, we are only able to observe a stochastically smeared or folded version Y of this quantity. As a result, the observed distribution of Y is a "blurred" version of the true, physical distribution of X and the task is to use the observed values of Y to estimate the distribution of X.
The main challenge in unfolding is the ill-posedness of the problem in the sense that a simple inversion of the forward mapping from the true space into the smeared space is unstable with respect to small perturbations of the data (Engl et al., 2000;Kaipio and Somersalo, 2005;Panaretos, 2011). As such, the trivial maximum likelihood solution of the problem often exhibits spurious high-frequency oscillations. These oscillations can be tamed by regularizing the problem which is done by taking advantage of additional a priori knowledge about plausible solutions.
An additional complication is the non-Gaussianity of the data which follows from the fact that both the true and the smeared observations are realizations of two interrelated Poisson point processes denoted by M and N , respectively. As such, unfolding is an example of a Poisson inverse problem (Antoniadis and Bigot, 2006;Reiss, 1993) where the intensity function f of the true process M is related to the intensity function g of the smeared process N via a Fredholm integral operator K; that is, g = Kf , where K represent the response of the detector. The task at hand is then to estimate and make inferences about the true intensity f given a single observation of the smeared process N . Due to the Poisson nature of the data, many standard techniques based on a Gaussian likelihood, such as Tikhonov regularization, are only approximately valid for unfolding. Furthermore, estimators properly taking into account the Poisson distribution of the observations are rarely available in a closed form making the problem computationally challenging.
At present, the unfolding methodology used in LHC data analysis is far from being well-established (Lyons, 2011). The two main approaches are the expectation-maximization (EM) algorithm with an early stopping (D'Agostini, 1995;Vardi et al., 1985;Lucy, 1974;Richardson, 1972), and a certain variant of Tikhonov regularization (Höcker and Kartvelishvili, 1996). In high energy physics (HEP) terminology, the former is called the D'Agostini iteration and the latter, somewhat misleadingly, SVD unfolding (with SVD standing for singular value decomposition). In addition, a HEPspecific heuristic, called bin-by-bin unfolding, which provably accounts for smearing effects incorrectly through a multiplicative efficiency correction, is widely used. Recently, Choudalakis (2012) proposed a Bayesian solution to the problem, but this seems to have seldom been used in practice, thus far.
The main problem with the D'Agostini iteration is that it is difficult to give a physical interpretation to the regularization imposed by early stopping of the iteration. SVD unfolding, on the other hand, ignores the Poisson nature of the observations and does not enforce the positivity of the solution. Furthermore, both of these methods suffer from not dealing with two significant issues satisfactorily: (1) the choice of the regularization strength and (2) quantification of the uncertainty in the solution. The delicate problem of choosing the regularization strength is handled in most LHC analyses using non-standard heuristics or, in the worst case scenario, by simply fixing a certain value "by hand". When quantifying the uncertainty of the unfolded spectrum, the analyses rarely attempt to take into account the uncertainty related to the choice of this regularization strength. Each year, the experimental collaborations working with LHC data publish dozens of papers using such unsatisfactory unfolding techniques. Recent examples include studies of the characteristics of jets (Chatrchyan et al., 2012d), the transverse momentum distribution of W bosons (Aad et al., 2012a) and charge asymmetry in top-quark pair production (Chatrchyan et al., 2012a), to name a few.
In this paper, we propose a novel unfolding technique aimed at addressing the above-mentioned issues within a principled framework. The main features of our method, which casts the problem as Bayesian estimation of series expansion coefficients of the intensity, subject to a regularising prior, are: • Empirical Bayes selection of the regularization parameter using a Monte Carlo expectation-maximization algorithm McClure, 1985, 1987;Saquib et al., 1998;Casella, 2001); • Frequentist uncertainty quantification, including the uncertainty of the regularization parameter, using the parametric bootstrap.
To the best of our knowledge, neither of these techniques has been previously used to solve the HEP unfolding problem. Our method also properly takes into account the Poisson distribution of the observations, enforces the positivity constraint of the unfolded spectrum and imposes a curvature penalty on the solution with a straightforward physical interpretation. The unfolding problem is closely related to image reconstruction in emission tomography (Shepp and Vardi, 1982;Vardi et al., 1985;Green, 1990) and to image deblurring in optics (Richardson, 1972) and astronomy (Lucy, 1974) -once discretized, all of these are described by a similar Poisson regression problem. There are however at least three important differences between these problems and unfolding. First, in tomography and image processing, the unknown is a two-or three-dimensional image, while in HEP unfolding one is typically interested in a one-dimensional intensity spectrum. This makes the scale of the problem at least an order of magnitude smaller enabling the use of computationally intensive statistical methods, such as the ones described in this paper. Second, uncertainty quantification of the solution is crucial in high energy physics which is rarely the case on other domains using similar models; and third, images are in principle naturally discretized using pixels, while for HEP spectra other basis expansions can be more appropriate.
Classical, well-understood techniques for choosing the regularization strength in inverse problems include the Morozov discrepancy principle (Morozov, 1966) and cross-validation (Stone, 1974). Bardsley and Goldes (2009) study these techniques in the context of Poisson inverse problems, while Veklerov and Llacer (1987) provide an alternative approach based on statistical hypothesis testing. On the contrary, empirical Bayes selection of the regularization parameter, one of the key elements of our unfolding procedure, has received relatively less attention in the literature. Among the few recent contributions, Johnstone and Silverman (2005) demonstrated the good performance of the marginal maximum likelihood estimator (MMLE) in choosing the threshold levels in wavelet smoothing for direct observations under Gaussian noise. The approach we follow bears similarities to that of Saquib et al. (1998) where the MMLE is used to select the regularization parameter in tomographic image reconstruction with Poisson data. In spite of their demonstrated good performance on many real-world datasets, empirical Bayes techniques have not become widely-used in tomography due to their high computational cost (Leahy and Qi, 2000;Green, 2012). In our case, however, the smaller scale of the problem makes the computations tractable on a modern desktop computer.
The second key element of our methodology is frequentist uncertainty quantification based on the parametric bootstrap. The use of the credible intervals of the empirical Bayes posterior would provide the most straightforward way of giving confidence statements for this problem, but due to the data-driven choice of the hyperparameter, these intervals do not enjoy the same subjective interpretation as standard Bayesian intervals. Moreover, in HEP, frequentist confidence statements are generally preferred over Bayesian uncertainty quantification (Lyons, 2013). For these reasons, we explore the use of simple, albeit computationally expensive, bootstrap resampling for constructing frequentist confidence bands for the unknown intensity. This also enables us to take into account the uncertainty regarding the choice of the regularization parameter which is usually ignored in related frequentist procedures (Berk et al., 2013;Efron, 2013). A sensible alternative to our methodology would be to use hierarchical Bayes by placing a hyperprior on the unknown regularization parameter (Kaipio and Somersalo, 2005). Such an approach would enable an automatic choice of the regularization strength along with standard Bayesian uncertainty quantification, but is dependent on the choice of the hyperprior. In effect, our proposed methodology carries over the benefits of hierarchical Bayes to the frequentist setting without the need to worry about the choice of the hyperprior.
The paper is structured as follows. Section 2 provides the necessary background on the experimental data produced at the LHC and the role of unfolding in the analysis of these data. We then formulate in Section 3 a forward model for the unfolding problem using Poisson point processes. The proposed methodology of empirical Bayes unfolding is explained in detail in Section 4 which forms the backbone of this paper. This is followed by simulation studies in Section 5 and a real-world data analysis scenario in Section 6 consisting of the unfolding of the Z boson invariant mass spectrum measured at the CMS experiment at the LHC. We then close the paper with some concluding remarks in Section 7.
2 LHC data and unfolding 2.1 Experimental data at the LHC The Large Hadron Collider is a 27 km long circular proton-proton collider located in an underground tunnel at CERN in Geneva, Switzerland. With proton-proton collisions of up to 8 TeV 1 center-of-mass energy, the LHC is the world's most powerful particle accelerator. The protons are accelerated in bunches of billions of particles and bunches moving in opposite directions are led to collide at the center of four gigantic particle detectors called AL-ICE, ATLAS, CMS and LHCb. In the current experimental configuration, these bunches collide every 50 ns at the heart of the detectors resulting in some 20 million collision events per second in each detector out of which the few hundred most interesting ones are stored for further analysis.
Out of the four detectors, ATLAS and CMS are multipurpose experiments capable of performing a large variety of physics analyses ranging from Figure 1: Illustration of the detection of particles at the CMS experiment ( Barney, 2004). Each type of a particle leaves its characteristic trace in the various subdetectors of the experiment. This enables identification of different particles as well as the measurement of their energies and trajectories. Copyright: CERN, for the benefit of the CMS Collaboration. the discovery of the Higgs boson to precision studies of quantum chromodynamics. The other two detectors, ALICE and LHCb specialize in studies of lead-ion collisions and b-hadrons, respectively. In what follows, we focus on describing the data collection and analysis in the CMS experiment, which is also the source of the data of our unfolding demonstration in Section 6, but similar principles also apply to ATLAS and, to some extent, to other high energy physics experiments.
The CMS experiment (Chatrchyan et al., 2008), an acronym for Compact Muon Solenoid, is situated in an underground cavern along the LHC ring near the village of Cessy, France. The detector, weighing a total of 12 500 tons, has a cylindrical shape with a diameter of 14.6 m and a length of 21.6 m. The construction, operation and data analysis of the experiment is conducted by an international collaboration of over 4000 scientists, engineers and technicians. When two protons collide at the center of CMS, their energy is transformed into matter in the form of new particles. A small fraction of these particles are exotic, short-lived particles, such as the Higgs boson or the top quark, which are at the center of the scientific interest of the high energy physics community. Such particles decay almost instantly into more familiar, stable particles, such as electrons, muons and photons. Using various subdetectors, the energies and trajectories of these particles are recorded in order to study the properties and interactions of the exotic particles created in the collision.
The layout of the CMS detector is illustrated in Figure 1. The detector is immersed in a 3.8 T magnetic field created using a superconducting solenoid magnet. This magnetic field bends the trajectory of any charged particle traversing the detector, and since the higher the momentum of the particle, the less it bends, this enables the measurement of its momentum. CMS consists of three layers of subdetectors: the tracker, the calorimeters and the muon detectors. The innermost detector is the silicon tracker, which consists of an inner layer of pixel detectors and an outer layer of microstrip detectors. When a charged particle passes through these semiconducting detectors, it leaves behind electron-hole pairs and hence creates an electric signal. These signals are combined into a particle track using a Kalman filter in order to reconstruct the trajectory of the particle. The next layer of detectors are the calorimeters, which are devices for measuring the energies of particles. The CMS calorimeter system is divided into an electromagnetic calorimeter (ECAL) and a hadron calorimeter (HCAL). Both of these devices are based on the same general principle: they are made of extremely dense materials with the aim of stopping the particles passing through. In the process, a portion of the energy of these particles is converted into light in a scintillating material and the amount of light, which depends on the energy of the incoming particle, is measured using photodetectors inside the calorimeters. The ECAL measures the energy of particles that interact mostly via the electromagnetic interaction, in other words, electrons, positrons and photons. The HCAL, on the other hand, measures the energies of hadrons, i.e., particles composed of quarks. These include, e.g., protons, neutrons and pions. The HCAL is also instrumental in measuring the energies of jets, i.e., collimated streams of hadrons produced by quarks and gluons, and in detecting the so-called missing transverse energy, an energy imbalance caused by non-interacting particles, such as neutrinos, escaping the detector.
The outermost layer of the CMS detector consists of muon detectors, whose task is to identify and measure the momenta of muons. Accurate detection of muons was of central importance in the design of CMS since muons provide a clean signature for many exciting physics processes. This is because there is a very low probability for other particles, with the exception of non-interacting neutrinos, to penetrate through the CMS calorimeter system. For example, the four-muon decay channel played an important role in the discovery of the Higgs boson at CMS (Chatrchyan et al., 2012b).
The information of all CMS subdetectors is combined (Chatrchyan et al., 2009) to identify the stable particles, i.e., muons, electrons, positrons, photons and various types of hadrons, produced in each collision event, see Figure 1. For example, a muon will leave a track in both the silicon tracker and the muon chamber, while a photon produces a signal in the ECAL without an associated track in the tracker. The information of these individual particles is then used to reconstruct higher-level physics objects, such as jets and missing transverse energy.

The role of unfolding in LHC data analysis
The need for unfolding arises because any quantity measured by the detectors outlined above is corrupted by stochastic noise. For example, let E be the energy of an electron hitting the CMS ECAL. Then the measured value of the energy follows to a good approximation the Gaussian distribution N (E, σ 2 (E)) where the variance satisfies (Chatrchyan et al., 2008) where S, N and C are fixed constants. The noise is not necessarily additive. Furthermore, for more sophisticated measurements, such as the ones combining information from several subdetectors or more than one particle, the distribution of the response is not usually available in a closed form. Indeed, most analyses rely on detector simulations to determine the response of their physical quantity of interest.
It should be pointed out that not all LHC physics analyses directly rely on unfolding. The common factor between the examples listed in Section 1 is that these are measurement analyses and not discovery analyses meaning that these are analyses studying in detail the properties of some already known phenomenon. In such a case, the experimental interest often lies in the detailed physical shape of some distribution for which nonparametric unfolding is the appropriate tool to use, while discovery analyses almost exclusively use parametric models in the smeared space. The importance of unfolding for discovery of new physics lies in the fact that many unfolded results are either directly or indirectly used as inputs to discovery analyses. An example of this are parton distribution functions (Forte and Watt, 2013) which quantify the internal structure of a proton. These functions are estimated via fits to unfolded spectra and are then used to derive theory predictions in various discovery analyses. They, for example, played an important role in the recent discovery of the Higgs boson (Aad et al., 2012b;Chatrchyan et al., 2012b) and are vital in further searches of new physics, such as dark matter and extra dimensions (Chatrchyan et al., 2012c).
The need to unfold the measurements usually arises for the purposes of: • Comparison of experiments with different responses: The only direct way of comparing the spectra measured in two different experiments, such as ATLAS and CMS, is to compare the unfolded measurements.
• Input to a subsequent analysis: Certain tasks, such as the estimation of parton distributions functions and fine-tuning of Monte Carlo event generators, typically require unfolded input spectra.
• Comparison with future theories: When unfolded spectra are published, theorists can directly use them to compare with any new theoretical predictions which might not have existed at the time of the original measurement. This use case is sometimes considered controversial since alternatively one could publish the response of the detector and the theorists could use it to smear their new predictions.
• Exploratory data analysis: The unfolded spectrum could reveal hidden structure in the data which is not considered in any of the existing theoretical predictions.
According to the CERN Document Server (https://cds.cern.ch/), the CMS experiment published in 2012 a total of 103 papers out of which 16 made direct use of unfolding and many more indirectly relied on unfolded results. Unfolding was most often used in studies of quantum chromodynamics (4 papers), forward physics (4) and properties of the top quark (3). Most of these results relied on the questionable bin-by-bin heuristic (8), while the EM algorithm (3) and various forms of penalization (6) were also used. We expect similar statistics to also hold for the other LHC experiments.

Problem formulation
In most situations in high energy physics, the data generation mechanism can be modelled as a Poisson point process (see, e.g. Reiss (1993)). Let E be a compact interval on R, f a non-negative function in L 2 (E) and M a discrete random measure on E. Then M is a Poisson point process on state space E with intensity function f if and only if: In other words, the number of points M (B) observed in the set B ⊂ E is Poisson distributed with mean B f (s) ds and the number of points in disjoint sets are independent random variables.
For the problem at hand, the Poisson process M represents the true, particle-level observables generated in the proton-proton collisions. The smeared, detector-level observables are represented by another Poisson process N . The process N is assumed to have a state space F , which is a compact interval on R, and a non-negative intensity function g ∈ L 2 (F ). The intensities of the two processes are related by a bounded linear operator K : In what follows, we assume K to be a Fredholm integral operator, that is, where the kernel k ∈ L 2 (F × E) is assumed to be known. The unfolding problem is then to estimate the true intensity f given a single observation of the smeared Poisson process N . This Poisson inverse problem (Antoniadis and Bigot, 2006;Reiss, 1993) is ill-posed in the sense that in virtually all practical cases the pseudoinverse K † of the forward operator K is an unbounded -and hence discontinuouslinear operator (Engl et al., 2000). This means that the naïve approach of first estimating g using, for example, a kernel density estimateĝ and then estimating f usingf = K †ĝ is unstable with respect to fluctuations ofĝ. The resulting naïve estimator has a huge variance which typically exhibits itself as large, unnatural oscillations in the estimates.
To better understand the physical meaning of the kernel k, let us consider the unfolding problem at the point level. Denoting by X i the true observables, the Poisson point process M can be written as where δ X i is the Dirac measure at X i ∈ E , the variables τ, X 1 , X 2 , . . . are independent random variables such that τ ∼ Poisson(λ(E)), and the X i are identically distributed with the probability density f (·)/λ(E), where When the particles corresponding to X i traverse the detector, the first thing that can happen is that they might not be observed at all due to the limited efficiency and acceptance of the device. Mathematically, this corresponds to thinning of the Poisson process. Let Z i ∈ {0, 1} be an indicator variable showing whether the point X i is observed (Z i = 1) or not (Z i = 0). We assume that τ, (X 1 , Z 1 ), (X 2 , Z 2 ), . . . are independent and that the pairs (X i , Z i ) are identically distributed. Then the thinned true process is given by For each observed point X * i ∈ E, the detector measures a noisy value Y i ∈ F . We assume that the smeared observations Y i are i.i.d. with probability density From this, it follows that the smeared observations Y i constitute a Poisson point process whose intensity function g is given by We hence identify that the kernel k of Equation (2) is given by 4 Empirical Bayes unfolding

Outline of the proposed methodology
In this section, we propose a novel combination of statistical methods for solving the high energy physics unfolding problem formalized in Section 3. The proposed methodology is based on the following four key ingredients: 1. Discretization of the unknown particle-level intensity using a B-spline basis expansion, that is, where B j (s), j = 1, . . . , p, are the B-spline basis functions.
3. Empirical Bayes estimation of the scale δ of the regularizing smoothness prior p(β|δ) using a Monte Carlo expectation-maximization algorithm.

Frequentist uncertainty quantification and bias correction using the parametric bootstrap.
This methodology enables a principled solution of the unfolding problem, including the choice of the regularization strength and uncertainty quantification, without having to resort to heuristics or approximations. We explain below each of these steps in detail and argue why this particular choice of techniques provides a natural framework for solving the problem at hand.

Discretization of the problem
Poisson inverse problems are almost exclusively studied in a form where for the observable process N and the unobservable process M are discretized. Usually the first step is to discretize the observable process using a histogram. In many applications this has to be done due to the discrete nature of the detector. In our case, the observations are, at least in principle, continuous, but we still carry out the discretization due to computational reasons. Indeed, in many analyses, there can be millions of observed collision events and treating each of these individually would not be computationally feasible.
In order to discretize the smeared process N , let {F i } n i=1 be a partition of the smeared space F into n ordered intervals and let y i denote the number of points falling on interval F i , that is, y i = N (F i ), i = 1, . . . , n. This can be seen as recording the observed points in a histogram with bin contents y = y 1 , . . . , y n T and is indeed the form of discretization most often employed in HEP. This discretization is convenient since it now follows from N being a Poisson process that the y i are independent and Poisson distributed with means In the true space E, there is no need to settle only for histograms. Instead, we consider a basis expansion of the true intensity f , that is, where {φ j } p i=1 is a sufficiently large dictionary of basis functions. Substituting the basis expansion of f into Equation (10), we find that the means µ i are given by where we have denoted Consequently, unfolding reduces to estimating β in the Poisson regression problem y|β ∼ Poisson(Kβ) for an ill-conditioned matrix K = (K i,j ).
Since spectra in high energy physics are typically smooth functions, splines (de Boor, 2001;Schumaker, 2007;Wahba, 1990) provide a particularly attractive way of representing the unknown intensity f . Let min E = s 0 < s 1 < s 2 < · · · < s L < s L+1 = max E be a sequence of knots in the true space E. Then an order-m spline with knots s i , i = 0, . . . , L + 1, is a piecewise polynomial whose restriction to each interval [s i , s i+1 ), i = 0, . . . , L, is an order-m polynomial (i.e., a polynomial of degree m − 1) and which has m − 2 continuous derivatives at each interior knot s i , i = 1, . . . , L. An order-m spline with L interior knots has p = L + m degrees of freedom. In this work, we use exclusively order-4 cubic splines which consist of third degree polynomials and are twice continuously differentiable. Note also that an order-1 spline gives us the histogram representation of f .
There exist various bases {φ j } p j=1 for expressing splines of arbitrary order. We use B-splines B j , j = 1, . . . , p, that is, spline basis functions of minimal local support, because of their numerical stability and conceptual simplicity. O'Sullivan (1986,1988) was among the first authors to use regularized B-spline estimators in statistical applications, with the approach later popularized by Eilers and Marx (1996). In the HEP unfolding literature, penalized maximum likelihood estimation with B-splines goes back to the work of Blobel (1985) and recent contributions using similar methodology include Dembinski and Roth (2011) and Milke et al. (2013). We use the Matlab Curve Fitting Toolbox to efficiently evaluate and perform basic operations on B-splines. These algorithms rely on the recursive use of lower-order B-spline basis functions, for details, see de Boor (2001).
The non-negativity of the intensity function f is enforced by constraining β to be in R p p}. This restricts f to be non-negative since each of the B-spline basis functions B j , j = 1, . . . , p, is non-negative.

Bayesian estimation of the spline coefficients
In contrast to most work on unfolding, we take a Bayesian approach to estimation of the spline coefficients β. That is, we estimate β using the Bayesian posterior where the likelihood is given by the Poisson regression model (14), The prior p(β|δ), which regularizes the otherwise ill-posed problem, depends on a scale parameter δ, which is analogous to the regularization parameter in the classical inverse problems literature.
We decided to use the Bayesian approach for two reasons. First, it provides a natural interpretation for the regularization via the prior density p(β|δ), which should be chosen in such a way that most of its probability mass lies in physically plausible regions of the parameter space R p + . Second, the Bayesian framework enables a principled, data-driven way of choosing the regularization strength δ using empirical Bayes estimation as explained below in Section 4.4.
In order to regularize the problem, we consider the truncated Gaussian smoothness prior where the elements of the p×p matrix Ω are given by Ω The interpretation of this prior is that the total curvature of f , characterized by f 2 2 , should be small. In other words, f should be a relatively smooth function, which is true for most intensities encountered in high energy physics. The strength of the regularization is controlled by the hyperparameter δ -the larger the value of δ, the smoother f is required to be.
The prior as defined by Equation (19) does not enforce any boundary conditions for the unknown intensity f . In this case, the matrix Ω has rank p − 2 and hence the prior is potentially improper (this depends on the orientation of the null space of Ω). Although the posterior would still be a proper probability density, the rank deficiency of Ω is undesirable since the empirical Bayes approach requires a proper prior distribution. Furthermore, without any boundary constraints, the unfolded intensity has an unnecessarily large variance near the boundaries.
To address these issues, we use Aristotelian boundary conditions (Calvetti et al., 2006), where the idea is to condition the smoothness penalty on the boundary values f (s 0 ) and f (s L+1 ) and then place additional hyperpriors for these values. Since f (s 0 ) = β 1 B 1 (s 0 ) and f (s L+1 ) = β p B p (s L+1 ), we can equivalently condition on (β 1 , β p ). As a result, the prior model becomes where p(β 2 , . . . , β p−1 |β 1 , β p , δ) ∝ exp −δβ T Ωβ . We model the boundaries using once again truncated Gaussians: where γ L , γ R > 0 are fixed constants. The full prior can then be written as where the elements of the p × p matrix Ω A are given by The augmented matrix Ω A is positive definite and hence Equation (23) defines a proper probability density.
Once the hyperparameter δ has been estimated using empirical Bayes (see Section 4.4), we plug its estimateδ into Bayes' rule (15) to obtain the empirical Bayes posterior p(β|y,δ). We then use the mean of this posterior as a point estimatorβ of the spline coefficients β, that is,β = E β|y,δ , yielding the estimatorf (s) = p j=1β j B j (s) of the unknown intensity f . Of course, in practice, the posterior p(β|y, δ) is not available in a closed form because of the intractable integral in the denominator of Bayes' rule (15). Hence, we need to resort to Markov chain Monte Carlo (MCMC) (Robert and Casella, 2004) sampling from the posterior and the posterior mean is then computed as the empirical mean of the Monte Carlo sample. Unfortunately, the most elementary MCMC samplers are not well-suited for solving the problem at hand: Gibbs sampling is not computationally tractable since the full posterior conditionals do not belong to any of the standard families of probability distributions and the Metropolis-Hastings sampler with multivariate proposals is difficult to implement since the posterior can have very different scales for different components of β.
To be able to efficiently sample from the posterior, we adopt the singlecomponent Metropolis-Hastings sampler (also known as the Metropoliswithin-Gibbs sampler) proposed by Saquib et al. (1998). Denoting β −k = β 1 , . . . , β k−1 , β k+1 , . . . , β p T , the basic idea of the sampler is to approximate the full posterior conditionals p(β k |β −k , y, δ) of the Gibbs sampler using a more tractable density . One then samples from this approximate full conditional and performs a Metropolis-Hastings acceptance step to correct for the approximation error. In our case, we take a second-order Taylor expansion of the non-quadratic part of the log full conditional resulting in a Gaussian approximation of the full conditional. When the mean of this Gaussian is non-negative, we sample from its truncation to the non-negative real line, and if the mean is negative, we replace the Gaussian tail by an exponential distribution. Further details on the MCMC sampler can be found in Section III.C of Saquib et al. (1998). During each iteration of the Monte Carlo expectation-maximization algorithm used in the empirical Bayes estimation of δ (see Section 4.4), we verify the convergence and mixing of the MCMC sampler by monitoring the acceptance rates of the Metropolis-Hastings proposals and the autocorrelation times κ j , j = 1, . . . , p, of the Markov chain. The latter measure how often the sampler on average produces an independent observation from the posterior and is estimated using Geyer's initial convex sequence estimator (ICSE) (Geyer, 1992) computed using the R package mcmc (Geyer and Johnson, 2013). The autocorrelation times κ j enable us to define the effective sample sizes ESS j = S/κ j , j = 1, . . . , p, where S is the size of the MCMC sample. ESS j measures the effective number of independent observations obtained for the jth component of the Markov chain (Kass et al., 1998, p. 99). For the MCMC iteration producing the final point estimateβ, we also monitor the trace plots, histograms, estimated autocorrelation functions and cumulative means of each component β j , j = 1, . . . , p, of the Markov chain.

Empirical Bayes selection of the regularization strength
The Bayesian approach to solving inverse problems is particularly attractive since it admits selection of the regularization strength δ using marginal maximum likelihood estimation. For a comprehensive introduction to this and related empirical Bayes methods, see, e.g., Chapter 5 of Carlin and Louis (2009). The main idea in empirical Bayes is to regard the marginal distribution p(y|δ) appearing in the denominator of Bayes' rule (15) as a parametric model for the data y and then use standard frequentist point estimation techniques to estimate the hyperparameter δ.
The marginal maximum likelihood estimator (MMLE) of the hyperparameter δ is defined as the maximizer of p(y|δ) with respect to δ. That is, we estimate δ usinĝ Computing the MMLE is non-trivial since we cannot evaluate the highdimensional integral in (25) either in a closed form or using standard numerical integration methods. Monte Carlo integration, where one samples {β (s) } S s=1 from the prior p(β|δ) and then approximates is also out of question. This is because, in the high-dimensional parameter space, most of the β (s) 's fall on regions where the likelihood p(y|β (s) ) is numerically zero. Hence we would need an enormous sample size S to get even a rough idea of the marginal likelihood p(y|δ). Luckily, it is possible to circumvent these issues by using the expectationmaximization (EM) algorithm (Dempster et al., 1977;McLachlan and Krishnan, 2008) to find the MMLE. In the context of Poisson inverse problems, this approach was originally proposed by McClure (1985, 1987) for tomographic image reconstruction and later studied and extended by Saquib et al. (1998), but has received little attention since then. When applied to the unfolding problem, the standard EM prescription reads as follows. Let (y, β) be the complete data, in which case the complete-data log-likelihood is given by l(δ; y, β) = log p(y, β|δ) = log p(y|β) + log p(β|δ), where we have used p(y, β|δ) = p(y|β)p(β|δ). In the E-step of the algorithm, one computes the expectation of the complete-data log-likelihood over the unknown spline coefficients β conditional on the observations y and the current hyperparameter δ (t) : = E log p(y, β|δ) y, δ (t) where the constant does not depend on δ. In the subsequent M-step, one maximizes the expected complete-data log-likelihood Q(δ; δ (t) ) with respect to the hyperparameter δ. This maximizer is then used as the hyperparameter on the next step of the algorithm: By Theorem 1 of Dempster et al. (1977), each step of this iteration is guaranteed to increase the incomplete-data likelihood p(y|δ), that is, p(y|δ (t+1) ) ≥ p(y|δ (t) ), t = 0, 1, 2, . . . With this construction, the incomplete-data likelihood conveniently coincides with the marginal likelihood and hence the EM algorithm enables us find the MMLEδ of the hyperparameter δ. The expectation in Equation (31), again involves an intractable integral, but can be computed using Monte Carlo integration. We simply need to sample {β (s) } S s=1 from the posterior p(β|y, δ (t) ) and then replace the expectation by its Monte Carlo approximation: The posterior sample can be obtained using the single-component Metropolis-Hastings sampler described in Section 4.3. The resulting variant of the EM algorithm is called a Monte Carlo expectation-maximization (MCEM) algorithm (Wei and Tanner, 1990). Due to the inevitable Monte Carlo error on each E-step, the MCEM algorithm loses the monotonicity property of the standard EM algorithm and theoretical analysis of its convergence becomes involved. However, in certain special cases, the iteration has been shown to eventually reach an arbitrarily small neighborhood of the maximizer with a high probability (Chan and Ledolter, 1995).
This MCEM algorithm has a rather intuitive interpretation. First, on the E-step, we use the current iterate δ (t) to produce a sample of β's from the posterior. Since this sample summarizes our current best understanding of β, we then tune the prior by varying δ on the M-step to match this sample as well as possible, and the value of δ that matches the posterior sample the best will then become the next iterate δ (t+1) .
One could also wonder why Monte Carlo integration works for the expectation of Equation (32) while it did not work for directly computing the marginal likelihood p(y|δ) in Equation (26). There are at least two reasons for this. First, in the MCEM algorithm, the β's are sampled from the posterior and hence most of them correspond to reasonable unfolded intensities. This means that they should also lie within the region where the bulk of the prior probability mass is located, thus making the sample mean in Equation (33) well-behaved. On the contrary, in Equation (26), the sample is generated from the prior resulting mostly in intensities that do not match the data very well. Second, the sum in (26) is over plain densities instead of log-densities as in Equation (33). This makes the MCEM computations considerably more robust against small probability density function values.
When p(β|δ) is given by the Aristotelian smoothness prior (23), the M-step of the MCEM algorithm is available in a closed form. Taking normalization into account, the prior density is given by where the normalization constant C(δ) depends on the hyperparameter δ and satisfies Hence where the constant does not depend on δ. Plugging this into Equation (34), we find that the maximizer on the M-step is given by The resulting iteration for finding the MMLEδ is summarized in Algorithm 1. The MCMC sampler is started from the empirical mean of the posterior sample of the previous iteration in order to facilitate the convergence of the Markov chain. In this work, we run the MCEM algorithm for a fixed number of steps T , but one could easily devise more elaborate stopping rules for the algorithm. Note, however, that the optimal choice of this stopping rule and the MCMC sample size S are, to a large extent, open problems (Booth and Hobert, 1999).

Uncertainty quantification and bias correction
The final ingredient of our procedure is uncertainty quantification and bias correction of the estimated intensityf . In contrast to most other applications of Poisson inverse problems, uncertainty quantification is of vital importance in our problem setting. It turns out that, because of our use of empirical Bayes, uncertainty quantification off is not entirely straightforward. For example, credible intervals based on the empirical Bayes posterior p(β|y,δ) lose their subjective Bayesian interpretation because of the datadriven choice of the hyperparameter δ. Also, such intervals do not take into account uncertainty regarding the choice of δ and their frequentist properties are poorly understood.
There has been a fair amount of work on correcting the naïve empirical Bayes confidence intervals (EBCI) obtained using the posterior p(β|y,δ) to account for the uncertainty ofδ (see Section 5.4 of Carlin and Louis (2009)), including the bootstrap technique of Laird and Louis (1987). This work, however, is aimed at achieving coverage with respect to the hierarchical sampling model p(y, β|δ) = p(y|β)p(β|δ), while in our case standard frequentist coverage with respect to p(y|β) would arguably be a more desirable goal. This is because in our case the prior p(β|δ) is introduced simply to regularize the ill-posedness of the problem and does not take part in the actual physical process generating the data y.
We propose quantifying the uncertainty off using a parametric bootstrap technique which is distinct from that of Laird and Louis (1987) by aiming for confidence intervals with standard frequentist coverage. The approaches we propose are similar to those of Cowling et al. (1996) but extend their results to the case of an indirectly observed Poisson point process. Our starting point is to regard the estimatorβ as a frequentist point estimator of β, that is,β =β(y) = E β|y,δ(y) . We then resample the data y and plug in the resampled observations y * to obtain the resampled estimateŝ β * =β(y * ) = E β|y * ,δ(y * ) .
In our problem setting, one can envisage several different resampling schemes to obtain the bootstrapped observations y * . In particular, we consider the following two parametric resampling procedures: ∼ Poisson(Kβ), whereβ = E β|y,δ(y) , our empirical Bayes point estimate of the spline coefficients β.
Of these, the former corresponds to Method 1 of Cowling et al. (1996) and the latter to their Method 2. Irrespective of the resampling method used, we rerun the MCEM algorithm for each y * to find the bootstrapped hyperparameterδ * =δ(y * ). By doing this, we are able to also take into account the uncertainty regarding the choice of the hyperparameter δ. The resampled spline coefficients are then found as the mean of the bootstrapped posteriorβ * = E β|y * ,δ * resulting in the bootstrapped unfolded intensitŷ f * (s) = p j=1β * j B j (s). This procedure is then repeated R times to obtain a . . . sample of bootstrapped intensities F * = {f * (r) } R r=1 . The resulting bootstrap procedure is illustrated in Figure 2.
Various techniques have been proposed for constructing confidence bands for f based on the bootstrap sample F * , see Efron and Tibshirani (1993) and Davison and Hinkley (1997). Lettingf * α (s) denote the 100 · αth percentile of the sample F * evaluated at s ∈ E, we form pointwise confidence bands for f using the following two standard techniques: Basic bootstrap interval: For every s ∈ E, an approximate 1 − 2α con- Percentile interval: For every s ∈ E, an approximate 1 − 2α confidence interval for f (s) is given by [f * α (s),f * 1−α (s)].
Choosing between resampling schemes 1 and 2 and basic and percentile intervals is tricky since there exists no clear consensus on their relative merits and superiority (Cowling et al., 1996;Davison and Hinkley, 1997;Efron and Tibshirani, 1993). Scheme 2 will produce bootstrapped estimatesf * which follow closely the actual sampling distribution off . As such, we found that E f * |y ≈f which invalidates the use of the bootstrap to recover the bias off . Furthermore, when scheme 2 is used, there is usually little difference between the basic intervals and the percentile intervals. Under scheme 1, on the other hand,f * will follow the sampling distribution off conditional on the observed value of the estimator, hence enabling the bootstrap to probe the bias off . If a large bias is present inf and scheme 1 is used, the percentile intervals will perform poorly as they will be "upside down", while the basic intervals will implicitly account for the bias. We thus recommend the combination of scheme 1 and basic intervals be used iff is suspected to be significantly biased, while for large sample sizes with small biases, the conceptually simpler combination of scheme 2 and percentile intervals can also be used. Of course, if sufficient computational resources are available, the best would be to construct both of these combinations and see if they agree. If they do, either can be used, but if there is a disagreement, then the combination of scheme 1 with basic intervals is likely to be more trustworthy due to its ability to (partially) account for the bias.
In the case of significant bias, it is also possible to apply a bootstrap bias correction to the point estimatef . The standard bootstrap estimate of the bias off at s ∈ E is bias * f (s) = 1 R R r=1f * (s) −f (s) which gives rise to the the bias-corrected point estimatef BC (s) =f (s) − bias * f (s) . Note that, given the discussion above, bias correction only makes sense when resampling scheme 1 is used. We conclude this section by noting that although using the bootstrap is computationally intensive, the computational cost can be alleviated through the use of parallel computing. Indeed, the bootstrap procedure outlined above is fully parallelizable since no communication is required between the individual bootstrap replications. We used the Matlab Parallel Computing Toolbox to parallelize all the bootstrap computations reported below and generally obtained a roughly three-fold speed-up of the computations on a quad-core desktop computer setup.

Experiment setup
We first demonstrate the empirical Bayes unfolding methodology using simulated data. The data were generated using a two-component Gaussian mixture model on top of a uniform background and smeared by convolving the particle-level intensity with a Gaussian density. Specifically, the true process M had the intensity f (s) = λ tot π 1 N (s| − 2, 1) + π 2 N (s|2, 1) + π 3 1 |E| , s ∈ E, where λ tot = E(τ ) = E f (s) ds > 0 is the expected number of true observations, |E| denotes the Lebesgue measure of E and the mixing proportions π i sum up to one and were set to π 1 = 0.2, π 2 = 0.5 and π 3 = 0.3. The true space E and the smeared space F were both taken to be the interval [−7, 7]. The true points X i were smeared with additive Gaussian noise of zero mean and unit variance. Points smeared beyond the boundaries of F were discarded from further analysis. With this setup, the smeared intensity is given by the convolution Note that this setup corresponds to the classically most difficult class of deconvolution problems since the Gaussian error has a supersmooth probability density (Meister, 2009).
The smeared space F was discretized using n = 40 histogram bins of uniform size, while the true space E was discretized using order-4 B-splines with L = 26 uniformly placed interior knots resulting in p = L + 4 = 30 unknown basis coefficients. With these choices, the condition number of the smearing matrix K was cond(K) ≈ 2.6 · 10 8 indicating that the problem is severely ill-posed. The boundary hyperparameters were set to γ L = γ R = 5. All experiments reported in this paper were implemented in Matlab and the computations were carried out on a desktop setup with a quad-core 2.7 GHz Intel Core i5 processor.

Results
We first consider a relatively easy large-sample problem where λ tot = 20000. The MCEM algorithm was started using the initial hyperparameter δ (0) = 1 · 10 −5 and was run for 20 iterations. The MCMC sampler was started from the non-negative least-squares spline fit to the smeared data, i.e., β init = min β≥0 K β − y 2 2 , where the elements ofK are given by Equation (13) with the smearing kernel k(t, s) = δ 0 (t−s). This problem is significantly less ill-posed than the unfolding problem -the condition number ofK was only 25. For each EM iteration, the single-component Metropolis-Hastings algorithm was used to obtain 500 post-burn-in observations from the posterior. After convergence of the EM algorithm, the final point estimateβ was obtained using a sample size of 1 000. The whole procedure was then repeated with R = 200 bootstrap replications obtained using resampling scheme 2. Running the MCEM iteration once to find the point estimateβ took 3 minutes, while the running time of the whole algorithm was 3 h 36 min with the bootstrap computations parallelized on the four cores of the quad-core setup. Figure 3(a) shows the true intensity f , the smeared intensity g and the unfolded intensityf with 95 % pointwise percentile intervals. The unfolded intensity beautifully captures the two peaks of the true intensity despite the severely corrupted observations. The only artifacts are the small wiggles on both tails of the intensity. Moreover, the percentile intervals cover the true f for all values of s ∈ E except for a short interval near s = 0.7. This is best seen in Figure 4(a) where the confidence intervals are shown after subtraction of the true intensity f and normalizing for the expected sample size λ tot . To enable comparison between the bootstrap and the empirical Bayes intervals, we have also plotted the naïve empirical Bayes confidence intervals in Figure 4(b). These intervals seem to cover equally well, but are consistently longer then the percentile intervals which is likely to result in overcoverage. Figure 5(a) shows the convergence of the hyperparameter estimates of the MCEM algorithm. The algorithm reduced the regularization strength from the initial value to the final estimateδ = 2.5 · 10 −7 and converged after   Figure (a) shows the unfolded intensity obtained using empirical Bayes unfolding along with 95 % pointwise percentile intervals. Figure (b) illustrates the ill-posedness of the problem by showing the unfolded intensities obtained using non-negative least-squares estimation and posterior mean estimation with a uniform prior. Also shown is the least-squares fit to the smeared data. approximately 10 iterations. During the iteration, the autocorrelation time of the MCMC sampler averaged over the components of β increased from 4.7 to 8.4 indicating that it was easier to sample from the more regularized posterior. A typical proposal acceptance rate was 98 %. For the final MCMC run producing the point estimateβ, a more careful performance analysis was made for each component of the sampler. Figure 6 shows the diagnostic plots for the components β 5 and β 21 after the removal of the burn-in. These plots indicate that the chain has converged and mixes reasonably well although the performance of the chain is typically slightly better in the interior of the space (β 21 ) than closer to the boundaries (β 5 ).
To illustrate the importance of regularization in solving this ill-posed problem, we also ran the MCMC sampler with the uniform prior p(β) ∝ 1, β ∈ R p + . In the absence of regularization, the single-component Metropolis-Hastings algorithm had significant issues exploring the parameter space. Indeed, for a sample size 1 000 (after a burn-in of 500 observations), the average autocorrelation time was 60.8 and the largest autocorrelation time 191.5 corresponding to only 5.2 effective observations from the posterior. This slow mixing was also apparent in the the trace plots and cumulative means of the chain. Unsurprisingly, the posterior mean computed based on this sample exhibits many undesired oscillations as seen in Figure 3(b). The figure also depicts the non-negative least-squares solutions corresponding to the design matrices K andK. The latter was used as the starting point β init of the MCMC iteration and is relatively well-behaved, but once one tries to undo the smearing, the solution falls apart.
In order to consider a more difficult test case, we repeated the experiment with the expected sample size λ tot = 1 000. In this case, we found that the MCEM algorithm converged more slowly and that the hyperparameter es-  Figure (a) illustrates the convergence of the Monte Carlo EM algorithm and shows that the algorithm converges faster for larger sample sizes. Figure (b) shows the convergence of the mean integrated squared error (MISE) as the expected sample size λ tot grows. Note that convergence is only obtained for MISE/λ 2 tot . The error bars indicate approximate 95 % confidence intervals, and the dotted straight line was added as a reference to illustrate that the convergence appears to be slightly slower than that given by a power law.
timates δ (t) exhibited larger Monte Carlo variation. We hence increased the number of MCEM iterations to 30 and sampled 1 000 observations from the posterior on each EM iteration. In addition, we used bootstrap resampling scheme 1 which enables us to probe the bias of the estimator. Otherwise the parameters of the experiment were the same as above. With these changes, obtaining the point estimateβ took 9 minutes and the full running time was 9 h 56 min. Figure 5(a) illustrates that the MCEM iteration increased the regularization strength toδ = 1.8 · 10 −4 and converged after approximately 20 iterations. During the iteration, the mean autocorrelation time increased from roughly 3.5 to 4.6 indicating that in this case it was slightly more difficult to sample from the more regularized posterior. The diagnostic plots for the final sampling did not indicate any problems with the convergence and mixing of the sampler. The resulting unfolded intensityf is represented by the dashed curve in Figure 7. The estimate is clearly biased near the peaks and the trough of the true intensity, but this can be mitigated with bootstrap bias correction. The bias-corrected estimatef BC , shown as the solid curve, captures the shape of the true intensity significantly better than the original estimate, but, as always, this reduction in the bias comes at the cost of increased variance visible particularly near the boundaries of the space. The figure also shows the 95 % pointwise basic bootstrap intervals which seem to cover the true intensity reasonably well, albeit potentially at the price of some slight undercoverage 2 (as suggested by Figure 8(a) where the estimates are plotted after subtracting the true intensity f and normalizing for the expected sample size λ tot ). Figure 8(b) shows also the corresponding naïve empirical Bayes confidence intervals. These intervals are longer than the basic bootstrap intervals, but as explained in Section 4.5, their statistical interpretation is unclear. Note also between Figures 4 and 8 the improvement in the point estimatef and the reduction in the length of the confidence intervals when moving from the expected sample size λ tot = 1000 to λ tot = 20 000.
To further study how empirical Bayes unfolding behaves as a function of the sample size, we repeated our first experimental setup on a logarithmic grid of expected sample sizes from λ tot = 5 000 up to λ tot = 100 000. For each sample size, we unfolded 100 independent smeared observations y and estimated the mean integrated squared error (MISE) off as the sample mean of the integrated squared errors ISE = E (f (s) − f (s)) 2 ds. As λ tot → ∞, we expect the MISE to diverge, but MISE/λ 2 tot should converge to zero, and this is indeed what we observe in Figure 5(b). On a log-log scale, the MISE estimates appear to slightly deviate from a straight line indicating that the convergence speed is likely to be close to a power law but slightly slower. 6 Unfolding of the Z boson invariant mass spectrum 6.1 Description of the data In this section, we illustrate empirical Bayes unfolding using real data from the CMS experiment at the Large Hadron Collider. In particular, we unfold the Z boson invariant mass spectrum published in Chatrchyan et al. (2013). The Z boson, which is produced in copious quantities at the LHC, is a mediator of the weak interaction. The particle is very short-lived and decays almost instantly into other elementary particles. The decay mode considered here is the decay of a Z boson into into a positron and an electron, Z → e + e − . The original purpose of these data was to calibrate and measure the resolution of the CMS electromagnetic calorimeter but they also serve as an excellent testbed for unfolding since the true intensity of this spectrum is known with remarkable precision from previous experiments. The electron and the positron produced in the decay of the Z boson are first detected in the CMS silicon tracker after which their energies E i , i = 1, 2, are measured by stopping the particles at the ECAL, see Section 2.1. From this information, one can compute the invariant mass W of the electron-positron system defined by the equation where p i , i = 1, 2, are the momenta of the two particles and the equation is written using the natural units where the speed of light c = 1. Since where m e is the rest mass of the electron, one can reconstruct the invariant mass W using only the ECAL energy deposits E i and the opening angle between the two tracks in the silicon tracker.
The invariant mass W is preserved in particle decays. Furthermore, it is invariant under Lorentz transformations and has therefore the same value in every frame of reference. This means that the invariant mass of the Z boson, which is simply its rest mass m, is equal to the invariant mass of the electron-positron system, W = m. It follows that measurement of the invariant mass spectrum of the electron-positron pair enables us to measure the mass spectrum of the Z boson itself.
Due to the time-energy uncertainty principle, the Z boson does not have a unique rest mass m. Instead, the mass follows the Cauchy distribution, also known in particle physics as the Breit-Wigner distribution, whose density is given by where m Z = 91.1876 GeV is the mode of the distribution (often simply called the mass of the Z boson) and Γ = 2.4952 GeV is the full width of the distribution at half maximum (Beringer et al., 2012). Since the contribution of background processes to the electron-positron channel near the Z peak is negligible (Chatrchyan et al., 2013), the underlying true intensity f (m) is proportional to p(m).
The dominant source of smearing in measuring the Z boson invariant mass m is the measurement of the energy deposits E i in the ECAL. The resolution of these energy deposits is in principle described by Equation (1). However, when working on a small invariant mass interval around the Z peak, it is possible to ignore the energy dependence of the resolution. Moreover, the left tail of the Gaussian resolution function is typically replaced with a more slowly decaying tail function in order to account for energy losses in the ECAL. It is therefore customary to model the smearing of the invariant mass by convolving the true intensity f (m) with the so-called Crystal Ball (CB) function (Oreglia, 1980;Chatrchyan et al., 2013) where σ, α, γ > 0 and C is a normalization constant chosen so that the function is a probability density. The Crystal Ball function is a Gaussian density with mean ∆m and variance σ 2 where the left tail is replaced with a power-law function. The parameter α controls the location of the transition from exponential decay into power-law decay and the parameter γ controls the decay rate of the power-law tail.
The dataset we use is a digitized version of the lower left hand plot of Figure 11 in Chatrchyan et al. (2013). These data correspond to an integrated luminosity 3 of 4.98 fb −1 collected at the LHC in 2011 at the 7 TeV center-of-mass energy and include 67 778 electron-positron events with the measured invariant mass between 65 GeV and 115 GeV. The data are discretized using a histogram with 100 bins of uniform width. The chosen electrons and positron have narrow particle showers in the central parts of the ECAL and as such correspond to "high quality" electron-positron pairs. For more details on the event selection, see Chatrchyan et al. (2013) and the references therein.
In order to estimate the parameters of the Crystal Ball function, we divided the dataset into two independent samples by drawing a binomial random variable independently for each bin with the number of trials equal to the observed bin contents. Consequently, the bins of the resulting two histograms are marginally mutually independent and Poisson distributed. Each observed event had a 70 % probability of belonging to the sample y used for unfolding and a 30 % probability of belonging to the sample used for CB parameter estimation.
The CB parameters (∆m, σ 2 , α, γ) were estimated using maximum likelihood with the subsampled data on the full invariant mass range 65-115 GeV. The maximum likelihood estimates were (∆m,σ 2 ,α,γ) = (0.58 GeV, (0.99 GeV) 2 , 1.81, 1.60) indicating that the measured invariant mass is on average 0.58 GeV too high and has an experimental resolution of approximately 1 GeV. As a cross-check of the fit, the estimated Crystal Ball function was used to smear the Breit-Wigner shape of the Z boson invariant mass to obtain the corresponding expected smeared histogram, which was found to be in good agreement with the observations.

Unfolding setup and results
To carry out the empirical Bayes unfolding of the Z boson invariant mass, we used the subsampled n = 30 bins on the interval F = [82.5 GeV, 97.5 GeV]. The resulting histogram y had a total of 42 475 electron-positron events.
To account for events that are smeared into the observed interval F from the outside, we let the true space E = [81.5 GeV, 98.5 GeV], that is, we extended it by approximately 1σ on both sides with respect to F . The true space E was discretized using order-4 B-splines with L = 34 uniformly placed interior knots resulting in p = 38 unknown spline coefficients. It was found out that such overparameterization with p > n facilitated the mixing of the MCMC sampler. With these choices, the condition number of the smearing matrix was cond(K) ≈ 9.0 · 10 3 . The boundary hyperparameters were set to γ L = γ R = 70. The MCEM algorithm was initialized with δ (0) = 1 · 10 −6 and was run for 20 iterations. During each MCEM iteration, the single-component Metropolis-Hastings algorithm was used to obtain 500 post-burn-in observations and the final point estimateβ was computed using a sample size of 5 000 observations. As above, the MCMC sampler was initialized with the non-negative least-squares fit to the smeared data. However, since E F , we extended y to match the size of E by replicating the leftmost and the rightmost observations when computing the least squares fit. To form the bootstrap confidence intervals, R = 200 bootstrap replications were computed using resampling scheme 1. Running the MCEM iteration once to find the point estimateβ took 5 minutes. With the bootstrap, the running time of the whole algorithm was 6 h 13 min with the bootstrap computations parallelized on the four cores.
The convergence of the MCEM algorithm was confirmed using a plot similar to Figure 5 Carlo variation. During the MCEM iteration, the proposal acceptance rate remained at roughly 98 % and the average autocorrelation time increased from 3.0 to 8.6 indicating reasonable performance of the sampler throughout the whole iteration. As earlier, plots similar to Figure 6 were produced for each component of β for the final MCMC run in order to verify the appropriate convergence and mixing of the sampler.
In Figure 9, the bias-corrected unfolded intensityf BC of the Z boson invariant mass, along with 95 % pointwise basic bootstrap intervals, is compared with the Breit-Wigner shape of the true mass peak. We observe that empirical Bayes unfolding captures reasonably well the overall shape of the Breit-Wigner distribution with few undesired artifacts. The figure also shows a histogram estimate of the smeared intensity given by the observed event counts y divided by the 0.5 GeV bin size. We see that the unfolding algorithm is able to correctly reconstruct the location and width of the Z mass peak which are both distorted by the smearing in the ECAL. Moreover, thanks to the smoothness penalty and the Aristotelian boundary conditions, the intensity is estimated reasonably well in the 1 GeV regions in the tails of the intensity where no smeared observations were available.
However, in a closer examination, we observe that, starting from the top of the Z mass peak, the unfolded intensity is first slightly too wide on both sides of peak and then slightly too narrow. This artifact is likely to be a residual bias which is not accounted for by the bootstrap bias correction. The end result of this effect is better seen in Figure 10(a) which shows the unfolded intensity after subtraction of the true intensity f and normalization for an estimate of the expected total number of eventsλ tot = n i=1 y i . The figure shows both the original point estimatef (dashed curve) and the bias-corrected estimatef BC (solid curve) along with the 95 % pointwise basic bootstrap intervals. Although it cannot be directly deduced from this figure, it seems likely that, because of the remaining bias, the confidence intervals do not attain their nominal 95 % frequentist coverage across the whole spectrum. See Section 7 for further discussion on this observation. Note also that the bias correction has improved the point estimate only at the top of the Z boson mass peak but not at the sides of the peak.
To conclude this section, we show in Figure 10(b) the 95 % naïve empirical Bayes confidence intervals for the Z boson invariant mass. These intervals are again wider than the bootstrap intervals and hence seem to enjoy better coverage. Nevertheless, the interpretation of these intervals remains unclear.

Concluding remarks
We have studied a novel approach to solving the high energy physics unfolding problem involving empirical Bayes selection of the regularization strength and frequentist uncertainty quantification using the bootstrap. We have shown that empirical Bayes provides a principled way of choosing the hyperparameter δ with excellent practical performance in a wide variety of cases. As such, it provides an appealing alternative to classical methods for choosing the regularization strength, such as cross-validation or the Morozov discrepancy principle. Given the good performance of the approach, we anticipate empirical Bayes methods to also be valuable in solving other inverse problems beyond the unfolding problem.
It is nevertheless possible to find true intensities where empirical Bayes unfolding will not yield a good reconstruction. This happens when the smoothness penalty, i.e., penalizing for large values of f 2 2 , is not the appropriate way of regularizing the problem. For instance, if the true intensity f contains sharp peaks or rapid oscillations, the solution would potentially be biased to an extent where the bootstrap bias correction would be unlikely to sufficiently alleviate the situation. Naturally, in such a case, a more suitable choice of the family of regularizing priors {p(β|δ)} δ>0 should fix the situation. This highlights the fact that all the inferences considered here are contingent on the chosen family of priors and should always be interpreted with this in mind.
The other main component of our approach is frequentist uncertainty quantification of the solution using bootstrap resampling. We have shown that the bootstrap confidence intervals can serve as good estimates of the uncertainty of the solution, especially when there is little to moderate bias. However, with the Z boson dataset studied in Section 6, it is likely that these intervals do not attain their nominal confidence level. There are several possible explanations for this. First, we did not take into account the uncertainty stemming from the estimation of the smearing matrix K. Taking this uncertainty into account should widen the confidence intervals and hence improve coverage. The study of effective approaches to incorporating this uncertainty into the bootstrap procedure part of ongoing work. Second, the main problem in the unfolded Z boson invariant mass shown in Figure 9 is the presence of a bias in the form of small wiggles around the true intensity. The bootstrap is unable to probe this bias since the ill-posedness of K "smears away" these oscillations when we compute the product Kβ. In some sense, the bootstrap is blind to these artifacts and is hence unable to account for them either in the confidence intervals or the bias correction. To alleviate this problem, one could consider more elaborate bootstrap schemes. Perhaps one could, for example, sample the bootstrapped observations from the Poisson(K β ) distribution, where K is a regularized version of K.
Quite surprisingly, we found that in all our experiments the naïve empirical Bayes confidence intervals were longer than the bootstrap intervals, even though the former do not take into account the uncertainty regarding the choice of the hyperparameter δ. In practice, this is likely to mean that when there is little bias, the empirical Bayes intervals will overcover, but with larger bias, they might provide better coverage than the shorter boot-strap intervals. This means that the empirical Bayes intervals could also potentially serve as useful measures of uncertainty, especially since they are significantly cheaper to compute than the bootstrap intervals. Nevertheless, the indisputable advantage of the bootstrap intervals is that they enjoy a clear-cut frequentist interpretation, while the meaning of the empirical Bayes intervals is at best unclear. Interestingly, the recent theoretical results by Petrone et al. (2012) could potentially be used to prove the asymptotic coverage of the empirical Bayes intervals. Alternatively, one could perhaps use confidence distributions (Xie and Singh, 2013) to form a link between the empirical Bayes posterior and frequentist coverage. Such results would significantly help to demystify the meaning of the empirical Bayes intervals, at least in the asymptotic sense.