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

Update kolmogorov quantile newton ub to 1 #1002

Merged
merged 4 commits into from
Jul 31, 2023
Merged

Update kolmogorov quantile newton ub to 1 #1002

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
37ba5a4
change kolmogorov quantile Newton ub from max to 1
ryanelandt Jul 25, 2023
b612afa
remove no longer needed newton step size logic
ryanelandt Jul 25, 2023
6bf5125
add infinity skew test back
ryanelandt Jul 25, 2023
e7ada7f
do not use Inf for Newton Bounds
ryanelandt Jul 25, 2023
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
4 changes: 2 additions & 2 deletions include/boost/math/distributions/kolmogorov_smirnov.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -382,7 +382,7 @@ inline RealType quantile(const kolmogorov_smirnov_distribution<RealType, Policy>
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.

return tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), boost::math::tools::max_value<RealType>(), get_digits, m);
k, RealType(0), RealType(1), get_digits, m);
} // quantile

template <class RealType, class Policy>
Expand All @@ -407,7 +407,7 @@ inline RealType quantile(const complemented2_type<kolmogorov_smirnov_distributio

return tools::newton_raphson_iterate(
detail::kolmogorov_smirnov_complementary_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), boost::math::tools::max_value<RealType>(), get_digits, m);
k, RealType(0), RealType(1), get_digits, m);
} // quantile (complemented)

template <class RealType, class Policy>
Expand Down
4 changes: 2 additions & 2 deletions include/boost/math/distributions/skew_normal.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -676,8 +676,8 @@ namespace boost{ namespace math{

// refine the result by numerically searching the root of (p-cdf)

const RealType search_min = range(dist).first;
const RealType search_max = range(dist).second;
const RealType search_min = support(dist).first;
const RealType search_max = support(dist).second;

const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
Expand Down
9 changes: 1 addition & 8 deletions include/boost/math/tools/roots.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -275,14 +275,7 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t&
if (fabs(delta * 2) > fabs(delta2))
{
// Last two steps haven't converged.
T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
if ((result != 0) && (fabs(shift) > fabs(result)))
{
delta = sign(delta) * fabs(result) * 1.1f; // Protect against huge jumps!
//delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216
}
else
delta = shift;
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
Copy link
Collaborator

Choose a reason for hiding this comment

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

Wait, we may need to protect against huge jumps in other people's code right? We only know that this is now unnecessary for Komogorov-Smirnov.

Copy link
Contributor Author

@ryanelandt ryanelandt Jul 25, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

I looked into this. Here is the history...

#145 @jzmaddock opened an issue where skew-normal quantiles return non-finite values. This issue occurs because the min/max bounds for the 1D Newton search problem are -Inf/+Inf. Attempting to bisect the search point at either ends gets one stuck at either -Inf/+Inf.
91c193d @jzmaddock Fixed #145 by adding an early version of huge step protection. A test was also added to test_roots to prevent possible regressions.
3ab69d0 The test in test_roots was removed as part of adding the complex Newton solver, possibly due to a merge conflict oversight.
#422 @evanmiller added the kolmogorov-smirnov distribution. The min/max bounds for the 1D Newton search problem are 0/max_value(). This new functionality requires huge step protection because evaluating the functor at large values causes the program to hang.

Copy link
Contributor Author

@ryanelandt ryanelandt Jul 25, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

Huge step protection seems to cover up issues with Newton search bounds for quantile problems.

This PR attempts to resolve this with the following steps:

  1. Change the skew_normal quantile search to go from -max_value/+max_value instead of -Inf/+Inf.
  2. Change the upper search bound of the kolmogorov-smirov distribution to RealType(1). This fixes the hang issue, and causes no test issues, I don't know if it is right though.
  3. Add the test that the skew-normal quantile is finite back in

Copy link
Member

Choose a reason for hiding this comment

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

I think that changing the bounds to [0, 1] is correct per this table of sample values http://www.matf.bg.ac.rs/p/files/69-[Michael_J_Panik]_Advanced_Statistics_from_an_Elem(b-ok.xyz)-776-779.pdf since it doesn't look like we support D_n,m.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Apologies for being late to the party... it'll still be another week at least till I have any time.

Removing the large step protection is probably a bad idea, but I would have to go searching for a test case to tell you why ;)

These root finders are rather problematic to test because they are fundamentally quite resilient even to bugs/errors in their code because of the way they work. I worry about that, but I also want them to be resilient!

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Removing the large step protection is probably a bad idea, but I would have to go searching for a test case to tell you why ;)

There's a test case that requires it in the current code-base, although it isn't causing a test failure. Discussed further here: #1000 (comment).

// reset delta1/2 so we don't take this branch next time round:
delta1 = 3 * delta;
delta2 = 3 * delta;
Expand Down
4 changes: 4 additions & 0 deletions test/test_roots.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -654,6 +654,10 @@ BOOST_AUTO_TEST_CASE( test_main )

test_beta(0.1, "double");

// bug reports:
boost::math::skew_normal_distribution<> dist(2.0, 1.0, -2.5);
BOOST_CHECK(boost::math::isfinite(quantile(dist, 0.075)));

#if !defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && !defined(BOOST_NO_CXX11_LAMBDAS)
test_complex_newton<std::complex<float>>();
test_complex_newton<std::complex<double>>();
Expand Down