Repository navigation
lbeta is unstable when the arguments differ a lot #1611
Description
Activity
Nice catch.
Is there a known fix here? We're currently using this:
lgamma(a) + lgamma(b) - lgamma(a + b);Fixed along #1592 in #1497 (if it is preferable, I might split this fix from the #1497 PR into a separate one , but so far it is IMHO all related as failures in
neg_binomial_2_lpmfmotivated the improvements forbinomial_coefficient_logandlbeta).It is a neat trick I found in the R codebase - you first compute the
lbetavia Stirling approximation tolgammawhich allows you to analytically manipulate the expression and avoid the cancellation for the "bulk" of the value. Then you compute the difference betweenlgammaand the Stirling approximation and add this to the result to restore precision.The code in R for the difference between
lgammaand Stirling was opaque and used some tricks I didn't really get, but it turns out just using 3 additional terms of the Stirling series as found on Wikipedia is enough for all the tests I could come up with :-)Assuming it doesn't take an inordinate amout of effort, having it separate makes it easier to review, and will help if in the future this needs to be revisited for whatever reason.
- Thanks for the explanation. This numerics stuff is tricky.…On Jan 13, 2020, at 2:15 PM, Martin Modrák ***@***.***> wrote: Fixed along #1592 in #1497 (if it is preferable, I might split this fix from the #1497 PR into a separate one , but so far it is IMHO all related as failures in neg_binomial_2_lpmf motivated the improvements for binomial_coefficient_log and lbeta). It is a neat trick I found in the R codebase - you first compute the lbeta via Stirling approximation to lgamma which allows you to analytically manipulate the expression and avoid the cancellation for the "bulk" of the value. Then you compute the difference between lgamma and the Stirling approximation and add this to the result to restore precision. The code in R for the difference between lgamma and Stirling was opaque and used some tricks I didn't really get, but it turns out just using 3 additional terms of the Stirling series as found on Wikipedia is enough for all the tests I could come up with :-) — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.
So, I separated this into PR #1612, should be ready for review.
Description
When the arguments are on wildly different scales,
lbetasuffers from catastrophic cancellation.Example + Expected output
https://www.wolframalpha.com/input/?i=LogGamma%2810%5E20%29+%2B+LogGamma%283%29+-+LogGamma%2810%5E20+%2B+3%29
https://www.wolframalpha.com/input/?i=LogGamma%283*10%5E15%29+%2B+LogGamma%2812895%29+-+LogGamma%283*10%5E15+%2B+12895%29
Current Version:
v3.0.0