Poisson Quantile function returns wrong results in certain cases

[Caleb.Welton@…]

Sample use case:

main() {
  double pi = 3.1415927;
  boost::math::poisson poisson(pi);

  //  std::cout << quantile(poisson, 0.00) << " ";  
  std::cout << quantile(poisson, 0.05) << " ";
  std::cout << quantile(poisson, 0.30) << " ";
  std::cout << quantile(poisson, pi/6) << " ";
  std::cout << quantile(poisson, 0.7) << " ";
  std::cout << quantile(poisson, 1-1e-10) << " ";
  // std:: cout << quantile(poisson, 1) << std::endl;  
}

Produces the results (error), 0, 1, 3, 4, 20, (error)

Similar code in R or MATLAB results in: (0),1,2,3,4,20,(Inf)

R> qpois(c(0,0.05,0.3,pi/6,0.7,1-1e-10,1),pi)
[1]   0   1   2   3   4  20 Inf

MATLAB> f := stats::poissonQuantile(PI):
MATLAB> f(0), f(1/20), f(0.3), f(PI/6), f(0.7), f(1-1/10^10), f(1)
0,1,2,3,4,20,Inf

The errors on the edge cases are less troubling than the seemingly incorrect results at low probabilities.

[John Maddock]

I'll look into the edge cases (0 and 1) shortly - the differences in the other values you see are intended behaviour, given:

poisson pd(constants::pi<double>());

then the default behaviour of:

quantile(pd, 0.05)

is to return the largest integer i such that cdf(pd, i) <= 0.05, which really is the value 0, not 1 as the other programs suggest.

Whether that is actually the correct thing to do depends upon the actual test you are conducting (I assume R and Mathlab calculate the quantile as a real-value and then round to the nearest int). You can configure Boost.Math to do whatever you want here, see: http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/policy/pol_tutorial/understand_dis_quant.html

[John Maddock]

(In [83512]) Fix edge case in Halley iteration where the derivative is flatlining. Fix edge cases in poisson distribution quantile. Fixes #8314. Fixes #8308. Update and regenerate docs.