Tackling limited simulation and small signals

TL;DR

This paper introduces an analytic Poisson likelihood method to better handle statistical uncertainties in small simulation samples, improving coverage and computational efficiency.

Contribution

The paper presents a novel analytic Poisson likelihood technique that improves statistical uncertainty modeling in limited simulation data.

Findings

Better coverage properties than existing methods

Valid for small data samples

Maintains good computational performance

Abstract

We present a new, analytic, Poisson likelihood derived, technique to account for the statistical uncertainties inherent in simulation samples of limited size. This method has better coverage properties than other techniques, is valid for small data samples, and maintains good computational performance.

Equations

$$\mathcal{L}(\vec{\theta }|k)=\mathrm{Poisson}(\mathrm{k};\lambda (\vec{\theta }))=\frac{\lambda (\vec{\theta }{)}^{\mathrm{k}}{\mathrm{e}}^{-\lambda (\vec{\theta })}}{\mathrm{k}!},$$

$${\mathcal{L}}_{\text{AdHoc}}(\vec{\theta }|k)=\frac{{\left({\sum }_i{w}_i(\vec{\theta })\right)}^k{e}^{-\left({\sum }_i{w}_i\left(\vec{\theta }\right)\right)}}{k!}.$$

$${\mathcal{L}}_{\mathrm{General}}(\vec{\theta }|k)={\int }_0^{\infty }\frac{{\lambda }^k{e}^{-\lambda }}{k!}\mathcal{P}\left(\lambda |\vec{w}(\vec{\theta })\right)d\lambda .$$

$$\mathcal{P}\left(\lambda |\vec{w}(\vec{\theta })\right)=\frac{\mathcal{L}(\lambda |\vec{w}(\vec{\theta }))\mathcal{P}(\lambda )}{{\int }_0^{\infty }\mathcal{L}({\lambda }^{\prime }|\vec{w}(\vec{\theta }))\mathcal{P}({\lambda }^{\prime })d{\lambda }^{\prime }},$$

$$\mu \equiv \sum _{i=1}^m{w}_i\text{and}{\sigma }^2\equiv \sum _{i=1}^m{w}_i^2,$$

$$\mu =wm\text{,}{\sigma }^2=w^{2}m\text{,}w={\sigma }^2/\mu \text{, and }m={\mu }^2/{\sigma }^2.$$

$$\mathrm{Poisson}(M=m;\bar{m})=\frac{e^{-\bar{m}}{\bar{m}}^m}{m!}.$$

$$\mathcal{L}(\lambda |\vec{w}(\vec{\theta }))=\mathcal{L}(\lambda |\mu ,\sigma )=\frac{e^{-\lambda \mu /{\sigma }^2}{\left(\lambda \mu /{\sigma }^2\right)}^{{\mu }^2/{\sigma }^2}}{({\mu }^2/{\sigma }^2)!}.$$

$$\mu =w_{\mathrm{Eff}}{m}_{\mathrm{Eff}}\text{and}{\sigma }^2=w_{\mathrm{Eff}}^2{m}_{\mathrm{Eff}}.$$

$$\mathcal{L}(\bar{m}|m_{\mathrm{Eff}})$$

$$\mathcal{L}(\lambda |\vec{w}(\vec{\theta }))=\mathcal{L}(\lambda |\mu ,\sigma )=\frac{e^{-\lambda \mu /{\sigma }^2}{\left(\lambda \mu /{\sigma }^2\right)}^{{\mu }^2/{\sigma }^2}}{\Gamma ({\mu }^2/{\sigma }^2+1)}.$$

$$\alpha \equiv \frac{{\mu }^2}{{\sigma }^2}+1\text{and}\beta \equiv \frac{\mu }{{\sigma }^2}.$$

$$\mathcal{P}(\lambda |\vec{w}(\vec{\theta }))$$

$$=\frac{e^{-\lambda \beta }{\lambda }^{\alpha -1}{\beta }^{\alpha }}{\Gamma (\alpha )}$$

$$=\mathcal{G}(\lambda ;\alpha ,\beta ),$$

$${\mathcal{L}}_{\text{Eff}}(\vec{\theta }|k)$$

