Repository navigation
binomial_coefficient_log produces wrong results when N >> n #1592
Description
Activity
This also suggests that binomial and beta densities could be subject to imprecisions :-(
- Thanks for catching and debugging this. We should definitely fix. We can use any approximation that works.…On Jan 7, 2020, at 8:01 AM, Martin Modrák ***@***.***> wrote: Description For some values, the results from binomial_coefficient_log don't match those computed by R or Mathematica. In #239, it is noted that we don't know, where the approximation we use for the case when N >> n comes from. I also don't understand it. R uses a quite a different approximation - their implementation is double attribute_hidden lfastchoose(double n, double k) { return -log(n + 1.) - lbeta(n - k + 1., k + 1.); } Stirling approximation also produces the correct results in the cases I've found. Example In Stan: binomial_coefficient_log(1e20, 100) = 4141.4308104325282 binomial_coefficient_log(845000, 8.45e-10) = 1.1925051146203358e-008 Expected Output Compare to Mathematica, Rs lchoose https://www.wolframalpha.com/input/?i=log%2810%5E20+choose+100%29 Stan: 4141.4308104325282 Math: 4241.430810432527877842402915999072777157961556194658677330 R: 4241.4308104325273 https://www.wolframalpha.com/input/?i=log%28845000+choose+845*10%5E-12%29 Stan: 1.1925051146203358e-008 Math: 1.2019540397111151911538109607610600259536643522634943 × 10^-8 R: 0 #Not great Current Version: v3.0.0 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.
So I made a mistake in the initial description - Stirling approximation works reasonably only for some of the problematic cases (description now updated). Searching online, the links I found either say "use lgamma" (which we know is not enough) or "steal code from R", which we can't do as it is GPL. And the idea in R's code is (for me) non-trivial to adapt. The main workhorse in R is a function called
lgamma_corwhich computes the difference between Stirling's approximation tolgammaand actuallgammavalue. They then proceed by computing the Stirling's approximation forlbetaand adding the corrections separately. Unfortunately, the code forlgamma_coris opaque to me, so I can't reimplement it to avoid licensing issues. It uses something calledchebyshev_evalwhich is probably related to Chebyshev polynomials, which have their boost implementation but I am totally out of my depth here and don't understand the surrounding math at all :-(Also, if I read the code correctly, Stan now relies on autodiff for derivatives of
binomial_coefficient_log(I believe Stan'slchoosedelegates tobinomial_coefficient_login math and there is norevversion of the function. But autodiffing through a complex approximation (as the one in R) might be undesirable, so analytic derivatives might need to be developed and tested as well. We do have analytic derivatives forlbetaimplementation we have, but I would suspect they would fail where the current implementation fails.- linked a pull request that will close this issueImproved numerical stability of binomial_coefficient_log #1614
on Mar 13, 2020
Description
For some values, the results from
binomial_coefficient_logdon't match those computed by R or Mathematica.In #239, it is noted that we don't know, where the approximation we use for the case when
N >> ncomes from. I also don't understand it. R uses a quite a different approximation - their implementation isStirling approximation produce good results for some of the cases, but also terrible results for others... The R's implementation seems to work neatly, but depends on quite a bunch of other code to compute the
lbeta.Example
In Stan:
Expected Output
Compare to Mathematica, Rs
lchooseand Stirling approximation implemented in Rvia
In both cases, there is a mismatch in second digit.
https://www.wolframalpha.com/input/?i=log%2810%5E20+choose+100%29
https://www.wolframalpha.com/input/?i=log%28845000+choose+845*10%5E-12%29
Here R-cor is computed using the
lfastchooseformula directly, e.g.:for some reason,
lchoosecalls a different computation method for those parametersCurrent Version:
v3.0.0