A unified perspective on modified Poisson likelihoods for limited Monte Carlo data

TL;DR

This paper unifies and extends methods for incorporating Monte Carlo uncertainties into Poisson likelihoods, proposing new probabilistic models that improve parameter estimation in counting experiments.

Contribution

It introduces a unified framework for existing approaches, derives three new probabilistic likelihoods based on gamma-Poisson mixtures, and demonstrates their advantages through extensive testing.

Findings

Probabilistic methods outperform Frequentist ones in various scenarios.

Using the exact compound Poisson distribution can be suboptimal for parameter estimation.

New formulas effectively incorporate Monte Carlo fluctuations into likelihoods.

Abstract

Counting experiments often rely on Monte Carlo simulations for predictions of Poisson expectations. The accompanying uncertainty from the finite Monte Carlo sample size can be incorporated into parameter estimation by modifying the Poisson likelihood. We first review previous Frequentist methods of this type by Barlow et al, Bohm et al, and Chirkin, as well as recently proposed probabilistic methods by the author and Arg\"uelles et al. We show that all these approaches can be understood in a unified way: they all approximate the underlying probability distribution of the sum of weights in a given bin, the compound Poisson distribution (CPD). The Probabilistic methods marginalize the Poisson mean with a distribution that approximates the CPD, while the Frequentist counterparts optimize the same integrand treating the mean as a nuisance parameter. With this viewpoint we can motivate three new probabilistic likelihoods based on generalized gamma-Poisson mixture distributions which we derive in analytic form. Afterwards, we test old and new formulas in different parameter estimation settings consisting of a "background" and "signal" dataset. The probablistic counterpart to the Ansatz by Barlow et al. outperforms all other existing approaches in various scenarios. We further find a surprising outcome: usage of the exact CPD is actually bad for parameter estimation. A continuous approximation performs much better and in principle allows to perform bias-free inference at any level of simulated livetime if the first two moments of the CPD of each dataset are known exactly. Finally, we also discuss the situation where new Monte Carlo simulation is produced for a given parameter choice which leads to fluctuations in the computed likelihood value. Two of the new formulas allow to include this Poisson uncertainty directly into the likelihood which substantially decreases these fluctuations.

Figures

Tables

Table 1. Table 1: Summary table of different likelihood approaches that modify the Poisson likelihood. The number of Monte Carlo events N 𝑁 N , the number of datasets N s ​ r ​ c subscript 𝑁 𝑠 𝑟 𝑐 N_{src} and the number of observed data k 𝑘 k are to be evaluated per bin. The Frequentist approaches (Barlow/Beeston and Chirkin) typically require multiple iterations (”niter”) of a numerical solver when there are more than one source dataset (Barlow/Beeston) or events with different weights (Chirkin).

Name	different weights	extra Poisson uncertainty	Per-dataset mean adjustment	$N\to \infty$ $\Rightarrow L\to L_{\mathrm{Poisson}}$ (”probabilistic”)	Computational complexity
Bohm/Zech Bohm2012	✗	✗	✗	✗	$𝒪\left(N\right)$
Chirkin Chirkin2013b	✓	✗	✓(not manifest)	✗	$𝒪\left(N\cdot \mathrm{niter}\right)$ ($N>1$)
Barlow/Beeston Barlow1993	✗	✗	✓(not manifest)	✗	$𝒪\left(N+N_{src}\cdot \mathrm{niter}\right)$ ($N_{\mathrm{src}}>1$)
Argüelles et al. Argueelles2019	✓	✗	✗	✓	$𝒪\left(N\right)$
Glüsenkamp Gluesenkamp2018	✓	✗	✓(not manifest)	✓	$𝒪\left(N\cdot k^2\right)$
$L_{\mathrm{Gen},\left(1\right)}$	✓	✓	✓(not manifest)	✓	$𝒪\left(N\cdot k^2\right)$
$L_{\mathrm{Gen},\left(2\right)}$	✓	✗	✓	✓	$𝒪\left(N_{\mathrm{src}}\cdot k^2+N\right)$
$L_{\mathrm{Gen},\left(2\right),\mathrm{eff}}$	✓	✗	✓(not manifest)	✓	$𝒪\left(N\right)$
$L_{\mathrm{Gen},\left(3\right)}$	✓	✓	✓	✓	$𝒪\left(N_{src}\cdot k^2+N\right)$

Equations

$$L(\theta )=\prod _{\mathrm{bins}i}p(k;\theta )=\prod _{\mathrm{bins}i}\frac{e^{-{\lambda }_i}\cdot {\lambda }_i^k}{k!}\mathrm{with}{\lambda }_i=\sum _{j=1}^N{w}_j(\theta )$$

$$\mathrm{E}[Z]$$

$$\mathrm{Var}[Z]$$

$$\stackrel{[}{\mathrm{Var}[N]=\mathrm{E}[N]}]={\mu }_N\cdot (\mathrm{Var}[W]+(\mathrm{E}[W]{)}^2)$$

$${\hat{\mu }}_Z$$

$${\widehat{\mathrm{Var}}}_Z$$

$$\stackrel{[}{{\widehat{\mathrm{Var}}}_N={\hat{\mu }}_N}]={\hat{\mu }}_N\cdot \left(\frac{1}{N}\sum _j(w_j-{\hat{\mu }}_W{)}^2+{\hat{\mu }}_W^2\right)$$

$$={\hat{\mu }}_N\cdot \left(\frac{1}{N}\sum _j{w}_j^2\right)=\sum _j{w}_j^2$$

$$L_{\mathrm{bin}}$$

$$\frac{k}{(1-t)}=\sum _j{\lambda }_j\cdot p_j{\hat{w}}_j=\sum _j\frac{k_{mc,j}\cdot p_j{\hat{w}}_j}{1+p_j{\hat{w}}_j\cdot t}$$

$$\ln (L_{\mathrm{bin}})$$

$$=-\sum _j\frac{(1+p_j{\hat{w}}_j)\cdot k_{mc,j}}{1+p_j{\hat{w}}_j\cdot t}+k\cdot \ln (\frac{k}{1-t})+\sum _j{k}_{mc,j}\cdot \ln (\frac{k_{mc,j}}{1+p_j{\hat{w}}_j\cdot t})-\sum _j\ln (k_{mc,j}!)-\ln (k!)$$

$$=-k-k_{mc,\mathrm{tot}}+k\cdot \ln (\frac{k}{1-t})+\sum _j{k}_{mc,j}\cdot \ln (\frac{k_{mc,j}}{1+p_j{\hat{w}}_j\cdot t})-\sum _j\ln (k_{mc,j}!)-\ln (k!)$$

$$\sum _j\frac{(1+p_j{\hat{w}}_j)\cdot k_{mc,j}}{1+p_j{\hat{w}}_j\cdot t}$$

$$=k_{mc,\mathrm{tot}}+\sum _j\frac{p_j{\hat{w}}_j\cdot (1-t)\cdot k_{mc,j}}{1+p_j{\hat{w}}_j\cdot t}=k_{mc,\mathrm{tot}}+k$$

$$\frac{d}{dt}\ln (L_{\mathrm{bin}})=0=\frac{k}{(1-t)}-\sum _j\frac{k_{mc,j}\cdot p_j{\hat{w}}_j}{1+p_j{\hat{w}}_j\cdot t}$$

$$t=\frac{k_{mc,1}\cdot p_1\widehat{w_1}-k}{p_1\widehat{w_1}\cdot (k+k_{mc,1})}$$

$$L_{\mathrm{bin},eq.}=\frac{e^{-(k+k_{mc,1})}}{k!k_{mc,1}!}\cdot {\left(\frac{k+k_{mc,1}}{1+1/(p_1{\hat{w}}_1)}\right)}^k\cdot {\left(\frac{k+k_{mc,1}}{1+p_1{\hat{w}}_1}\right)}^{k_{mc,1}}$$

$$\frac{k}{(1-t)}=\sum _{i=1}^{k_{mc,tot}}\frac{w_i}{1+w_i\cdot t}$$

$$L_{\mathrm{bin}}=\frac{e^{-{\sum }_i{w}_i{\lambda }_i}\cdot ({\sum }_i{w}_i{\lambda }_i{)}^k}{k!}\cdot \prod _{i=0}^{k_{mc,tot}}\frac{e^{-{\lambda }_i}\cdot ({\lambda }_i{)}^1}{1!}$$

$$L_{\mathrm{bin}}=\frac{e^{-\lambda }\cdot {\lambda }^k}{k!}\cdot \frac{e^{-\lambda \cdot \beta }\cdot (\lambda \cdot \beta {)}^{\alpha }}{\Gamma (\alpha +1)}$$

$$L_{bin}(\theta )$$

$$={\int }_0^{\infty }\frac{e^{-\lambda }{\lambda }^k}{k!}\cdot \left[\delta (\lambda -w_1)\ast \dots \ast \delta (\lambda -w_N)\right](\lambda )d\lambda$$

$$\to L_{bin}(\theta )$$

