Bug Fix: Inaccurate Tails of Incomplete Gamma #481
Conversation
akleeman
force-pushed
the
akleeman/incomplete_gamma_tails
branch
from
064989b to
87ef448
Compare
akleeman
force-pushed
the
akleeman/incomplete_gamma_tails
branch
from
87ef448 to
8c8dbaa
Compare
akleeman
force-pushed
the
akleeman/incomplete_gamma_tails
branch
from
8c8dbaa to
d8409e0
Compare
| } else { | ||
| dof += 1; | ||
| } | ||
| } |
There was a problem hiding this comment.
Any way you can incorporate this function into the test to automatically dump out the values for use in the comparisons inside the test_chi_squared_matches_gsl unit tests? Might make any future modifications/iterations easier
There was a problem hiding this comment.
The problem with that is that we'd then need to pull GSL in as a third party library and their license (GNU public) is too restrictive for albatross. I suppose it's possible they'd allow use of the library as a test/comparison not as part of the core algorithms ... but that seems like a pretty gray area.
There was a problem hiding this comment.
Ugh never mind then.. we can stick it here https://github.com/swift-nav/shared-scratch in the future if we ever need to tweak it in the future
There was a problem hiding this comment.
Can you add a note to that comment so we remember why we couldn't just add that output function there? Or just copy and paste into https://github.com/swift-nav/shared-scratch and add a link from orion-engine there
There was a problem hiding this comment.
oops, sorry already merged it. but I'll follow up.
There was a problem hiding this comment.
follow up PR for this is fine
There was a problem hiding this comment.
| // 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))); |
There was a problem hiding this comment.
Were these bounds experimentally determined?
There was a problem hiding this comment.
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.
This PR fixes a bug in which the bounds to the incomplete gamma function quadrature were not sufficiently large to properly evaluate the tails. As a result evaluations of the chi squared cdf for very large values would (for specific degrees of freedom) not asymptotically approach 1.
To evaluate the discrepancy I used the GSL implementation of the incomplete gamma function to provide an alternative
chi_squared_cdfmethod. Then compared the values to the one implemented in albatross. Here is the error (y-axis) as a function of the expected value (from GSL). Notice both x and y axes are log scale:The differences were not too large, mostly only ever reaching
0.001. BUT, this can have serious consequences! Consider outlier detection, we often use the chi squared cdf as a threshold for marking an inlier. Usually that threshold is really close to 1, something like1 - 1e-6, which can be interpreted as "throw out anything which looks like a one in a million event. Notice in the plot above that the largest errors are actually as you approach 1. What this means is that when we've given a measurement which is really bad, the chi squared cdf may only ever get as large as1 - 1e-4, so even for insanely large residuals we may never throw out an outlier!To fix this I was able to increase the integration bounds used in the quadrature versions of the incomplete gamma:
There's still a strange blip near 1. In those cases the chi squared percentile we compute is typically larger than GSL (by at most 3e-7). So in those cases we should continue to throw out bad values, even for inlier thresholds around
0.99999999