Quick fix for notable problems with neg_binomial_2 and neg_binomial_2_log

[martinmodrak]

Summary

Minimal changes to avoid some (but not all) pathological behaviour of the neg_binomial_2 variants.

neg_binomial_2 no longer overrides the logp accumulator when phi > 1e5. This cutoff for phi is too low, but I ended up not increasing it as larger cutoff breaks the current finite diffs tests (purely because the test is wrong, the returned derivatives are good, finite diffs are waaaaay off).

neg_binomial_2_log now delegates to Poisson for phi > 1e5 to avoid returning positive log. prob.

In both cases, I am aware that the computation is not stable in some cases, but those fixes should improve the situation for current release.

Tests

• Test for neg_binomial_2_lpmf where a vector of values, some below and some above the phi cutoff are passed.
• Test for neg_binomial_2_log_lpmf comparing against precomputed values from Mathematica (cannot be computed for values of phi above the 1e10 cutoff).
• Comparing the results of neg_binomial_2_lpmf and neg_binomial_2_log_lpmf across a broad range of values.
• Comparing the derivatives of neg_binomial_2_lpmf against a complex-step derivative across a range of values.

Side Effects

When propto=TRUE and some phi values trigger the Poisson code path, the normalizing constant is not properly handled (for neg_binomial_2_lpmf, this was the case before this PR, but now this will also affect neg_binomial_2_log_lpmf)

Checklist

• Math issue neg_binomial_2_log_lpmf is not stable for large phi #1495 and partially The way neg_binomial_2_lpmf delegates to Poisson is broken #1496

• Copyright holder: Institute of Microbiology of the Czech Academy of Sciences

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[martinmodrak]

@syclik Those are IMHO the only fixes worth pushing to the release. It's getting late here, so will likely respond to any inquiries tomorrow morning.

[bbbales2]

This all looks good. I added a couple questions on one of the tests. I agree looking in R it seems like the threshold to turn this into a Poisson should be much higher. Maybe we can work out a better condition in a later pull.

I think it makes sense to add simple autodiff tests since it's easy and makes it clear that the values and gradients are consistent with each other. Something super rough like:

#include <test/unit/math/test_ad.hpp>
#include <limits>

TEST(mathMixScalFun, neg_binomial_2_log_lpmf_derivatives) {
  auto f1 = [](const auto& eta, const auto& phi) { return stan::math::neg_binomial_2_log_lpmf(0, eta, phi); };
  auto f2 = [](const auto& eta, const auto& phi) { return stan::math::neg_binomial_2_log_lpmf(6, eta, phi); };

  stan::test::expect_ad(f1, -1.5, 4.1);
  stan::test::expect_ad(f1, 2.0, 1.1);
  stan::test::expect_ad(f2, -1.5, 4.1);
  stan::test::expect_ad(f2, 2.0, 1.1);
}

and then the same for neg_binomial_2_lpmf.

[bbbales2]

What do the complex step tests cover that the others don't?

These fancy tests are for neg_binomial_2_log_lpmf. Should they cover neg_binomial_2_lpmf as well?

[martinmodrak]

What do the complex step tests cover that the others don't?

They are able to test gradients numerically to (much) higher precision than what is possible with finite diffs. They are also able to cover large phi values where the precomputed test is not used, since Mathematica (at least the free cloud version) gives up on computing exact answers for large phi.

These fancy tests are for neg_binomial_2_log_lpmf. Should they cover neg_binomial_2_lpmf as well?

Well, they should :-) And I do have them written in #1497 . The thing is I am aiming for really a minimal PR and the tests fail (unless weakened noticeably), because there are all sorts of minor issues with the numerics that require more time to be handled. I don't think there is any particularly defensible line on what tests to include (and hence what issues to address). I went with "avoid obvious bugs (positive log prob., overwriting logp)".