$$L_{bin}(\theta )=D_k\cdot \prod _{i=1}^N{\left(\frac{1}{1+w_i}\right)}^{1+\alpha /N}$$

$$D_k=\frac{1}{k}\sum _{j=1}^k\left[\left(\sum _{i=1}^N(1+\alpha /N)\cdot {\left(\frac{1}{1+1/w_i}\right)}^j\right)D_{k-j}\right]$$

$$L_{\mathrm{bin},eq.}$$

$$=\frac{(1/w{)}^{N+\alpha }\cdot \Gamma (k+N+\alpha )}{\Gamma (N+\alpha )\cdot k!\cdot (1+1/w{)}^{k+N+\alpha }}$$

$$L_{bin}(\theta )$$

$$L_{bin,exact}=\int \frac{e^{-\lambda }\cdot {\lambda }^k}{k!}\cdot p_{CPD}(\lambda )d\lambda$$

$$p_{CPD}(\lambda )=\sum _{k_{mc}=0}^{\infty }\mathbf{P}(k;\lambda )\cdot p_{\mathrm{fix}}(\lambda )=\sum _{k_{mc}=0}^{\infty }\frac{e^{-\mu }\cdot {\mu }^{k_{mc}}}{k_{mc}!}\cdot \delta (\lambda -k_{mc}\cdot w)$$

$$L_{bin,exact,equal}$$

$$=\sum _{k_{mc}=0}^{\infty }\frac{e^{-k_{mc}w}\cdot (k_{mc}w{)}^k}{k!}\cdot \frac{e^{-\mu }\cdot {\mu }^{k_{mc}}}{k_{mc}!}$$

$$L_{\mathrm{Gen},(1)}(\theta )$$

$$={\int }_0^{\infty }\frac{e^{-\lambda }{\lambda }^k}{k!}\cdot \left[\left(\sum _{i=0}^{\infty }{\mathrm{G}}_1({\lambda }_1;i+\alpha /N,1/w_1)\cdot {\mathrm{PG}}_1(i)\right)\ast \dots \ast \left(\sum _{j=0}^{\infty }{\mathrm{G}}_N({\lambda }_N;j+\alpha /N,1/w_N)\cdot {\mathrm{PG}}_N(j)\right)\right](\lambda )d\lambda$$

$$={\int }_0^{\infty }\frac{e^{-\lambda }{\lambda }^k}{k!}\cdot \left[{\mathrm{GPG}}_1({\lambda }_1)\ast \dots \ast {\mathrm{GPG}}_N({\lambda }_N)\right](\lambda )d\lambda$$

$$\sum _{i=0}^{\infty }$$

Full text

1. 1 Introduction

2. 2 Previous likelihood approaches for limited Monte Carlo

3. 2.1 Statistics of weighted Monte Carlo

4. 2.2 Frequentist Approaches

5. 2.2.1 Barlow/Beeston (1993)

6. 2.2.2 Chirkin (2013)

7. 2.2.3 Bohm/Zech (2012)

8. 2.3 Probabilistic approaches

9. 2.3.1 Glüsenkamp (2018)

10. 2.3.2 Argüelles et al. (2019)

11. 2.4 Summary

12. 3 Exact CPD for equal weights

13. 4 Further generalizations

14. 4.1 Generalization 1 - discrete marginalization of $k_{mc}$ with a Poisson-gamma mixture

15. 4.2 Generalization (2) - Modeling the CPD for each dataset

16. 4.3 Generalization (3) - Marginalization of the unknown Poisson mean

17. 4.4 Overview

18. 5 Incorporating extra prior information

19. 5.1 Reducing bias on $\mu_{n}$

20. 5.2 Reducing bias for $\mathrm{var}[W]$

21. 5.3 Handling empty bins

22. 5.3.1 Strategy 1

23. 5.3.2 Strategy 2

24. 6 Comparisons using Toy MC

25. 6.1 Situation 1: Likelihood ratio bias and coverage

26. 6.1.1 Equal weight test

27. 6.1.2 Background dataset with limited statistics

28. 6.1.3 Other scenarios

29. 6.2 A new MC set every minimization step

30. 7 Conclusion

31. A Derivation generalization (1)

32. A.1 Equal weights

33. A.2 General weights

34. B Derivation generalization (3)

35. B.1 Single Source Dataset

36. C Overview of some generic formulas

37. C.1 Poisson-gamma mixture form

38. C.1.1 Simple generalization

39. C.1.2 Extended generalization

40. C.2 Negative Binomial form

41. C.2.1 Simple generalization

42. C.2.2 Extended generalization

43. D Supplementary figures

1 Introduction

Monte Carlo (MC) simulations are used throughout science for statistical inference. In high-energy physics, modern experiments simulate physical processes with Monte Carlo generators to approximate intractable functions that map certain physical inputs $x$ to observable outputs $y$. In experiments related to astro-particle physics, for example, cosmic-ray simulations with MC-event generators like CORSIKA Heck1998 generate MC events following a certain cosmic-ray primary flux model with a specific composition and energy spectrum. This flux model defines the generation function $f_{\mathrm{gen}}(x;\theta_{0})$ which depends on spectral and composition parameters $\theta_{0}$ and is used to draw MC events from. The input $x$ here could for example be the true primary energy $E_{\mathrm{prim}}$ of a given cosmic-ray nucleus. After a full simulation involving physical interactions in the atmosphere and the experimental detector response, the end products of the shower development are usually recorded in bins of an observable quantity by counting. In the high-energy neutrino detector IceCube Halzen2010 , such an observable quantity is for example the total charge deposited in all Photo-multipliers from cosmic-ray induced muons Aartsen2016a or neutrinos Aartsen2016 . Initially, all recorded MC events have equal weights, which is illustrated in fig. 1 (a) in the central panel.

In order to save computing power, a MC simulation is often re-weighted with a function $f_{x;\mathrm{new}}(\theta)$, which might yield individual event weights different from unity (fig. 1 b) given by $w_{i}=\frac{f_{\mathrm{new}}(x;\theta)}{f_{\mathrm{gen}}(x;\theta_{0})}$. One can then bin the MC samples of the observable quantity and obtain an approximation for the expectation value in each bin $i$ as $\lambda_{i}=\sum_{j}w_{j}(\theta)$ (1 (a,b) in the lower panel). The MC events therefore serve as a mapping from parameters $\theta$ to expectation values $\lambda_{i}$, which can be used to write down a Poisson likelihood function for i.i.d. (independent and identically distributed) data as

[TABLE]

. As can be seen in fig. 1 (a,b) lower panel, a suitable change in $\theta$ results in a better match to the data which can be achieved via maximum-likelihood optimization Fisher1990 of $L$ (eq. 1) over $\theta$. It can also happen that re-weighting is not possible, and one has to simulate a whole new set of simulations for a given set of parameters $\theta$. This happens for example in the IceCube experiment for a specific algorithm to reconstruct neutrino directions and energies Chirkin2013 where the MC events now represent simulated Cerenkov photons emitted by a process parametrized by $\theta$. In addition to the extra computational cost in such a scheme, each MC sample is different from a previous one, even for similar generation parameters $\theta$ (see fig. 1 c). This generates extra complications for minimization routines due to fluctuations in the resulting likelihood function.

The issues we discuss in this paper arise because the Monte Carlo event count is finite which induces a statistical uncertainty in the evaluation of the related likelihood functions which depend on the MC events and their distribution in the observable bins (see eq. 1). We focus on Poisson-type likelihoods of the form in (eq. 1). The investigation is a continuation of a previous paper Gluesenkamp2018 which also discussed the more involved multinomial likelihood for unbinned approximations which we omit here.

In the first section, we review previous Frequentist and probabilistic approaches for this problem from the last 25 years and indicate how they are related. In the second section, we introduce new generalized Poisson-gamma mixture distributions that model the underlying statistical more manifestly or take care of extra uncertainty so far neglected. In the third section, we discuss techniques to incorporate extra prior information into the likelihood and a scheme to incorporate empty bins. In the final section, we compare different methods in various benchmark scenarios. All results of this paper are motivated by the use case in high-energy physics, but they are generally applicable in situations where MC simulations are used to calculate Poisson likelihoods.

2 Previous likelihood approaches for limited Monte Carlo

In the following we will revisit the common approaches in the literature that deal with limited Monte Carlo statistics. We will change the notation of some of them from their original publications in order to get a unified understanding.

2.1 Statistics of weighted Monte Carlo

