The way neg_binomial_2_lpmf delegates to Poisson is broken

[martinmodrak]

Description

The neg_binomial_2_lpmf delegates to poisson_lpmf for large values of the phi parameter. Currently the cutoff is set at phi > 1e5. This introduces discontinuities into the function - I am nut sure whether this is because the cutoff is too optimistic or too pessimistic.

Example

neg_binomial_2_lpmf(181, 13, 1e5 - 1e-8) == -312.05065787932836
neg_binomial_2_lpmf(181, 13, 1e5 + 1e-8) == -312.19077641210203

neg_binomial_2_lpmf(121358, 13, 1e5 - 1e-8) == -933542.7362432084
neg_binomial_2_lpmf(121358, 13, 1e5 + 1e-8) == -988062.09132748179

I've written a failing test for the cutoff at: develop...martinmodrak:bugfix/1496-poisson-phi-cutoff

Expected Output

Values at cutoff - 1e-8 and cutoff + 1e-8 are almost constant, i.e. changes in phi no longer influence the lpmf value and delegating to Poisson is safe.

Current Version:

v3.0.0

[martinmodrak]

It gets even more funny - the cutoff is implemented as

    // if phi is large we probably overflow, defer to Poisson:
    if (phi__[i] > 1e5) {
      logp = poisson_lpmf(n_vec[i], mu__[i]);
    }

Note the =!!! So when a single element of a vector provokes cutoff, it overwrites the logp accumulator. Also, the way poisson_lpmf is called means that the propto information cannot be conserved and thus the final logp will wildly change around the cutoff under certain propto combinations. On the same branch I wrote a failing test for this behavior as well.

[bob-carpenter]

It gets even more funny

That's just a flat-out bug. The conditional should have been at the top of the loop and have the rest of the body of the loop in the alternative branch.

Thanks for tracking this down!

[bob-carpenter]

On Dec 9, 2019, at 12:42 PM, Martin Modrák ***@***.***> wrote: It gets even more funny - the cutoff is implemented as // if phi is large we probably overflow, defer to Poisson: if (phi__[i] > 1e5) { logp = poisson_lpmf(n_vec[i], mu__[i]); } Note the =!!!

That's just a flat-out bug.

…

So when a single element of a vector provokes cutoff, it overwrites the logp accumulator. Also, the way poisson_lpmf is called means that the propto information cannot be conserved and thus the final logp will wildly change around the cutoff under certain propto combinations. On the same branch I wrote a failing test for this behavior as well. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[martinmodrak]

One thing I discovered when working on this is that nb(n, mu, phi) = poisson(n, mu) + O(1/phi), so maybe delegating to poisson forces us to have too large a threshold. We could use further terms in the Taylor expansion at phi -> inf to obtain nb(n, mu, phi) = poisson(n, mu) + (mu^2 -n -2*mu*n +n ^2)/phi + O(1/phi^2).

This would (probably) let us lower the threshold and use the simpler derivatives for smaller phi... Currently, I need the threshold to be 1e10 (same as similar threshold in neg_binomial_lpmf) to make the tests pass.

Does anybody have any opinion on this? I am slightly reluctant to add this as the current test coverage of the neg. binomial seems not really great...

[martinmodrak]

Update: tested with the additional term, doesn't allow me to reduce the cutoff.

[martinmodrak]

So I probably spent more time than I have with this, will probably be able to continue only after Christmas, maybe somebody else could look into it.

It appears that the neg_binomial_2_lpmf is riddled with small numerical issues. I wrote more tests for the branch (as seen in the PR: #1497). There are some numerical differences between the values and derivatives as we compute them and what I get from Mathematica (the biggest relative errors are on the order of 10^-5, not sure if this can be ignored).

I also failed to get the finite differences test to pass, though I am not sure what could have triggered its failures, as it fails for phi values before cutoff where I didn't make big changes.

[bob-carpenter]

1e-5 relative errors of autodiff relative to someone's else's numerical eval or to finite diffs? The latter should be fine.

…

On Dec 10, 2019, at 10:01 AM, Martin Modrák ***@***.***> wrote: So I probably spent more time than I have with this, will probably be able to continue only after Christmas, maybe somebody else could look into it. It appears that the neg_binomial_2_lpmf is riddled with small numerical issues. I wrote more tests for the branch (as seen in the PR: #1497). There are some numerical differences between the values and derivatives as we compute them and what I get from Mathematica (the biggest relative errors are on the order of 10^-5, not sure if this can be ignored). I also failed to get the finite differences test to pass, though I am not sure what could have triggered its failures, as it fails for phi values before cutoff where I didn't make big changes. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[martinmodrak]

Moving the discussion from here to the PR #1497 - as it IMHO relates more to the code than to the issue.