$$=\frac{{\beta }^{\alpha }\Gamma \left(k+\alpha \right)}{k!{\left(1+\beta \right)}^{k+\alpha }\Gamma \left(\alpha \right)}$$

$$={\left(\frac{\mu }{{\sigma }^2}\right)}^{\frac{{\mu }^2}{{\sigma }^2}+1}\Gamma \left(k+\frac{{\mu }^2}{{\sigma }^2}+1\right){\left[k!{\left(1+\frac{\mu }{{\sigma }^2}\right)}^{k+\frac{{\mu }^2}{{\sigma }^2}+1}\Gamma \left(\frac{{\mu }^2}{{\sigma }^2}+1\right)\right]}^{-1}.$$

Full text

For binned data, the Poisson likelihood is taken to be the probability to observe events in a bin and is commonly used in high-energy physics and particle astrophysics experiments. Given an exact expectation rate, $\lambda(\vec{\theta})$, the probability of observing an integer $k$ events is

[TABLE]

where $\vec{\theta}$ is some set of physics parameters that determine $\lambda$. Given the stochastic nature of processes in particle physics, exactly determining $\lambda$ is often not possible and requires Monte Carlo (MC) simulations. Simulation is often expensive, and so reweighting is employed as it enables a single simulation to be used to describe many physical hypotheses [1].

In such scenarios, an ad hoc likelihood is commonly used, where $\lambda$ is simply taken to be the sum of the weights in each bin, namely

[TABLE]

A notable downside of this ad hoc likelihood is that it neglects the statistical uncertainty inherent in estimating $\lambda(\vec{\theta})$ from a simulation of limited size. For expensive simulations, or physical hypotheses far from the original simulation, this uncertainty can be non-negligible [2, 3, 4, 5]. One can account for this uncertainty by incorporating a simulation derived prior on $\lambda$, denoted as $\mathcal{P}\left(\lambda|\vec{w}(\vec{\theta})\right)$, that has non-zero variance. Thus, we can write the likelihood as the marginalization of the Poisson likelihood with the prior,

[TABLE]

We construct this prior based on the likelihood function of the simulation outcome and a prior on $\lambda$, $\mathcal{P}(\lambda)$,

[TABLE]

where in our implementation we have chosen we chosen $\mathcal{P}(\lambda)$ to be uniform.

Let us first consider the case where all simulation events in the bin have equal weight. In this scenario we can relate the number of events, $m$, and the weight of the events, $w$, to the quantities $\mu$ and $\sigma$, which are defined

[TABLE]

and satisfy the relationships

[TABLE]

The probability of obtaining $m$ events in the simulation bin can be modelled with the Poisson distribution assuming the true but unknown mean $\bar{m}$:

[TABLE]

This allows us to rewrite the likelihood of $\lambda$ in terms of $\mu$ and $\sigma$ as

[TABLE]

If the simulation event weights are not all equal, as is usually the case, then we can replace $w$ and $m$ with their “effective” counterparts $w_{\mathrm{Eff}}$ and $m_{\mathrm{Eff}}$. These then relate to $\mu$ and $\sigma$ as

[TABLE]

The replacement redefines the likelihood of our simulation outcome

[TABLE]

which, assuming $\lambda=w_{\mathrm{Eff}}{\bar{m}}$, can be rewritten as

[TABLE]

To simplify the notation, define

[TABLE]

Substituting Eq. (11) into Eq. (4) and assuming a uniform $\mathcal{P}(\lambda)$, we obtain

[TABLE]

where $\mathcal{G}(\lambda;\alpha,\beta)$ is the gamma distribution, with shape and inverse rate parameters $\alpha$ and $\beta$. Finally, this can be substituted for $\mathcal{P}\left(\lambda|\vec{w}(\vec{\theta})\right)$ in Eq. (3) so that

[TABLE]

Equation (16) is an effective likelihood, motivated by Poisson statistics, and derived with a Bayesian approach. It incorporates statistical uncertainties inherent in the MC approximation of the rate by encoding the distribution of weights in terms of $\mu$ and $\sigma^{2}$. The effective likelihood $\mathcal{L}_{\textmd{Eff}}$ can be easily substituted for $\mathcal{L}_{\textmd{AdHoc}}$. A more thorough exposition, along with a generalization for different priors, $\mathcal{P}(\lambda)$, is given in [6].