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

Bug Fix: Inaccurate Tails of Incomplete Gamma #481

Merged
merged 4 commits into from
Apr 29, 2024
Merged

Bug Fix: Inaccurate Tails of Incomplete Gamma #481

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
e694d5f
Add thorough test of chi_squared_cdf
akleeman Apr 23, 2024
1101bd1
Allow for non-integer degrees of freedom
akleeman Apr 23, 2024
4ea5866
Increase the bounds of quadrature
akleeman Apr 23, 2024
d8409e0
Test new bounds
akleeman Apr 23, 2024
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
1 change: 1 addition & 0 deletions BUILD.bazel
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,7 @@ swift_cc_test(
"tests/test_block_utils.cc",
"tests/test_call_trace.cc",
"tests/test_callers.cc",
"tests/test_chi_squared_versus_gsl.cc",
"tests/test_compression.cc",
"tests/test_concatenate.cc",
"tests/test_conditional_gaussian.cc",
Expand Down
1 change: 1 addition & 0 deletions include/albatross/Stats
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#ifndef ALBATROSS_STATS_H
#define ALBATROSS_STATS_H

#include <albatross/src/details/typecast.hpp>
#include <albatross/src/stats/gaussian.hpp>
#include <albatross/src/stats/gauss_legendre.hpp>
#include <albatross/src/stats/incomplete_gamma.hpp>
Expand Down
17 changes: 13 additions & 4 deletions include/albatross/src/stats/chi_squared.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,11 @@ namespace albatross {

namespace details {

inline double chi_squared_cdf_unsafe(double x, std::size_t degrees_of_freedom) {
return incomplete_gamma(0.5 * cast::to_double(degrees_of_freedom), 0.5 * x);
inline double chi_squared_cdf_unsafe(double x, double degrees_of_freedom) {
return incomplete_gamma(0.5 * degrees_of_freedom, 0.5 * x);
}

inline double chi_squared_cdf_safe(double x, std::size_t degrees_of_freedom) {
inline double chi_squared_cdf_safe(double x, double degrees_of_freedom) {

if (std::isnan(x) || x < 0.) {
return NAN;
Expand All @@ -53,7 +53,16 @@ inline double chi_squared_cdf_safe(double x, std::size_t degrees_of_freedom) {

} // namespace details

inline double chi_squared_cdf(double x, std::size_t degrees_of_freedom) {
template <typename IntType,
typename = std::enable_if_t<std::is_integral<IntType>::value>>
inline double chi_squared_cdf(double x, IntType degrees_of_freedom) {
// due to implicit argument conversions we can't directly use cast::to_double
// here.
return details::chi_squared_cdf_safe(x,
static_cast<double>(degrees_of_freedom));
}

inline double chi_squared_cdf(double x, double degrees_of_freedom) {
return details::chi_squared_cdf_safe(x, degrees_of_freedom);
}

Expand Down
31 changes: 5 additions & 26 deletions include/albatross/src/stats/incomplete_gamma.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,32 +75,11 @@ inline double incomplete_gamma_quadrature_recursive(double lb, double ub,

inline std::pair<double, double> incomplete_gamma_quadrature_bounds(double a,
double z) {

if (a > 800) {
return std::make_pair(std::max(0., std::min(z, a) - 11 * sqrt(a)),
std::min(z, a + 10 * sqrt(a)));
} else if (a > 300) {
return std::make_pair(std::max(0., std::min(z, a) - 10 * sqrt(a)),
std::min(z, a + 9 * sqrt(a)));
} else if (a > 90) {
return std::make_pair(std::max(0., std::min(z, a) - 9 * sqrt(a)),
std::min(z, a + 8 * sqrt(a)));
} else if (a > 70) {
return std::make_pair(std::max(0., std::min(z, a) - 8 * sqrt(a)),
std::min(z, a + 7 * sqrt(a)));
} else if (a > 50) {
return std::make_pair(std::max(0., std::min(z, a) - 7 * sqrt(a)),
std::min(z, a + 6 * sqrt(a)));
} else if (a > 40) {
return std::make_pair(std::max(0., std::min(z, a) - 6 * sqrt(a)),
std::min(z, a + 5 * sqrt(a)));
} else if (a > 30) {
return std::make_pair(std::max(0., std::min(z, a) - 5 * sqrt(a)),
std::min(z, a + 4 * sqrt(a)));
} else {
return std::make_pair(std::max(0., std::min(z, a) - 4 * sqrt(a)),
std::min(z, a + 4 * sqrt(a)));
}
// NOTE: GCEM uses a large conditional block to select tighter bounds, but in
// practice those bounds were not tight enough, particularly on the upper
// bound, so we've modified this function to be more conservative
return std::make_pair(std::max(0., std::min(z, a) - 12 * sqrt(a)),
std::min(z, a + 13 * sqrt(a + 1)));
Copy link
Contributor

Choose a reason for hiding this comment

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

Were these bounds experimentally determined?

Copy link
Collaborator Author

@akleeman akleeman Apr 29, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

Yeah, I basically just copied the format from the master bounds, but simplified it instead of changing the multipliers (12/13 now) based on the value of a I just set them higher than any of the previous if statements. I didn't exactly look for the lowest bound possible (so it's possible there's some wasted CPU here?) but I did increase them each until all the new tests I added passed, so they shouldn't be too crazy.

}

inline double incomplete_gamma_quadrature(double a, double z) {
Expand Down
Loading
Loading