You may disagree on the specific choice or with the whole concept of the quick fix, obviously. I am myself not 100% happy with leaving so many holes, and maybe waiting for a full fix is a better option in some sense.

[martinmodrak]

Thanks @bbbales2 for looking into it.

I think it makes sense to add simple autodiff tests since it's easy and makes it clear that the values and gradients are consistent with each other.

Yes, I added those, but between the distribution tests and the precomputed tests from Mathematica, I am not sure this adds much value.

I agree looking in R it seems like the threshold to turn this into a Poisson should be much higher. Maybe we can work out a better condition in a later pull.

I think I already did that in #1497 (you may witness me slowly starting to grasp this over time in the comments :-) ). The answer seems to be: "Delegating to Poisson is finicky. If you improve the numerics, you can avoid the Poisson branch completely and circumvent of all sorts of issues (ensuring continuity around the cutoff, making sure propto=true behaves consistently across the branches, ...)".

[martinmodrak]

The build fails on "Windows headers & unit" with

MathSetup— Restore files previously stashed1m 31s
java.nio.channels.ClosedChannelException

And nothing more... is this an issue with Jenkins or is there somewhere to investigate what that means?

[martinmodrak]

Also, I had to move the cutoff for phi to delegate to Poisson back to 1e5 for both versions, as the 1e10 triggered a lot of false-positive test failures that I am not willing to address in this PR (but those are handled in #1497)

[bbbales2]

I am not sure this adds much value

Just to be clear (since I agree it's kindof annoying for me to ask for these), it gives me confidence that all the basic mechanics are right. Things like double, double vs var, double vs double, var vs fvar, double work.

Also if I spot check the prim numerics once in R or something, then I can kinda have faith that things are generalizing, and that gives me confidence that your high precision tests are testing the right thing (and aren't just high precision testing something else).

The build fails on "Windows headers & unit" with

I restarted the tests. We can merge when the tests pass.

[martinmodrak]

it gives me confidence that all the basic mechanics are right. Things like double, double vs var, double vs double, var vs fvar, double work.

I agree this is important, I just thought this is what the "Distribution tests" (the tests under test/prob) are for. Should I be generally worried that distribution tests don't catch this type of issues and test this stuff explicitly? Also sorry, if I sounded annoyed.

I restarted the tests.

Thanks

[bbbales2]

I just thought this is what the "Distribution tests" (the tests under test/prob) are for. Should I be generally worried that distribution tests don't catch this type of issues and test this stuff explicitly?

Good point lol.

Also sorry, if I sounded annoyed.

Nah nah, you're good.

[bbbales2]

Oops, I shoulda caught this. This should be in mix, but given how long the tests have taken merge this anyway and we can just move this to mix with the next pull.

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	4.89	4.86	1.01	0.66% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.97	-3.19% slower
eight_schools/eight_schools.stan	0.09	0.09	0.99	-0.99% slower
gp_regr/gp_regr.stan	0.22	0.23	1.0	-0.2% slower
irt_2pl/irt_2pl.stan	6.1	6.07	1.01	0.59% faster
performance.compilation	87.64	87.17	1.01	0.54% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	7.31	7.3	1.0	0.1% faster
pkpd/one_comp_mm_elim_abs.stan	20.4	20.23	1.01	0.81% faster
sir/sir.stan	93.4	93.38	1.0	0.03% faster
gp_regr/gen_gp_data.stan	0.04	0.04	1.01	0.54% faster
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.97	2.96	1.0	0.19% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.31	0.31	1.01	0.9% faster
arK/arK.stan	1.76	1.74	1.01	0.74% faster
arma/arma.stan	0.8	0.8	0.99	-0.57% slower
garch/garch.stan	0.63	0.63	0.99	-0.88% slower
Mean result: 0.999620446184

---

Jenkins Console Log
Blue Ocean
Commit hash: c6febc4

---

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

[mcol]

I think this can be merged, right? I'm pretty sure @syclik was fine with it if it had tests and approval.