# Quantile of Poisson distribution

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?

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``````

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

Hi everybody

@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…

f.SetNpx(10000);

Thanks Pepe, that was the missing part