Quantile of Poisson distribution

Roofit and RooStats
1

Dear Rooters,

I am counting events that follow the Poisson distribution and I need to put a threshold to the number of counts. The threshold is defined as the value associated to the probability P that the number of counts C is bigger than the threshold; to my knowledge this translates into the quantile of the Poisson distribution.

I was given the following code for boost

boost::math::poisson_distribution<> pd(C); threshold = boost::math::quantile(complement(pd, P);
I would like to do it in ROOT, but in the documentation I couldn’t find a function named poisson_quantile or equivalent. I also searched the forum but I found no posts of immediate solution.

Could you help me? :blush:

2

It should work out of the box if boost is installed:

$ root
root [0] #include <boost/math/distributions/poisson.hpp> 
root [1] boost::math::poisson_distribution<> pd(42);
root [2] boost::math::quantile(complement(pd, 0.01))
(double) 58.000000
3

Hi,

You can also invert the existing poisson_cdf function in ROOT.
For example

TF1 f("f","ROOT::Math::poisson_cdf_c(x,42)",0,100);
int value = int( f.GetX(0.01) ) + 1;

Lorenzo

4

Hi everybody :slight_smile:

@behrenhoff: thanks, it works for me as well, at least when I compile. But I forgot to mention that we are trying to reduce the number of dependencies and I am pushing for keeping only ROOT for the math…

@moneta: apparently it works. But before it can be accepted I have to demonstrate that for the same values of C and P the results are consistent with the previous implementation, or explain the differences. I expect some discrepancies, let’s call them “rounding errors”, and I am not worried about that, but with ROOT I also get non-monotonic behavior in some circumstances. The upper boundary of the formula definition seems to be particularly important and I don’t understand why. You can see what I mean if you run this:

  double prob = 0.001;

  TCanvas c0("c0", "Comparison", 1280, 1024);
  TH1F h_boost("h_boost", "Boost", 200, 0.0, 200.0);
  TH1F h_croot("h_croot", "ROOT",  200, 0.0, 200.0);

  for (int i=1; i<201; ++i) {
    boost::math::poisson_distribution<> pd(i);
    TF1 f("f",Form("ROOT::Math::poisson_cdf_c(x,%d)",i), 0, 1000);
    h_boost.SetBinContent(i,boost::math::quantile(complement(pd, prob)));
    h_croot.SetBinContent(i,f.GetX(prob));
  }
  h_boost.Draw();
  h_croot.SetLineColor(2);
  h_croot.Draw("SAME");
  c0.Print("thresholds.png");

There are some values (around 50, 90, 115…) for which if I increase the counts by 1, my new threshold is lower and this makes no sense.

If I put 300 instead of 1000 in the definition of TF1 the plot looks way better… I think I am missing something… :question:

5

f.SetNpx(10000);

6

Thanks Pepe, that was the missing part :smiley: