quantile for poisson distribution - precision issue?

[tan@…]

The following (there are others) returns an "unexpected" result:

#define BOOST_MATH_DISCRETE_QUANTILE_POLICY integer_round_up #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real

#include <cstdlib> #include <cstdio>

#include <boost/math/distributions/poisson.hpp>

int main() {

double m=54.3; double x=66.0;

boost::math::poisson_distribution<> dist( m );

std::printf(

"\nqpois(ppois(66,54.3),54.3) = %d\n\n", boost::math::quantile(

dist, boost::math::cdf( dist, x ) ) );

return 0

}

The "expected" answer is 66 but 67 is being returned.

I am aware of Ticket 8308 (​https://svn.boost.org/trac/boost/ticket/8308). I do not think R is "rounding to the nearest", and as far as I am aware, it has been consistent with

quantile(p)=infimum {x|cdf(x)>=p}

which is what "we" needed/expected.

If we change the policy to give real output---on my Ubuntu 12.04 (Intel 14.0 C++ with g++ 4.7.3), the result is 6.600000000000001e+01. So I wonder if this is just a (finite) precision/rounding issue ... in that it is really 66 to some "fuzz".

I have (also) tried this with revision 86054 of the trunk.

Thank you!

[anonymous]

I'm not surprised that that doesn't work as you expected - as you said it's a rounding issue - the calculated cdf is necessarily an approximation, and if you then try to invert that the chances of getting exactly 66 (or some other integer) back are basically nil. I guess we could check to see if the real-valued result is very near to an integer and then double check to see if that integer still fulfills the specified requirements, but I'm betting that might not work in this case either...

[John Maddock]

Can you try the attached patch? I think it fixes this and gives somewhat more sane results in the corner cases. Note however that it's still not possible to round trip in the general case - particularly when the probability is very close to 1 (within a few epsilon) as in that situation there may be multiple values for the random variable which all give the same CDF value.

[anonymous]

Dear John,

Thank you very much for the fixes. Yes, things look good here (I will run through even more tests later).

By "in the general case", do you mean when x (in the example above) is real, or that the "round trip" is also "bad" even when x is integral? Only the latter case is important to us.

Best regards, HS

[anonymous]

By "in the general case", do you mean when x (in the example above) is real, or that the "round trip" is also "bad" even when x is integral? Only the latter case is important to us.

When x is integral and large there may be multiple (all integral) values of x which all yield the same cdf, this usually occurs when the cdf is a couple of epsilon less than 1, but I found a few cases for very small cdf values as well. There's nothing that can be done about those cases - other than moving to a type with more precision, or, in the case where the cdf is near one, using complements instead.

[John Maddock]

(In [86205]) Improve quantiles of discrete distributions to round trip integers more often. Fixes #9183.

[John Maddock]

(In [86259]) Apply patch for issue #9183. Refs #9183.