In order to better understand the approaches in the literature, let us first take a look at the actual random variable describing weighted events in a bin $i$. For a given set of generation parameters $\theta_{0}$ and re-weighting parameters $\theta$, the weight of a Monte Carlo event ending up in bin $i$ follows a continuous and in general intractable distribution with random variable $W_{i}(\theta_{0},\theta)$, where $\theta_{0}$ includes a parameter $N_{tot}$ that specifies the total number of events simulated 111$N_{tot}$ could also be itself Poisson distributed, which would result in every $N_{i}$ be Poisson distributed even without the limit of the binomial distribution.. For a given $N_{tot}$, there is a resulting Poisson distributed number of MC events $N_{i}$ for bin $i$ with unknown underlying mean $\mu_{i}=N_{tot}\cdot p$. The probability $p$ is some unknown probability for an event to land in the specific bin $i$. As long as the number of bins is sufficiently large, say $n_{\mathrm{bins}}>10$, and $p$ sufficiently small, this follows from the usual Poisson limit of the binomial distribution. The actual observed weights are samples from the PDF corresponding to $W_{i}$. While the underlying PDF of $W_{i}$ is unknown, one could imagine repeating the MC experiment with the same parameters $\theta_{0}$ and $\theta$ a large number of times in order to visualize the PDF corresponding to $W_{i}$. Figure 2 shows such a PDF of $W_{i}$ for two different bins. The PDF is shown with and without detector-related smearing in the observable, which could for example be a reconstructed energy. This illustrates that the intractable weight distribution not only depends on the parameters $\theta$ and the corresponding bin, but also on detector effects.

From eq. 1 we see that we are really interested in a random variable describing the sum of weights $Z_{i}$ with $Z_{i}=\sum_{j=1}^{N_{i}}W_{i,j}$. From the preceding discussion we also see that $N_{i}$ is not fixed but itself Poisson distributed with unknown mean $\mu_{i}$. The resulting random variable $Z_{i}$ is called a compound Poisson distribution (CPD). To our knowledge this explicit statistical description of the problem was first pointed out in Bohm2014 . The problematic issue is that not only the PDF of $W_{i}$ is unknown (see fig. 2), but also the value of the underlying mean $\mu_{i}$. To get a better understanding of the CPD, figure 3 shows two example CPDs and how they are formed from the underlying weight distribution.

All approaches described later in this section can be interpreted as different ways to approximate the CPD given the observed MC samples.

The mean and variance of a compounding Poisson random variable $Z$ is given by

[TABLE]

where $Z=\sum_{j}W_{j}$, $W_{j}$ are independent and identical realizations of a continuous random variable $W$ and $N$ is a Poisson random variable defined on $\mathbb{N}^{+}_{0}$. Relationships 2 and 4 follow from the law of total expectation and total variance Weiss2006 , and in the last step we can exploit that mean and variance for a Poisson distribution are the same $\mu_{N}$. We can estimate the mean and variance using the sample mean $\widehat{\mu}_{Z}$ and sample variance $\widehat{\mathrm{Var}}_{Z}$ as

[TABLE]

where we just replace every term in eqs. 2 and 4 by the respective sample estimate, using $w_{j}$ to denote the observed weights. We used the biased sample variance of weights, but could have used the unbiased variance estimator instead. The resulting estimates are incidentally also the sample mean and sample variance of the sum of $N$ weighted Poisson random variables with $\lambda=1$ each, and they are widely used for plotting purposes to indicate the uncertainty of weighted MC histograms. The mean $\widehat{\mu}_{N}$ is unknown and by approximating it with $N$ there is always some sample uncertainty neglected. In the following we will sometimes use $N$ to denote number of Monte Carlo events or sometimes explicitly call it $k_{mc}$.

2.2 Frequentist Approaches

The approaches in this section are called ”Frequentist” because they can be thought of as ”profile likelihoods” Murphy2000 with additional nuisance parameters.

2.2.1 Barlow/Beeston (1993)

The systematic treatment of limited MC statistics for Poisson likelihood applications goes back at least 25 years to Barlow and Beeston Barlow1993 . The authors augment the Poisson likelihood with additional Poisson factors for MC datasets, effectively treating MC as additional data. Instead of a per-bin likelihood as described in eq. 1, the modified likelihood per bin looks like

[TABLE]

Each factor $j$ enumerates $N_{src}$ individual MC sources222We use slightly different notation from the original publication Barlow1993 . We write $k_{mc,j}$ instead of $a_{ji}$, $\lambda_{j}$ instead of $A_{ji}$, and directly use weights as described in section 6 in Barlow1993 ., for example a signal and a background dataset. The weights of Monte Carlo events are averaged per dataset as $\widehat{w}_{j}$. The parameters $\lambda_{j}$ are nuisance parameters and not of physical interest, and they should be optimized to maximize the values of the individual Poisson factors. The $p_{j}$ are global all-bin strength factors and can also be optimized as parameters, or fixed if relative strengths are known - in fact they could be incorporated into the weights directly. They are shared for all bins, and an overall scaling by $p_{j}\rightarrow c\cdot p_{j}$ of source $j$ would just increase the weights of source $j$ in all bins by the scaling factor $c$ in this case, making $p_{j}$ redundant if the weights are trusted in an absolute manner. This is actually the case in many modern applications. The nuisance optimization itself effectively handles some of the uncertainty from the finite MC event count. The authors Barlow1993 further show how the $j$ parameters $\lambda_{j}$ can be reduced to a single one per bin by taking the logarithm and finding the extremal values of $\mathrm{ln}(L_{\mathrm{bin}})$ in $\lambda_{j}$, which results in coupled equations for the $\lambda_{j}$ where $j-1$ of them can be eliminated. The resulting important relations eq. (25) and eq. (26) from Barlow1993 look like

[TABLE]

in our notation. Taking the logarithm of eq. 9, we can write

[TABLE]

where we use eq. 10 and insert $\lambda_{j}=\frac{k_{mc,j}}{1+p_{j}\widehat{w}_{j}\cdot t}$ to reduce the $j$ parameters $\lambda_{j}$ to a single parameter $t$ except for the first logarithm where we replace $p_{j}\widehat{w}_{j}\cdot\lambda_{j}$ with $\frac{k}{1-t}$. We further simplified eq. 12 to eq. 13 using the reformulation

[TABLE]

again making use of eq. 10. Assuming the $p_{j}$ are fixed, for example when they can be absorbed into the average weights $\widehat{w}_{j}$, we only have to find the value of $t$ that maximizes expression eq. 13 for every bin. We can form the derivative of eq. 13 with respect to $t$ and obtain

[TABLE]

which is just reproducing eq. 10. When there is only one source dataset this equation can be solved for $t$ exactly and one obtains

[TABLE]

which using eq. 10 can be used to substitute $\lambda_{1}$ in eq. 9 which yields the solution

[TABLE]

with no free parameters. This is also the exact solution when all weights are equal and the average weight $\widehat{w}_{1}$ is actually not averaged. In practice, one would insert $t$ into the log-likelihood expression eq. 13, but it is instructive to compare eq. 18 with other methods later. For more than one dataset we have to find the optimal value of $t$ numerically. The authors also discuss the case when no MC events are present in a bin. In that case they argue that $t$ should only depend on the largest $p_{j}$, which then can be used again to calculate the $\lambda_{j}$. While the method by Barlow/Beeston can not be directly interpreted as an approximation of the overall CPD, it can be interpreted as approximating individual CPD’s for each dataset, whose random variables are then summed. This will be exploited later in section 4 for a probabilistic generalization.

2.2.2 Chirkin (2013)

In 2013 Chirkin Chirkin2013b implicitly derived a certain generalization of eq. 9. The author starts with a multinomial likelihood and derives a formula that looks very similar to eq. 10 in the setting of what he calls ”weighted simulation without model errors”. In fact, the resulting equation is the same if we absorb the relative fractions $p_{j}$ into the weights and if we then treat each MC event independently as coming from its own ”source”, i.e. $k_{mc,j}=1$. From eq. 10 we then obtain

[TABLE]

which is similar to the equation obtained in section (5) in Chirkin2013b 333Using the notation of Chirkin2013b we have $d_{k}=k$, $t_{ki}=1/w_{i}$, $t_{d}=1$, $\xi_{k}=t$, $s_{ki}=1$.. The resulting corresponding likelihood formulation equivalent to eq. 9 is

[TABLE]

which now has a product over all individual MC events instead of whole datasets with averaged weights. The method by Chirkin can therefore be interpreted as a certain generalization of the Barlow/Beeston method. It omits weight averaging and gives better results for bias reduction in likelihood scans for a single weighted dataset as demonstrated in Gluesenkamp2018 , and sometimes also for multiple datasets (see later comparison). If all weights are equal the two methods are similar and their likelihoods are given by eq. 18 up to constant factors. The method can also be extended to include systematic uncertainties in a log-normal term that is dubbed ”model error” or handle bins with no MC using a constant ”noise” term. These terms have to be added by hand and are somewhat arbitrary.

2.2.3 Bohm/Zech (2012)

Bohm and Zech discuss another possibility to set up a pseudo-likelihood via

[TABLE]

