Simple Poissonian log-likelihood (`llh`) test statistic is inaccurate (undefined) for low (zero) observed counts

[thehrh]

... because of the use of (an approximation of) Stirling's approximation.

Ever since 442a7ef#diff-f7283bae353d06044b8ff9fe43c9fcc6fd92165586395923983197fc9871293e, utils.stats.llh has for some reason been approximating $\ln(n!) = n\text{ }\ln(n) - n$, which is undefined for $n=0$ and wildly inaccurate for low observed counts.

Here are some deviations between the left-hand side and our approximation on the right for various $n$:

• $n=0$: 0 vs. nan
• $n=1$: 0 vs. -1
• $n=5$: 4.79 vs. 3.05
• $n=10$: 15.10 vs. 13.03
• $n=50$: 148.5 vs. 145.6

Before that change, we had been using $\ln(\Gamma(n+1))$, which agrees with $\ln(n!)$ by definition for integer $n$.

Whether this leads to a noticeable deviation of the assumed log-likelihood expression from a Poissonian one or even gives rise to pseudoexperiments with nan llh values depends on how fluctuations are produced and how large these are with respect to the expectations of course.

Any ideas why we stopped using gammaln (https://en.wikipedia.org/wiki/Poisson_distribution#Evaluating_the_Poisson_distribution)?

In any case, I think we should enable the use of an (actually Poissonian) log-likelihood metric, if need be with a new name so that existing analyses that have been using llh as it stands today can still be reproduced.

[thehrh]

To illuminate this issue further, what the current implementation of utils.stats.llh allows is to skip computing the value at the truth hypothesis underlying an Asimov pseudoexperiment when evaluating the log-likelihood difference, see below. I suspect that this could have served as a motivation for adopting the approximation.

Consider the correct Poissonian log-likelihood for a pseudoexperiment with integer $n$ (one bin). The log-likelihood of the observation given any hypothesis, $\mathcal{L}^\prime$, is then defined as

$$\mathcal{L}^\prime = -\mu^\prime + n \ln(\mu^\prime) - \ln(n!).$$

Assume the value at the global maximum is

$$\mathcal{L} = -\mu + n \ln(\mu) - \ln(n!).$$

The difference is then

$$\mathcal{L} - \mathcal{L}^\prime = \mu^\prime - \mu + n\left(\ln(\mu) - \ln(\mu^\prime)\right) + \ln(n!) - \ln(n!).$$

Clearly, if we evaluated the log-likelihood at the maximum and somewhere else, the last two terms, which only depend on the observation, cancel, whether the the approximation is used or not. However, with the approximation, when $n=0$, evaluating the log-likelihood for any hypothesis would cause trouble.

Now consider the Asimov case, with $n=\mu$ (observation = expectation under some truth, which also defines the likelihood maximum). Here, from the last expression above, the exact log-likelihood difference would read

$$\mathcal{L} - \mathcal{L}^\prime = \mu^\prime - \mu + \mu\left(\ln(\mu) - \ln(\mu^\prime)\right).$$

Under the approximation ($\ln(n!)$ becomes $\mu\ln(\mu) - \mu$) we have

$$\mathcal{L}^\prime = -\mu^\prime + \mu\ln(\mu^\prime) - \mu\ln(\mu) + \mu = -\mu^\prime + \mu - \mu\left(\ln(\mu) - \ln(\mu^\prime\right) = -\left(\mathcal{L} - \mathcal{L}^\prime\right),$$

the negative of the exact log-likelihood difference!

In summary:

• the current implementation results in only approximate log-likelihood values $\mathcal{L}^\prime$ for pseudoexperiments (as per top comment), and in nan when an observation $n=0$
• the current implementation results in exact log-likelihood differences between any two hypotheses for pseudoexperiments, unless an observation $n=0$, in which case each term is nan
• if you do not observe fewer log-likelihood values in your histogram than there were pseudoexperiments generated, I believe no pseudoexperiment should have had a vanishing bin count
  • depending on how small the observations actually are, it could just be that the true Poissonian log-likelihood distribution would look different; after all, the (approximation of) Stirling's approximation underestimates $ln(n!)$ for small $n$
  • TODO: check if best-fit values are indeed unaffected due to deviation being restricted to term that is exclusively observation-dependent; check if true Poissonian log-likelihood is indeed overestimated
• for an Asimov pseudoexperiment, the current implementation results in the single log-likelihood value $\mathcal{L}^\prime$ (given any hypothesis) matching the negative of the exact log-likelihood difference between the truth and that hypothesis

NB: expectations $\mu=0$ or $\mu^\prime=0$ are always replaced by some small positive value SMALL_POS = 1e-10.