[On hold] New formulas for derivatives of neg_binomial_lpmf after alpha cutoff

[martinmodrak]

Summary

Currently on-hold. If PR #1497 gets merged, neg_binomial_lpmf should probably be rewritten using the same approach

The formulas for derivatives of neg_binomial_lpmf after delegating to poisson (for alpha > 1e10) are replaced with Taylor expansions of the exact derivatives for alpha -> Inf obtained in Mathematica via:

     nb[n_,alpha_,beta_]:= LogGamma[n + alpha] - LogGamma[n + 1] - 
        LogGamma[alpha ] + alpha * Log[beta/ (1 + beta)] - n * Log[1 + beta];
      nbdalpha[n_,alpha_,beta_]= D[nb[n, alpha, beta],alpha];
      nbdbeta[n_,alpha_,beta_]= D[nb[n, alpha, beta],beta];
      Series[nbdalpha[n, alpha, beta],{alpha, Infinity, 1}]
      Series[nbdbeta[n, alpha, beta],{alpha, Infinity, 0}]

(the lowest order of the series that let's me pass the tests was chosen)

Tests

Tests are adapted from #1497 where new tests for neg_binomial_2_lpmf were written. Namely

• Created a set of test values and derivatives in Mathematica (Wolfram Cloud) and test against those
• Testing the derivatives against complex step derivatives. This is important for large alpha, as Mathematica refuses to compute exact derivatives in this case.
• Over a grid of beta and n values test whether the change in the derivatives of neg_binomial_lpmf at alpha_cutoff +/- 1e-8 is negligible.
• Test whether there is no drastic change in neg_binomial_lpmf<propto=true, int, var, double> between alpha just below and just above the cutoff.

Side Effects

Not any.

Checklist

• Math issue The derivatives of neg_binomial_lpmf are incorrect for large alpha #1578

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

  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

[stan-buildbot]

(stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan, 1.02)
(stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan, 0.98)
(stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan, 1.0)
(stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan, 1.04)
(stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan, 1.02)
(stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan, 0.99)
(stat_comp_benchmarks/benchmarks/arK/arK.stan, 0.99)
(performance.compilation, 1.02)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan, 0.99)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan, 1.0)
(stat_comp_benchmarks/benchmarks/sir/sir.stan, 0.97)
(stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan, 0.92)
(stat_comp_benchmarks/benchmarks/garch/garch.stan, 1.0)
(stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan, 0.94)
(stat_comp_benchmarks/benchmarks/arma/arma.stan, 1.0)
Result: 0.99186497357
Commit hash: 095974e

[martinmodrak]

@bob-carpenter May I ask for review on this? The tests here are similar to those written for #1497 (although the problem addressed is different).

[bob-carpenter]

Looks good. My only requests were around doc, which shouldn't be first person and should be written so that someone who knows C++ but not Mathematica can understand it if you want to leave it in the code.

Also, have you checked that the IP is OK for us to include Mathematica code in comments and use the output of their approximations? I can't see why it wouldn't be, but it'd be good to check before we go too far down the road of using a proprietary software package.

[bob-carpenter]

Comments shouldn't be first person. I'd suggest just eliminating line 19.

Line 20 should say what the value is for---in this case, the threshold beyond which a Poisson approximation is used. It's also OK to say why it's 1e10, but I couldn't follow the logic of why it's set to 1e10.

[syclik]

@bob-carpenter, would it be ok if L20 said something like:
// Set to 1e10 based on empirical evaluation. At 1e9, the approximation isn't close enough

[bob-carpenter]

I don't think we should assume the readers of comments understand mathematica. So this should be put in Stan notation or just eliminated. I couldn't understand what LogGamma is or Log or D. If these comments stay, the functions used need to be explained.

[stan-buildbot]

(stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan, 1.0)
(stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan, 1.01)
(stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan, 1.01)
(stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan, 0.99)
(stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan, 0.99)
(stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan, 1.0)
(stat_comp_benchmarks/benchmarks/arK/arK.stan, 1.0)
(performance.compilation, 1.0)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan, 1.0)
(stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan, 1.0)
(stat_comp_benchmarks/benchmarks/sir/sir.stan, 1.0)
(stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan, 1.01)
(stat_comp_benchmarks/benchmarks/garch/garch.stan, 1.0)
(stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan, 1.01)
(stat_comp_benchmarks/benchmarks/arma/arma.stan, 1.0)
Result: 1.00075771834
Commit hash: ce88c83

[bob-carpenter]

@bob-carpenter, would it be ok if L20 said something like: // Set to 1e10 based on empirical evaluation. At 1e9, the approximation isn't close enough

Yes, that's clear.

[martinmodrak]

Thx @syclik and @bob-carpenter for the notes. I should've marked this PR as "on hold" because if #1497 turns out to be on the right track, I would probably rewrite neg_binomial_lpmf in the same spirit (eliminating the Poisson branch completely). So sorry for wasting your time here.

[syclik]

@martinmodrak, do you want to close this PR? #1497 was merged and from reading through the PR, it seems like you are favoring a different implementation.

[martinmodrak]

Thanks for reminidng me, yes this should be closed.