Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix inverse_discrete_quantile for large guess #1007

Merged
merged 1 commit into from
Aug 24, 2023
Merged

Fix inverse_discrete_quantile for large guess #1007

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 6 additions & 10 deletions include/boost/math/distributions/detail/inv_discrete_quantile.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -302,15 +302,13 @@ inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::va
//
while(result != 0)
{
cc = result - 1;
cc = floor(float_prior(result));
if(cc < support(d).first)
break;
pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
if(pp == p)
result = cc;
else if(c ? pp > p : pp < p)
if(c ? pp > p : pp < p)
break;
result -= 1;
result = cc;
}

return result;
Expand All @@ -336,15 +334,13 @@ inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::val
//
while(true)
{
cc = result + 1;
cc = ceil(float_next(result));
if(cc > support(d).second)
break;
pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
if(pp == p)
result = cc;
else if(c ? pp < p : pp > p)
if(c ? pp < p : pp > p)
break;
result += 1;
result = cc;
}

return result;
Expand Down
12 changes: 12 additions & 0 deletions test/test_binomial.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -716,6 +716,18 @@ void test_spots(RealType T)

check_out_of_range<boost::math::binomial_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.

// TODO: Generic ibeta_power_terms has accuracy issue when long
// double is not precise enough, causing overflow in this case.
if(!std::is_same<RealType, real_concept>::value ||
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happens here when T is a double, does that still work? Ideally we shouldn't be having spurious internal overflow as a) it breaks types with no infinity (this is potentially true for long double too) and b) it sets the hardware overflow flag which can trip up downstream consumers.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The problem is in beta.hpp#478 onwards. Don't see a trivial fix just yet :(

inv_distrete_quantile.hpp#150 also needs adjusting to else if((adder != 0) && (a + adder != a)) to speed up this new test case.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With built-in double and long double it works fine because Lanczos approximation is enabled.

sizeof(boost::math::concepts::real_concept_base_type) > 8)
{
using namespace boost::math::policies;
typedef policy<discrete_quantile<integer_round_outwards> > Policy;
binomial_distribution<RealType, Policy> dist(9079765771874083840, 0.561815);
// Accuracy is not too important here; the main purpose is to
// make sure it is not stuck.
BOOST_CHECK_CLOSE(quantile(dist, 0.0365346), 5101148604445670400, 1e12);
}

} // template <class RealType>void test_spots(RealType)

Expand Down