Repository navigation
binomial_coefficient_log producing wrong results #239
Description
Activity
- added this to the milestone
on Feb 19, 2016 This is a weird bug and possibly actually in rstan. If you add print statements after the end of the loop in generated quantities like
print(""); print("N[K] = ", N[K]); print("D[K] = ", D[K]); print("bcl of N[K] and D[K] = ", binomial_coefficient_log(N[K], D[K]));you get the wrong answer
N[K] = 100000 D[K] = 91116 bcl of N[K] and D[K] = 27676.8but if you call expose_stan_functions() on a file that just has
functions { real bcl(real x, real y) return binomial_coefficient_log(x,y); } model {}and then call
bcl(100000, 91116)you get the right answer, which is 29979.16. So, perhaps it has something to do with repeated calls to binomial_coefficient_log???
More weirdness. It looks as if the second argument is corrupted somehow
Chain 1, Iteration: 1 / 1 [100%] (Sampling) N[K] = 100000 D[K] = 91116 bcl of N[K] and D[K] = 27676.8 N[K] == 100000 is 1 D[K] == 91116 is 0but even if its last digit is off 27676.8 is nowhere close to the right answer.
Okay, the wrong answers occur when it takes the
elsebranch inbinomial_coefficient_logbut theelsebranch seems to give the right answer when you type it into R manually. And that does not explain why calling it viaexpose_stan_functions()yields the correct answer.I found it. It's integer division:
return n * log(N - n) + (N + 0.5) * log(N/(N-n)) + 1/(12*N) - n - 1/(12*(N-n)) - lgamma(n + 1.0);The problem is that
log(N/(N - n))is evaluatingN/(N - n)with integer division. I betexpose_stan_functions()instantiates the code withdoubleinstead ofintand it's fine.I haven't traced through this code, but we could use Stirling's approximation instead and it should be fine:
return (N + 0.5) * log(N) - (n + 0.5) * log(n) - (N - n + 0.5) * log(N - n) + NEG_LOG_SQRT_TWO_PI;Any thoughts? I'll implement a fix tomorrow.
Nice catch!
No preference as to approximation.
- Bob
On Mar 29, 2016, at 11:54 PM, Daniel Lee [email protected] wrote:
I found it. It's integer division:
return n * log(N - n) + (N + 0.5) * log(N/(N-n))
- 1/(12_N) - n - 1/(12_(N-n)) - lgamma(n + 1.0);
The problem is that log(N/(N - n)) is evaluating N/(N - n) with integer division. I bet expose_stan_functions() instantiates the code with double instead of int and it's fine.
I haven't traced through this code, but we could use Stirling's approximation instead and it should be fine:
return (N + 0.5) * log(N) - (n + 0.5) * log(n) - (N - n + 0.5) * log(N - n) + NEG_LOG_SQRT_TWO_PI;Any thoughts? I'll implement a fix tomorrow.
—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHubCan't we just replace n * log(N - n) + (N + 0.5) * log(N / (N - n)) ...
with
(n - 1) * log(N - n) + (N + 0.5) * log(N) ...
? But all those 1 / terms would to become 1.0 / .yes, we could. I don't know where that approximation comes from. I can follow the Stirling approx.
I'm guessing it doesn't matter so much since N >= 1000 or N-n >= 1000. Forward-mode is implemented incorrectly also, so I'll fix that too.
Originally submitted by @akodd as stan-dev/stan#1778
Summary:
Function
binomial_coefficient_logis not implemented correctly and causesbinomial_logto produce incorrect result.Description:
Probability mass function
binomial_logproduced incorrect results for some combinations of successes and trials.Reproducible Steps:
The following R code illustrates the issue:
Current Output:
Choosing only one example we get:
Note that R and Stan versions differ only in binomial coefficient.
I also verified the choose function on Wolfram for this one case http://www.wolframalpha.com/input/?i=log(1048+choose+10) to match R version.
A workaround for
binomial_logwould be a user defined function (_log suffix) that adds user supplied results of the choose function (can passed as data) to D_log(p)+(N-D)_log1m(p).Additional Information:
Current Version:
v2.9.0