Poisson interval calculation with asymptotic approximation

[AlkaidCheng]

Check duplicate issues.

• Checked for duplicates

Description

In the source code for computing the Poisson interval (used in histograms I believe). First there is a hard-coded threshold of n = 100 where the asymptotic approximation is used (which was not explained clearly in the documentation as well as where the formula come from). The formula seems to be only valid if nSigma=1 (https://root.cern.ch/doc/v626/RooHistError_8cxx_source.html#l00112):

  // use assymptotic error if possible
  if(n > 100) {
    mu1= n - sqrt(n+0.25) + 0.5;
    mu2= n + sqrt(n+0.25) + 0.5;
    return kTRUE;
  }

I believe the formula is derived by inverting the Score Test for the Poisson distribution. The general formula is probabaly:
$$n + \dfrac{Z^2}{2} \pm Z \sqrt{n +\dfrac{Z^2}{4}}$$

Reproducer

>> import ctypes
>> import ROOT
>> 
>> inst = ROOT.RooHistError.instance()
>> ym = ctypes.c_double()
>> yp = ctypes.c_double()
>> nSigma = 3
>> n_list = [100, 101]
>> errlo = {}
>> errhi = {}
>> for nSigma in [1, 2, 3]:
>>     print(f'{nSigma} Sigma')
>>     errlo = []
>>     errhi = []
>>     for n in [100, 101]:
>>         inst.getPoissonInterval(n, ym, yp, nSigma)
>>         errlo.append(n - ym.value)
>>         errhi.append(yp.value - n)
>>     print('ErrorLo [N = 100, N = 101] : ', errlo)
>>     print('ErrorHi [N = 100, N = 101] : ', errhi)
1 Sigma
ErrorLo [N = 100, N = 101] :  [9.983254822886579, 9.56230589874906]
ErrorHi [N = 100, N = 101] :  [11.033360941149681, 10.56230589874906]
2 Sigma
ErrorLo [N = 100, N = 101] :  [18.98411465660803, 9.56230589874906]
ErrorHi [N = 100, N = 101] :  [22.082468714289504, 10.56230589874906]
3 Sigma
ErrorLo [N = 100, N = 101] :  [27.35378945919335, 9.56230589874906]
ErrorHi [N = 100, N = 101] :  [33.829521810749526, 10.56230589874906]

ROOT version

6.38.00

Installation method

conda

Operating system

Linux

Additional context

No response

[guitargeek]

Thank you for the report. Indeed the approximation is not documented and the behavior for sigma != 1 and n > 100 is wrong.

I don't think we should fix the approximation though. Expressing exact Poisson intervals in terms of the chi-square distribution is easy, I think we should just do that 🙂

PR to fix this issue:

• [RF] Replace asymptotic Poisson interval with exact Garwood intervals #21443

You have any comments on the approach?

[AlkaidCheng]

Thanks @guitargeek for the quick response! If we can use a consistent approach that would be ideal and the exact Garwood interval seems the best candidate here. Though I have read that exact approach (Garwood interval) can be a little conservative for large N (Agresti, A., & Coull, B. A. (1998): "Approximate is Better than 'Exact' for Interval Estimation of Binomial Proportions.").

In terms of the implementation, how fast is the new method of computing the Poisson interval (with gamma functions I believe)? It might still be good to have some sort of caching when we want to compute a large number of intervals. If the overhead is minimal, we can just go without caching.

[guitargeek]

It's true that there are different intervals too. And interesting paper that you referenced! I found the pdf:
https://math.unm.edu/~james/Agresti1998.pdf

Figure 5 is pretty good to illustrate the behavior change! Before, RooFit was using Exact intervals for n < 100, and then the Score intervals beyond. This undocumented change of convention at an arbitrary point is not good, so now we are using exact intervals all the way. The other alternative, using score intervals for all n, was not acceptable because you get significant under-coverage for small n, and generally we don't accept under-coverage in statistical analysis.

About the caching: the new implementation is 10x faster, so at this point the caching can't be justified without a clear motivation. Caching should not be done lightly, because it makes memory usage patterns less transparent to the user. Also, the caching of interval values was introduced over 20 years ago, and performance characteristics surely also have changed since then 🙂