where $\lambda$ is a nuisance parameter, $\alpha=\frac{(\sum_{i}w_{i})^{2}}{\sum_{i}w_{i}^{2}}$ and $\beta=\frac{\sum_{i}w_{i}}{\sum_{i}w_{i}^{2}}$. A single Poisson factor is used to encode the mean and variance of the sum of weighted MC in terms of effective counts to approximate the CPD. Because of the gamma function instead of the factorial this is not a real Poisson distribution and therefore not a real likelihood, so in general it can not be used for likelihood applications. The authors also only discuss it in the context of goodness of fit tests where they argue it is sensible. For equal weights, however, the gamma function turns into a factorial and the expression is similar to the previous two approaches and eq. 21 turns into eq. 18. It will still be interesting to compare the general pseudo-likelihood (eq. 21) to its probabilistic counterpart discussed later.

2.3 Probabilistic approaches

In contrast to the previous Frequentist approaches this section deals with solutions that integrate nuisance parameters instead of optimizing them.

2.3.1 Glüsenkamp (2018)

In a previous paper Gluesenkamp2018 we discussed a generalization of the Poisson likelihood of the form

[TABLE]

which can be thought of as directly expanding the delta factors into gamma distributions. Each gamma distribution can be derived via Bayesian inference and approximates the underlying weight distribution $W_{i}$ of the CPD. This means a convolution of these gamma factors is some approximation of the CPD (see fig 4) , the difference being a fixed number of terms and different individual PDFs. There is some prior freedom in this approach, which we parametrized by a parameter $\alpha$ which we here call $\alpha^{`}$. The parameter is shared among individual gamma distributions such that for equal weights it reduces to a single gamma factor with shape parameter $N+\alpha^{`}$. In Gluesenkamp2018 we observed optimal performance in terms of likelihood ratio bias with respect to infinite statistics for $\alpha^{`}=0$ for a single weighted dataset. We will later see that this is not true in general. The setting $\alpha=0$ yields the expected first and second of the CPD as $\mu_{G_{1}\ast\ldots\ast G_{N}}=\sum_{i}w_{i}$ and $\mathrm{var}_{G_{1}\ast\ldots\ast G_{N}}=\sum_{i}w_{i}^{2}$. It also showed practically identical likelihood ratio behavior as the method by Chirkin Chirkin2013b for a single source datasetGluesenkamp2018 . This is not surprising since a comparison of eq. 24 with eq. 20 shows that it can be interpreted as the probabilistic counterpart to the approach by Chirkin Chirkin2013b which enforces the sum constraint via optimization instead of probabilistic convolution (see fig. 4). The solution of eq. 24 was derived in Gluesenkamp2018 as an iterative sum via

[TABLE]

where

[TABLE]

with $D_{0}=1$. The formula allows for an efficient computation 444See also appendix C for the general form and different forms of writing. For equal weights, the result can be written as a simple Poisson gamma mixture as

[TABLE]

, where again we use $N=k_{mc}$. Comparing eq. 27 with eq. 18 we see that for equal weights the approach has similar dependence on the weights $w$ as the Frequentist approaches using the unique prior $\alpha=0$. Likelihood scans with equal weights therefore give numerically identical confidence intervals as the Frequentist methods.

2.3.2 Argüelles et al. (2019)

In Argüelles et al. Argueelles2019 the authors describe an Ansatz to encode the mean and variance of the CPD in a single gamma distribution which effectively approximates the CPD. It looks like

[TABLE]

where they add some additional effective parameters $a$ and $b$ to expand its possibilities. It can be thought of as a direct probabilistic counterpart of the pseudo-likelihood (eq. 21) described by Bohm and Zech Bohm2012 (see figure 4). In contrast to the pseudo likelihood, this Ansatz involves a marginalization with a proper PDF and therefore works for arbitrary weights. The authors argue that the choice of $a=1$, $b=0$ works best in their test in terms of proper coverage, in particular better than the choice $a=0$, $b=0$ which would just encode the sample mean and sample variance of the CPD in the gamma distribution. We will see later that this is not the case in general, in particular that it depends on the number of datasets combined and the MC statistics per dataset. In section 4.2 we describe an alternative parametrization that directly parametrizes the unknown Poisson mean and unknown weight moments which we later use to incorporate extra Prior information (section 5).

2.4 Summary

The relationship between the different approaches is highlighted in fig. 4. It rests on the interpretation that all approaches fundamentally approximate the CPD, an interpretation that has previously only been explictly made for the formulas put forward by Bohm/Zech Bohm2012 and Argüelles et al. Argueelles2019 and vaguely for our previous construction in Gluesenkamp2018 . For the Frequentist approaches by Barlow/Beeston Barlow1993 and Chirkin Chirkin2013b , the prevailing interpretation has been that the Monte Carlo data is simply treated as additional data, as indicated in fig. 4 by the upper formulas, respectively. However, by a simple change of variables from $\lambda_{i}$ to $\lambda_{i}^{*}$, one can see that one obtains a formula that almost resembles the probabilistic counterpart, in the case of Chirkin it resembles the probabilistic formula we proposed earlier in Gluesenkamp2018 . In this interpretation, the Frequentist approaches therefore also approximate the CPD via optimization in $\lambda_{i}$, which now resemble nuisance parameters where each parameter has its own gamma Prior. It rests on the property of the Poisson distribution that it can always be interpreted as a corresponding gamma distribution, whose shape parameter is shifted by one. It turns out that in order for an optimization to give similar results as the corresponding marginalization, the mode of the gamma distribution for optimization has to match the mean during marginalization. This is exactly fulfilled if the shape $\alpha$ is larger by one on the Frequentist side, since the mode of the gamma distribution is $\frac{\alpha-1}{\beta}$, while the mean is $\frac{\alpha}{\beta}$. The only difference remaining is the extra scaling factor $w_{i}$, that is being pulled out of the gamma factor. This, however, would exactly be the inverse of the Jacobian factor from the variable transform $\lambda_{i}^{*}=w_{i}\cdot\lambda_{i}$ and cancel, if one actually performed an integration.

3 Exact CPD for equal weights

What if we knew the PDF of the exact CPD, $p_{CPD}(\lambda)$, and integrate over it

[TABLE]

instead of using gamma factors in some form? In the special case of equal weights we can actually do this. In Bohm2014 it was argued that the we can approximate the CPD with a scaled Poisson distribution, but for equal weights the CPD is a scaled Poisson distribution. This results in $p(w)=\delta(\lambda-w)$ and $p_{\mathrm{fix}}(\sum_{i}w_{i})=p(w_{1})\ast\ldots\ast p(w_{N})=\delta(\lambda-k_{mc}\cdot w)$, which yields

[TABLE]

which is a superposition of the Poisson PDF to see $N$ weighted events with the probability distribution $p_{\mathrm{fix}}$ to see the weight sum of a given $N$ weighted events. So instead of approximating the CPD with gamma factors (see fig. 3) we can write down the exact CPD and integrate over it. The exact CPD for equal weights looks roughly equivalent to the example in figure 3 (b) where the weight distribution is very narrow. The final likelihood then can be written as

[TABLE]

which can be solved numerically. The parameter $\mu$ is an unknown mean, but can be determined by the average of many Monte Carlo runs. The result is a likelihood that does not depend on individual MC realizations anymore and involves the true PDF of the sum of weights without any approximations. It will serve as a crosscheck in the next section.

4 Further generalizations

We will now discuss some further constructions which generalize the previously discussed probabilistic approaches. An overview is shown in figure 5, which we will discuss in detail in the following.

4.1 Generalization 1 - discrete marginalization of $k_{mc}$ with a Poisson-gamma mixture

A first possible generalization starts with the interpretation of eq. 24, in which each gamma factor encodes a potential weight distribution. The number of weight distributions is fixed to the number of weights $N$. However, we can imagine to marginalize over the shape parameter of each gamma distribution with another discrete distribution, which models the Poisson uncertainty part of the CPD. If all weights are the same, this amounts to integrating out $N$. For practical reasons, we marginalize with a Poisson-gamma mixture distribution ($\mathrm{PG}$) instead of a Poisson distribution. We can write the result as

[TABLE]

where $\mathrm{PG}(i)=\mathrm{PG}(i;1+\beta/N;1/v)$, i.e. the shape parameter of the Poisson-gamma mixture is $1+\beta/N$ and the rate parameter is $\frac{1}{v}$. We also define the mixture of the gamma distribution with the Poisson-gamma distribution as $\mathrm{GPG}=\sum\limits_{i=0}^{\infty}\mathrm{G}_{N}(\lambda_{N};i+\alpha/N,1/w_{N})\cdot\mathrm{PG}_{N}(i)$. Here we have one mixture for each MC event. If two MC events have the same weight $w_{1}=w_{2}=w$, the convolution of two $\mathrm{GPG}$ distributions behaves similar as the convolution of two gamma distribution with the same rate parameter. We set $\beta=0$ and $v=1$ for simplicity and can write

[TABLE]

When all weights are equal, we can therefore combine all $\mathrm{GPG}$-distributions in eq. 35, and obtain

[TABLE]

which represents a direct generalization of the formula involving the exact CPD (eq. 31) for equal weights discussed in section 3. Instead of the delta factor we have a gamma factor, and instead of the Poisson distribution we have a Poisson-gamma mixture distribution. The expression eq. 35 therefore involves a certain generalization of a general compound distribution, where each of the continuous distributions $W_{i}$ is different.

The solution for $M$ different $\mathrm{GPG}$ distributions with different weights is derived in appendix A and yields

[TABLE]

with

[TABLE]

for the case $\alpha=\beta=0$ and $v=1$ which we use for all later comparisons. The general solution with arbitrary parameters is given in appendix C. Computationally, the calculation is only a little slower than eq. 25 and still scales as $\mathcal{O}(N\cdot k^{2})$. We tried the same calculation with a Poisson distribution as mixture instead of a Poisson-gamma distribution, but the result is computationally not as efficient. The Poisson-gamma distribution as the Posterior-Predictive distribution of observed counts is also well motivated, because it includes some uncertainty of the unknown underlying mean.

4.2 Generalization (2) - Modeling the CPD for each dataset

In the second generalization we model the CPD of each dataset independently, and afterwards join all CPDs via convolution, instead of approximating all datasets with a combined CPD (see fig. 5). Before we do so, we first introduce a slightly different parametrization of the CPD compared to the one used by Argueelles2019 . Encoding mean and variance es defined in eq. (2) and eq. (4) to be mean and variance of a general gamma distribution with parameters $\alpha$ and $\beta$ we obtain

[TABLE]

where we can observe that the unknown mean $\mu_{N}$ cancels out in $\beta$ but remains in $\alpha$. Choosing $\mu_{N}=k_{mc}+\frac{k_{mc}\cdot\sum_{j}w_{j}^{2}}{(\sum_{j}w_{j})^{2}}$ results in $\alpha=\frac{(\sum_{j}w_{j})^{2}}{\sum_{j}w_{j}^{2}}+1$ which is the choice that was made in Argueelles2019 . Instead of an additive nuisance parameter $a$ one can therefore choose a multiplicative parameter that is naturally given by $\mu_{N}$, which compared to the other choice has an intuitive interpretation as the unknown mean of the Poisson part of the CPD. We also do not have to use the sample based mean and variance for the weight distribution $W$, but can use other ways to encode the moments. For example a wider variance might be reasonable to model very few events (see section 5) or we can use asymptotic Monte Carlo simulations to get the exact moments for crosschecks (see section 6). The likelihood function with several such CPDs joined by convolution looks like

[TABLE]

where each gamma factor has its parameters $\alpha$ and $\beta$ encoded as described above. Compared to the single CPD Ansatz by Argüelles et alArgueelles2019 , the advantage is that large differences in statistics between datasets can be explicitly modeled in the respective $\mu_{N,j}$. It also has an interpretation as a direct probabilistic counterpart of the Barlow/Beeston Barlow1993 Frequentist Ansatz. In contrast to Barlow/Beeston, however, each CPD can take into account the full weight distribution, instead of averaging the weights. The analytic solution is known from Gluesenkamp2018 and given in appendix C. It scales computationally as $\mathcal{O}(N_{\mathrm{src}}\cdot k^{2}+N)$. We can also introduce an effective version

[TABLE]

which is just a single CPD for all weights, similar to Argüelles et al Argueelles2019 , but with the different parametrization of $\alpha$ and $\beta$ as introduced above. This allows to model per-dataset mean adjustments or weight distributions approximations as discussed later in section 5.

4.3 Generalization (3) - Marginalization of the unknown Poisson mean

Starting with the second generalization, we can go one step further and marginalize the unknown mean $\mu_{N,j}$ of the Poisson part of each CPD $j$ with another gamma distribution. The idea is to capture some of the uncertainty from ignorance about the unknown mean of the number of MC events per bin. We obtain

[TABLE]

because we can write the mixture of a Poisson distribution with a convolution of gamma distributions as a discrete convolution of Poisson-gamma mixtures (see appendix C. The scaling parameter $s$ has a natural choice of $s=1$. Choosing $s=2$ results in a gamma distribution with twice the variance but the same mean, and can be used to mimic the variance of a Poisson-gamma mixture distribution. The term $\alpha_{j}^{*}=\frac{\mathrm{E}[W_{j}]^{2}}{\mathrm{Var}[W_{j}]^{2}}$ is again decoupled from the mean $\mu_{N_{j}}$, similar to the previuos generalization. To simplify terminology we use $\mathrm{GG}$ to denote a gamma distribution whose shape is marginalized by another gamma distribution, and $\mathrm{PGG}$ to denote a Poisson-gamma mixture distribution whose shape is marginalized with another gamma distribution. Several ways to calculate the solution exist, but we found the convolutional form (eq. 49) to be the most efficient way. Its computational complexity approximately scales as $\mathcal{O}(N_{src}\cdot k^{2}+N)$. The calculation of an individual $\mathrm{PGG}$ distribution is shown in appendix B and yields

[TABLE]

with $\alpha_{j}^{*}=\frac{1}{k_{mc,j}}\frac{(\sum_{j}w_{j})^{2}}{\sum_{j}w_{j}^{2}}$ and $\beta_{j}=\frac{\sum_{j}w_{j}}{\sum_{j}w_{j}^{2}}$ where we inserted sample mean and sample variance of the weight distributions. If not otherwise stated, we use $s=1$. The unsigned Stirling numbers of the first kind can be precomputed once via the recursion relation Abramowitz1970

[TABLE]

and then only require a table lookup.

4.4 Overview

An overview of the properties of all approaches is given in table 1. The first column with weight uncertainty shows that most approaches approximate the underlying weight distribution to some extent. The exceptions are the approaches by Barlow/Beeston and by Bohm/Zech which approximate it essentially with a delta peak by averaging the weights. The second column shows which approaches include extra uncertainty for the unknown mean of the Poisson part of the CPD. Generalizations (1) and (3) are the only ones which do this, either by a discrete marginalization of the observed MC counts or by marginalizing the unknown mean with another probability distribution. The third column discusses the ability to adjust the expected mean for different datasets independently. The Frequentist approaches by Chirkin and by Barlow/Beeston in particular do this to some extent, as the profile likelihood fit implicitly varies parameters that can be interpreted as means of different Poisson components. However, in the second interpretation (see fig. 4), this is not the case anymore. The fourth column indicates that only the probabilistic approaches converge to the absolute value of the Poisson likelihood as more Monte Carlo simulation is used. The Frequentist approaches only capture the shape, but their absolute values are meaningless. The last column shows approximate computational complexities. While some approaches scale similarly, for example generalization (2) and (3), one has to remember that individually they require different mathematical operations, so in practice they might differ by a factor of a few.

5 Incorporating extra prior information

While the old and new formulas capture more uncertainty than the standard Poisson likelihood, the approximations are often biased, in particular if only a few MC events are present. This can be seen directly in section 4.2: the mean $\mu_{N}$ is approximated by the observed MC counts and the weight distribution mean $\mu_{W}$ and variance $\mathrm{var}[W]$ approximated by the sample mean and variance of observed weights. If only a few events are present in a given bin, these estimates can be biased. For $\mu_{N}$ this happens when the Monte Carlo livetime is so low that several bins are empty, which indicates that the true mean $\mu_{N}$ is smaller than unity per bin. For $\mu_{W}$ and $\mathrm{var}[W]$, this happens when the weight distribution is rather wide and only few events are present. Here we discuss a few approaches that can remedy these biases to some extent.

5.1 Reducing bias on $\mu_{n}$

If the Monte Carlo livetime is very low it can happen that some bins are empty. Using the observed number of Monte Carlo events in the remaining bins that do see events can lead to an over estimation of the underlying mean. One strategy to circumvent this is to calculate the average number of Monte Carlo events over all bins that are eventually filled by the simulation. If this average is smaller than $1$, we subtract the difference from the sample mean of observed events in every bin. Once the average numer of Monte Carlo events per bin is above $1$, no correction is performed, so it only affects very low statistics. Importantly, it has to be done for every Monte Carlo dataset separately. If one has two datasets with very different statistics, and calculates this adjustment based on the combined average number of events per bin, it does not work. This is the reason generalizations (2) and (3) make this per-bin mean adjustment natural, because every dataset $j$ has its own mean $\mu_{N},j$ that gets adjusted. However, it is in principle also possible to do this for most of the other approaches on a per event basis. This is what we call ”not manifest” in table 1.

5.2 Reducing bias for $\mathrm{var}[W]$

For very few events the sample variance might give a biased estimate of the true variance of the weight distribution. In particular, for a single event, the sample variance is zero, which is always an underestimation if the true weight distribution has a non-zero with. This happens whenever we have non-equal weights. Simply using the unbiased sample variance often does not help. One possibility to counteract this behavior is to assume that each weight represents a gamma distribution, and the total weight PDF is sum of individual gamma PDFs. The variance of a single gamma distribution with rate parameter $1/w$ is $w^{2}$, so the variance of a superposition of multiple gamma distributions is $\mathrm{var}[W]=\frac{1}{N}\sum_{i}w_{i}^{2}$. We will see in section 6 how this compares to the standard biased sample variance.

5.3 Handling empty bins

Empty bins must be handled with care. In particular, it is important to differentiate between datasets. For one dataset, a bin might be empty, while for another dataset it might be filled with many Monte Carlo events. If nothing is done about that, the resulting parameter estimation, which is likely connected only to one of the datasets, will be biased. In Barlow1993 , the strategy advocated was to add the dataset with the largest potential contribution based on all datasets that have no Monte Carlo events in the bin to add to the bin. However, this strategy adds discrete jumps to the likelihood, because for different parameters a different dataset might have higher contributions. Another strategy is discussed in Chirkin2013b , which is just to add a constant ”noise” contribution to all bins that are empty. In this strategy, however, no differentiation is made between datasets, which is bad for the reason stated above. The weight of the pseudo event is chosen to be the largest weight from all other bins, similar to the suggestion in Barlow1993 , but we do it for all datasets without any differentiation. The corresponding CPD approximations of these single ”pseudo” Monte Carlo events again typically suffer the problem of an over estimation of the mean, which is discussed in the section 5.1, and are almost by definition bias for the weight distribution. Therefore we typically combine the strategies for empty bins with the mean adjustment in section 5.1 and increased weight variance 5.2. Here we discuss two different strategies to incorporate empty bins. Both strategies basically fill up bins with ”pseudo counts” for each dataset, which avoids the problem of jumps in the likelihood. The difference between the strategies involves the type of bins which are filled up.

5.3.1 Strategy 1

In the first strategy we only fill up bins that have at least one Monte Carlo event from any dataset. Once a bin is identified, we add a single ”pseudo” Monte Carlo event for every other dataset that has no events in it.

5.3.2 Strategy 2

In the second strategy we fill up all bins that can potentially see Monte Carlo events, even those that see no events from any dataset yet. This is a more aggressive strategy, and leads to potentially larger bias for very low counts. The main aim of this strategy is to reduce likelihood fluctuations in the case when a new Monte Carlo simulation is run every parameter query (section 6.2). These fluctuations tend to be larger when totally empty bins are neglected, which happens by definition in strategy 1.

6 Comparisons using Toy MC

To compare the various approaches, we will use a toy Monte Carlo setting of a falling background energy spectrum with a peaked signal spectrum on top. An artificial detector is simulated with energy smearing and energy-dependent detection efficiency and events are observed in some observable space. We might be interested in determining parameters related to the ”signal” distribution like normalization and position, or parameters related to the ”background” distribution, like normalization or spectral index. Figure 6 shows

likelihood scans for the different standard approaches and the respective Monte Carlo and data distribution in observable space for the best-fit value of the Poisson likelihood. Assuming asymptotic behavior of the log-likelihood ratio $\lambda=-2\cdot\Delta\mathrm{LLH}(\hat{p},p_{0})\sim\chi_{\nu}^{2}$ PDG2018 , one can determine the significance of exclusion of a given parameter value $p_{0}$, given $\hat{p}$ the best fit value. In the following tests, instead of using the log-likelihood ratio in the usual way, we show the bias of the log-likelihood ratio $\lambda$ (in $\sigma$ equivalents) for a specific data realization compared to infinite Monte Carlo statistics for a given true parameter $p_{0}$. We also show the expected coverage vs. the actual coverage by evaluating 500 Toy experiments for a given simulated livetime.

6.1 Situation 1: Likelihood ratio bias and coverage

6.1.1 Equal weight test

The simplest test one can perform is a fit involving only the background dataset. We use a simple weighting function that just scales normalization (see fig. 1) and all weights for the background dataset are the same. In this particular case of equal weights, we can calculate the CPD exactly (see section 3). We perform a large number of Monte Carlo simulations at each statistical level to obtain a good estimator of the true mean $\mu_{N}$, while $\mu_{W}=w$ and $\mathrm{var}[W]=0$ are known since we deal with equal weights. We compare this ”exact” CPD likelihood with the new generalizations (1), (2), and (3), and a particular form of generalization (2) where we replace the unknown mean $\mu_{N}$ not by the observed number of MC events $k_{mc}$, but also by the true mean (”exact $\mathrm{E}[N]$”) obtain from a large number of MC realizations. This can be understood as a continuous approximation of the exact CPD. The comparison of these different likelihoods in terms of $\lambda$-bias and coverage is shown in fig. 7.

The first surprising observation is that the exact CPD ($L_{\mathrm{CPD,exact}}$) has a large bias at low average MC events per bin, although one would expect that the exact CPD works at any level of Monte Carlo statistics. This probably has to do with the multi-modal structure of the exact CPD. The continuous approximation ($L_{\mathrm{Gen,2,Z1,exactE[N]}}$) on the other hand is basically unbiased everywhere. We find that generalization (1), which is similar in nature to the exact CPD, performs not so good, which is related to the fact that the exact CPD itself performs badly, and we only plot it for the simple case without any modifications. Generalization (2) and generalization (3), which are again continuous approximations, perform better and are almost equivalent in this scenario. We can also observe how the mean adjustment is important to reduce bias, while the strategy for empty bins or using a more generous estimate for the weight variances do not change things too much in this particular scenario. Comparing the approach using the adjusted mean ($L_{\mathrm{Gen,2,M}}$) with the exact formula ($L_{\mathrm{Gen,2,Z1,exactE[N]}}$), we see that in the intermediate statistics regime there still is some bias in the likelihood ratio that is not removed by the mean adjustment strategy applied here. One could think about addressing this problem by some modification of the strategy, but we do not further test it here. Figure 8 shows the best performing cases of generalization (2) and generalization(3) in comparison to older likelihood approaches. In general, the new approaches show lower bias and better coverage over a wide range of MC statistics.

6.1.2 Background dataset with limited statistics

In this test setting we fix the background dataset at limited statistics, while the signal statistics are increased step-by-step. We fit the normalization of the signal peak. Such a situation with limited background statistics can happen if hard cuts are applied in the data selection process and computational resources are not sufficient to still have enough simulated background events surviving at final level. We compare various combinations of treating empty bins with incorporating prior information. Again we first compare various strategies with the new formulas as shown in fig. 9.

We can obtain the exact moments from high statistics simulations and use these moments at the respective level, which results in a unbiased likelihood ratio at all statistic levels for crosschecks (($L_{\mathrm{Gen,2,Z1,exactmoments}}$)). If we use the exact Poisson mean, but use sample estimates for the weight moments ($L_{\mathrm{Gen,2,Z1,exactE[N]}}$), a bias becomes visible. We can compare this result to the likelihood of generalization (2) with mean adjustment, which shows the mean adjustment is nearly as good as using the exact Poisson mean. We can also see that filling up the empty bins with pseudo counts is necessary here. All likelihood approaches that do not fill up empty bins (no ”ZB”) end up heavily biased even with infinite signal Monte Carlo, only because the background Monte Carlo stays limited. Additionally, usage of a more flexible weight variances estimator (”U”) can be helpful. For generalization (3), though, having a more unbiased estimator (”U”) seems to introduce a slight bias at high statistics, which we have not fully understood yet. Potentially the total variance of the estimators is too large here (compare also fig. 12 in appendix D without ”U”). Lastly, the effective version of generalization (2) does not seem to perform well in this scenario for very low statistics, but for higher statistics it approaches the standard formula of generalization (2). In figure 12, which we show in appendix D, we compare again generalization (2) against older approaches from the literature and also show the maximum likelihood estimator (MLE) bias and variance. We see again that all approaches from the literature fail when empty bins are neglected. Interestingly, the standard Poisson likelihood performs quite good in a low-statistic certain region, which has to do with very low MLE bias compared to all modifications.

6.1.3 Other scenarios

We also tested several other scenarios, whose results we show in appendix D. One scenario deals with more statistics and a smaller relative peak where we fit the peak normalization. The results are shown in fig 13. This situation is more akin to the scenario that was studied by the authors in Argueelles2019 (”$L_{\mathrm{ASY,eff}}$”). We see that our likelihoods with mean adjustment perform better, in particular the effective version of generalization (2) (”$L_{\mathrm{Gen2,eff}}$”) , which shows that a different parametrization is really crucial. In yet another scenario we increase the livetime of background and signal simultaneously and fit the spectral index (see fig. 14). Again, the results are similar.

6.2 A new MC set every minimization step

In certain applications we might simulate a new set of MC events for every successive parameter call. This is a usual technique in IceCube to determine the ice properties or to perform expensive particle reconstructions Chirkin2013 . Even when the parameters are the same, a new set of MC events will have different weights and a different simulated count per bin (see 1 fig. 1). In this situation, the uncertainty is directly visible in fluctuations of the corresponding likelihood function. In fig. 10 we compare various versions of generalization (1) and generalization (3), which directly model this uncertainty within the PDF, for a likelihood scan of the normalization of the background dataset.

In this scenario all weights are the same. We test different parameter values of the scaling parameter $S$ from generalization (3) (see section 4.3). A larger scaling parameter leads to a larger uncertainty being considered in the likelihood function, and correspondingly the fluctuations go down. However, we can observe that the likelihood bias goes up. A useful strategy might therefore be to start with a large parameter value of $S$ and decrease it during the optimization process. Similarly, we could use generalization (1) and modify the parameter $v$ (see section 4.1) and introduce another scaling factor for the respective shape parameters of the Poisson-gamma mixtures. This should have a similar effective behavior. Additionally, we see that filling up the empty bins (Z2) is a second ingredient which is important. In fig. 11 we compare the best performing generalization (3) with previous approaches. Here we observe that it is also possible to lower relative fluctuations by removing constant parts from the likelihood function which lowers the absolute values, and thereby also the relative fluctuations. The approaches described in Chirkin2013 Barlow1993 Gluesenkamp2018 yield exactly similar fluctuations in this case. However, they can not compete with a direct uncertainty modeling.

7 Conclusion

The limited statistics of Monte Carlo simulation leads to extra uncertainty that has to be taken into account for correct parameter estimation. In this paper, we discussed several solutions from the literature that modify the Poisson likelihood to tackle this problem. We showed that all the approaches of this type can be interpreted as approximating the underlying probability distribution of the sum of weights (the CPD) and then use this approximation to solve a marginalization integral over the Poisson mean. In probabilistic approaches, the integral is calculated exactly. In Frequentist approaches, the integrand is optimized over the mean, not marginalized. This is an alternative viewpoint to the prevailing interpretation of Frequentist approaches (Barlow1993 Chirkin2013b ), which is that these simply treat the Monte Carlo data similar to real data in the likelihood function.

In the new perspective, the central question is really how to approximate the mean $\mu_{N}$ of the Poisson part of the CPD, as well as the mean $\mu_{W}$ and variance $\mathrm{var}[W]$ of the continuous ”weight part” of the CPD. We introduced three new generalizations of existing probabilistic approaches that are well-motivated in this context. In generalization (1), we marginalize each gamma distribution from the probabilistic Ansatz in Gluesenkamp2018 with another Poisson-gamma mixture distribution. For equal weights, this distribution involves a marginalization of a Poisson-gamma mixture compound distribution, which can be viewed as a generalization of the exact CPD with more variance. In generalization (2), we model the CPD of each dataset independently with a gamma distribution, which are then convolved to form the total CPD. It can be viewed as a generalization of the approach by Argueelles2019 in combination with the convolutional formula introduced in Gluesenkamp2018 . It is also the direct probabilistic counterpart to the approach by Barlow et al Barlow1993 . Each gamma distribution is parametrized by the relevant CPD moments $\mu_{N}$, $\mu_{W}$, $\mathrm{var}[W]$ of a given dataset, which are usually taken to be the sample moments, but they don’t have to be. We also discuss an effective version of generalization (2) which encodes a single CPD for all datasets, similar to Argueelles2019 , but with a different parametrization that allows to incorporate extra Prior information in a more natural way. In generalization (3), we marginalize the unknown mean $\mu_{N}$ of each CPD from generalization (2) with yet another gamma distribution to incorporate extra sampling uncertainty directly into the PDF. All of the new probability distributions are analytically solvable.

In the last part of the paper we test these new formulas in a Toy MC parameter estimation setting with two datasets, a background dataset and signal dataset. We find that generalization (1) performs quite poor in these tests in terms of coverage and likelihood ratio bias. We also surprisingly find that the exact CPD performs worse than a continuous approximation with a gamma distribution, especially at very low statistics. We suspect that this has to do with the multi-modal form of the exact CPD. This explains why generalization (1) performs so poor for correct coverage, because it is just a slightly more general version of the exact CPD in the equal weights setting and has similar multimodal structure. Generalizations (2) and (3), which are based on continuous approximations of the overall CPD, perform better in all coverage tests.

If the moments are known exactly, a continuous approximation of the CPD as it is done in generalization (2) is in a certain sense optimal: it never undercovers at any level of Monte Carlo livetime. Achieving good coverage therefore reduces to finding good approximations to the moments $\mu_{N}$, $\mu_{W}$ and $\mathrm{var}[W]$ - the better these approximations, the better the coverage properties. Here, the usual sample moments $\mu_{N}=k_{mc}$, $\mu_{W}=\frac{1}{k_{mc}}\sum_{i}w_{i}$, $\mathrm{var}[W]=\frac{1}{k_{mc}}\sum_{i}w_{i}^{2}$ are often biased which leads to undercoverage.

We can remedy this problem to some extent by incorporating extra Prior information about the Monte Carlo into these moments. One of these is a ”mean adjustment” based on the average Monte Carlo count in all bins. The mean, and thereby also the mean adjustment, can be interpreted as a hyperparameter similar to the parameter $a$ in the effective likelihood discussed in Argueelles2019 , and we show how they are precisely related. In contrast to the choice of an effective likelihood with $a=1$ discussed in Argueelles2019 , the mean adjustment is conceptually easier to interpret and gives better results in all test scenarios.

We also discuss a certain strategy how to treat bins with zero Monte Carlo events for a given dataset. In contrast to previous suggestions from the literature (see Barlow1993 ) how to treat such scenarios, we find that a good choice is to fill up bins with pseudo counts such that every dataset has at least a single event in every bin, indifferent to the relative weights between datasets. This has the advantage that no jumps in the likelihood can appear.

In a final test we studied the scenario of new simulation Monte Carlo simulation every parameter step, as it happens for example in the application of direct Cerenkov photon simulation in the Icecube observatory Chirkin2013 . Here, we observed that generalization (1) and generalization (3), which take into account the Poisson uncertainty from repeated sampling steps, greatly reduce fluctuations in successive likelihood evaluations compared to all other approaches. Given that generalization (3) has better properties in terms of log-likelihood ratio bias, it seems to be better suited than generalization (1) in this use case.

Our overall recommendation for parameter estimation situations can be summarized as follows:

Normal use case of weighted simulation: Use generalization (2) with strategy 1 for empty bins and mean adjustment. If the computational cost is too high for the application, use the effective version of generalization (2). 2. 2.

Re-simulation every parameter query: Use generalization (3) with strategy 2 for empty bins and a scaling parameter $s>1$ to minimize fluctuations. Potentially decrease the scaling paramter during the minimization process to unity to decrease bias.

Besides the tests with exact moments using perfect knowledge from ”infinite Monte Carlo”, it was never possible to reduce the bias of the log-likelihood ratio completely at all levels of simulated livetime. Blindly trusting a formula can therefore lead to trust in a systematically wrong result if log-likelihood ratio values are directly translated into significances. The way out are ensemble-based test statistic constructions. These might again be problematic because Monte Carlo is expensive to produce. If ensemble-based constructions are possible, however, we recommend to always use the new formulas because they capture the underlying statistics of the CPD. In general, it is always a good idea to calculate the log-likelihood ratio with multiple likelihood functions, at least also the standard Poisson likelihood, and conclude from the differences between the different approaches if enough Monte Carlo statistics have been produced. Probabilistic approaches allow to additionally compare absolute likelihood values. If those absolute values disagree substantially, for example between the Poisson likelihood and the modified likelihood in use, it is a strong hint that one must be careful on the interpretation of the result. Implementations of the new formulas can be found on http://www.github.com/thoglu/mc_uncertainty.

Acknowledgements

We would like to thank Carlos Argüelles, Marty Cohen, Austin Schneider and Tianlu Yuan for stimulating discussions.

Appendix A Derivation generalization (1)

This section computes the generalization described in section 4.1.

A.1 Equal weights

We start with equal MC weights for simulated events in a given bin. In this case, the generalized per-bin likelihood function is given by suitable marginalization as Gluesenkamp2018

[TABLE]

where the result is a standard gamma-Poisson mixture $\mathrm{PG}(k;k_{mc}+\alpha,1/w)$ distribution. For more information on the parameter $\alpha$, which can be seen as a hyperparameter parametrizing a gamma prior, in particular its special value $\alpha=0$, see Gluesenkamp2018 . Let us now extend this construction by further marginalizing out the discrete counts $k_{mc}$ via

[TABLE]

and obtain a new likelihood function $L_{\mathrm{\textbf{P},finite,twice,eq.}}$ that has been marginalized twice, whose probability distribution we abbreviate as $\mathrm{PGPG}$, because we marginalize a Poisson-gamma mixture with another Poisson-gamma mixture. This second marginalization is discrete, and happens with another gamma-Poisson mixture $\mathrm{PG}(k_{mc};\langle k_{mc}\rangle+\beta,1/v)$ which is similar in structure to the Poisson mixture in eq. 71, but with different parameters 555We could have used a standard Poisson distribution for marginalization instead of a Poisson-gamma mixture. However, it turns out that such an Ansatz gives a more complicated result for general weights which are discussed in the next section. In particular, essential singularities appear which make the end result harder to compute.. We will later set the ”count weight” $v$ to $v=1$, but for now keep it in this general form. A weight of unity makes sense since we are now interested in marginalizing pure counts, which are unweighted by construction. Another nuisance parameter $\beta$ is introduced to handle Prior freedom, similar to the introduction of $\alpha$ in Gluesenkamp2018 . In eq. 72 the variable $k_{mc}$ denotes the observed number of MC events, while in the marginalized expression (eq. 73) the variable that is integrated out is $k_{mc}$ and the observed MC counts as $\langle k_{mc}\rangle$. Equation 73 can be solved using Egorychev rules (Egorychev1984 p. 19) via

[TABLE]

Using the Taylor expansion of $f(z_{1})$ and $g(z_{1})$ around $z_{1}=0$, we can combinatorially combine the relevant terms to contribute to the residue of the contour integral at $z_{1}=0$. The result is then readily written down as

[TABLE]

where the constants $K$ and $C$ are defined in eq. 77 as $K=\left(\frac{1}{1+1/w}\right)^{k}\left(\frac{1}{1+w}\right)^{\alpha}\cdot\left(\frac{1}{1+v}\right)^{\langle k_{mc}\rangle+\beta}$ and $C=\frac{v}{(1+v)(1+w)}$.

A.2 General weights

For general weights we start with eq. 35 and rewrite it as

[TABLE]

where we first pull out all the infinite sums, then rewrite the integral over the convolution of gamma factors with a Poisson factor in combinatorial form (see appendix C) and then again combine each summation with the respective terms. The result is a combinatorial sum over $\mathrm{PGPG}$ distributions, each of which we calculated the result in the previous section. We can continue and replace every $\mathrm{PGPG}$ distribution in eq. 85 by its respective contour integral (eq. 79), which results in

[TABLE]

where we use $\alpha_{j}^{*}=\frac{\langle k_{mc,j}\rangle\cdot\alpha}{N}$ and $\beta_{j}^{*}=\frac{\langle k_{mc,j}\rangle\cdot\beta}{N}$ to share the nuisance parameter among groups of events with different weights and the definition of $K_{j}$ and $C_{j}$ is the same as $K$ and $C$ for given a weight group in the previous section. Again we use $\langle k_{mc}\rangle$ to denote observed number of MC events, and $k_{mc}$ as a variable to sum over. We can now solve the constrained sum over $k_{j}$ using Egorychev’s rules for constrained summation indices Egorychev1984 . The constraint here is $\sum k_{j}=k$, which is encoded with complex weight variables $\tau$ of equal magnitude. Equation eq. 87 with these weight variables $\tau$ and an extra outer contour integral then can be re-written as

[TABLE]

using $K_{j}\equiv K_{j}^{{}^{\prime}}\cdot E_{j}^{k_{j}}$, $E_{j}=\frac{1}{1+1/w_{j}}$, $C_{j}=\frac{v_{j}}{(1+v_{j})(1+w_{j})}$ and $K_{j}^{{}^{\prime}}=(\frac{1}{1+w_{j}})^{\alpha^{*}}\cdot(\frac{1}{1+v_{j}})^{\langle k_{mc}\rangle+\beta_{j}^{*}}$ from the definitions in the previous section. The terms we call $f$ and $\varphi$ are highlighted to indicate the structure used for the main theorem described in Egorychev1984 (p. 19, 1.4.2). The term $\Delta_{k}$ represents an efficient solution Ma2014 of the contour integral in eq. 91, and is written as

[TABLE]

with $\Delta_{0}=1$. One can also write the sum going over all individual events $N$ instead over $M$ groups of events that share a similar weight, which results in $k_{mc,j}=1$ for each summand given by the last equal sign. The nuisance parameters $\alpha_{j}^{*}=\frac{\langle k_{mc,j}\rangle\cdot\alpha}{N}$ and $\beta_{j}^{*}=\frac{\langle k_{mc,j}\rangle\cdot\beta}{N}$ are typically set to zero since $\alpha=\beta=0$ is a natural prior choice (see also Gluesenkamp2018 ). The parameter $v_{j}$ is further set to unity, because it came about as the scaling parameter of the gamma distribution in the Poisson gamma mixture describing pure counts. The general case of convolutions of arbitrarily mixed gamma distributions and gamma-Poisson-gamma mixtures is described in appendix C.

Appendix B Derivation generalization (3)

B.1 Single Source Dataset

For a single source dataset we have a single gamma factor that approximates the CPD similar to the description in section 2.3.2. The corresponding integral looks like

[TABLE]

where the first gamma factor in eq. 97 represents the CPD approximation of the source dataset with $Q=\frac{1}{k_{mc}}\frac{(\sum_{j}w_{j})^{2}}{\sum_{j}w_{j}^{2}}$, $\beta=\frac{\sum_{j}w_{j}}{\sum_{j}w_{j}^{2}}$, and the second gamma factor marginalizes over the unknown rate where we leave a scaling parameter $\gamma$ for generality. We will later set $\gamma=1$. Additionally we define $A=\frac{\beta}{1+\beta}$ and $B=\frac{1}{1+\beta}$. We can also use the fact that $\frac{\Gamma(x+k)}{\Gamma(x)}=\sum\limits_{n=0}^{k}{k\brack n}x^{n}$ Graham1989 where ${k\brack n}$ are the unsigned Stirling numbers of the first kind and reformulate eq. 98 as

[TABLE]

Appendix C Overview of some generic formulas

The solutions of generalization (1) 4.1 can be written in various forms which we summarize in the following. We abbreviate the Poisson distribution as $\mathrm{P}$, the gamma distribution as $\mathrm{G}$, the gamma-Poisson mixture as $\mathrm{GP}$, and the Negative-Binomial distribution as $\mathrm{NB}$. When we marginalize the shape parameter of a gamma distribution with a $\mathrm{GP}$ distribution, we abbreviate that as $\mathrm{GPG}$, and a marginalization of a $\mathrm{GP}$ distribution shape parameter with another $\mathrm{GP}$ distribution is abbreviated as $\mathrm{GPGP}$.

C.1 Poisson-gamma mixture form

C.1.1 Simple generalization

The PMF of the generalized Poisson-gamma mixture has already derived in Gluesenkamp2018 . In Vellaisamy2009 it was shown that the convolution of Poisson-gamma mixtures is actually describing the same distribution, and so we have

[TABLE]

with iterative definition

[TABLE]

The combinatorial form (eq. 105) has been derived in Gluesenkamp2018 and in Vellaisamy2009 in the different context of a convolution of Poisson-gamma mixtures. However, the calculation based on the iterative sum (eq. 25) is much more efficient. The convolution of Poisson-gamma mixtures indicated in eq. 103 is an important step in the efficient calculation derived in appendix B. This distribution is very rich in structure, and its relationship to the Carlson-R function Carlson1963 and Lauricella function $F_{D}$ Lauricella1893 have been explored in Gluesenkamp2018 . The iterative formula has substantial computational advantages in various applications of these special functions.

C.1.2 Extended generalization

In appendix A we derived a generalization of the extended Poisson gamma mixture by taking a further expectation value (eq. 86) with respect another Poisson-gamma mixture. Here, we proceed more generally and only take the expectation value of $M$ Poisson-gamma terms while $N$ standard terms remain. Starting with the combinatorial expression similar to eq. 86 we can write

[TABLE]

[TABLE]

where we have have introduced additional M-dimensional vectors $\bm{\gamma}$, $\bm{\delta}$, $\bm{\varepsilon}$ that parametrize the additional $M$ marginalized factors resulting in the $\mathrm{PGPG}$ or $\mathrm{GPG}$ distributions, depending how we write it. If we put $M=N$, $\delta_{i}=\langle k_{mc,i}\rangle$, $\gamma_{i}=\beta_{i}=1/w_{i}$ and $\varepsilon_{i}=1$ we obtain eq. A.2 as a special case.

C.2 Negative Binomial form

Poisson-gamma mixtures are often re-parametrized as negative binomial distributions, which can be obtained by setting $\beta=\frac{1-p}{p}$ and using $\alpha=r$. For completion we give the simple and extended negative-binomial generalizations below.

C.2.1 Simple generalization

[TABLE]

with iterative definition

[TABLE]

If all $r_{i}$ expect one are zero, which effectively means there is only one pair of $r_{i}$ and $p_{i}$, we obtain the standard negative binomial distribution as a special case.

C.2.2 Extended generalization

In order to obtain an extended negative binomial distribution that is sensible, we take eq. C.1.2 and set $\beta_{i}=\frac{1-p_{i}}{p_{i}}$, $\gamma_{i}=\frac{1-p^{*}_{i}}{p^{*}_{i}}$, $\alpha_{i}=r_{i}$, $\delta_{i}=r^{*}_{i}$ and $\varepsilon_{i}=1$. The result gives more variance to the gamma factors with $p^{*}_{i}$ and $r^{*}_{i}$ which are marginalized with an extra Poisson-gamma factor. Using unity for $\varepsilon_{i}$ is the natural choice that is also used for the application in eq. A.2. We obtain

[TABLE]

with iterative definition

[TABLE]

Appendix D Supplementary figures