Review for math bugfix

Any chance this can get reviewed in time for 2.15? I suggested Thel on the github menu but no idea if that sends a notification or who has time to look at it. At this point it all passes tests, has tests, and avoided trying to improve doc (but the few new functions have doc).

https://github.com/stan-dev/math/pull/515

Hi!

I could offer to do a pass over it on Thursday. I have to say not being an expert in these functions, but with some luck I find one of my models which as I recall uses these functions.

From a very quick look over it I noticed that the documentation quoted 1E4 maximum steps while the function used 1E5 as default.

This looks like massive work and it should go into 2.15 if we can.

Sebastian

I think Mitzi was planning to tag ASAP. I don’t know if she’s on Discourse, but I can let her know since she’s sitting right here. I can look at this if necessary, but I’d really rather have someone else review this.

I would encourage us to wait until this goes in to release 2.15 because we currently have bugs (fixed in the PR) that will cause Stan models to silently hang indefinitely and at least one of the functions can terminate early without throwing and therefore return a wrong answer. Not sure I have a model that demonstrates this second problem but the hanging def. happens if you code a generalized gamma in Stan even with decent priors. I’m willing to fix problems found in review same day or next day this week.

planning to do a release soon, but happy to wait for someone to do a proper review of this so that it gets into the release.

by “planning to do a release soon” I mean this week. my understanding is that we need to get the release out ASAP because of R. What are the consequences of releasing 2.15.0 without this feature?

  1. Models that hang, this has come up on the users list, this comes up in models that I have implemented using the incomplete gamma function, not sure if there are other user-facing functions that trigger it but this affects user-facing functions (e.g.-gamma_cdf).

  2. Models will potentially produce wrong answers if they on grad_2F1. This includes student’s t cdf and some of the beta-binomial functions.

  3. There are also a couple of places where the sign of calculations were not correctly tracked but I don’t have an example of a user-facing function where it matters.

Here’s the function with the arbitrary exit regardless of convergence: https://github.com/stan-dev/math/blob/develop/stan/math/prim/scal/fun/grad_2F1.hpp

Here’s where grad_2F1 gets used, most relevant bit at the end (sorry). Looks like student’s t cdf, negative binomial a couple of functions.

krzysztof@moze:~/packages/math$ grep 'grad_2F1' -R ./stan/math
./stan/math/rev/scal/fun/grad_inc_beta.hpp:#include <stan/math/prim/scal/fun/grad_2F1.hpp>
./stan/math/rev/scal/fun/grad_inc_beta.hpp:        grad_2F1(dF1, dF2, a + b, var(1.0), a + 1, z);
./stan/math/fwd/scal/fun/grad_inc_beta.hpp:#include <stan/math/prim/scal/fun/grad_2F1.hpp>
./stan/math/fwd/scal/fun/grad_inc_beta.hpp:        grad_2F1(dF1, dF2, a + b,
./stan/math/prim/scal/fun/grad_inc_beta.hpp:#include <stan/math/prim/scal/fun/grad_2F1.hpp>
./stan/math/prim/scal/fun/grad_inc_beta.hpp:      if (C) grad_2F1(dF1, dF2, a + b, 1.0, a + 1, z);
./stan/math/prim/scal/fun/grad_2F1.hpp:    void grad_2F1(T& g_a1, T& g_b1, const T& a1, const T& a2, const T& b1,
./stan/math/prim/scal/fun/grad_2F1.hpp:      check_2F1_converges("grad_2F1", a1, a2, b1, z);
./stan/math/prim/scal/fun/grad_2F1.hpp:          domain_error("grad_2F1", "k (internal counter)", max_steps,
./stan/math/prim/scal.hpp:#include <stan/math/prim/scal/fun/grad_2F1.hpp>
krzysztof@moze:~/packages/math$ grep 'grad_inc_beta' -R ./stan/math
./stan/math/rev/scal/fun/grad_inc_beta.hpp:    inline void grad_inc_beta(var& g1, var& g2,
./stan/math/rev/scal.hpp:#include <stan/math/rev/scal/fun/grad_inc_beta.hpp>
./stan/math/fwd/scal/fun/grad_inc_beta.hpp:    void grad_inc_beta(fvar<T>& g1,
./stan/math/fwd/scal/fun/inc_beta.hpp:#include <stan/math/fwd/scal/fun/grad_inc_beta.hpp>
./stan/math/fwd/scal.hpp:#include <stan/math/fwd/scal/fun/grad_inc_beta.hpp>
./stan/math/prim/scal/fun/grad_reg_inc_beta.hpp:#include <stan/math/prim/scal/fun/grad_inc_beta.hpp>
./stan/math/prim/scal/fun/grad_reg_inc_beta.hpp:      grad_inc_beta(dBda, dBdb, a, b, z);
./stan/math/prim/scal/fun/grad_inc_beta.hpp:    inline void grad_inc_beta(double& g1, double& g2,
./stan/math/prim/scal.hpp:#include <stan/math/prim/scal/fun/grad_inc_beta.hpp>
krzysztof@moze:~/packages/math$ grep 'grad_reg_inc_beta' -R ./stan/math
./stan/math/rev/scal/fun/inc_beta.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/rev/scal/fun/inc_beta.hpp:          grad_reg_inc_beta(d_a, d_b, avi_->val_, bvi_->val_,
./stan/math/fwd/scal/fun/inc_beta.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/fwd/scal/fun/inc_beta.hpp:      grad_reg_inc_beta(d_a, d_b, a.val_, b.val_, x.val_,
./stan/math/prim/scal/fun/grad_reg_inc_beta.hpp:    void grad_reg_inc_beta(T& g1, T& g2, const T& a, const T& b,
./stan/math/prim/scal/prob/beta_lccdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/beta_lccdf.hpp:          grad_reg_inc_beta(g1, g2, alpha_dbl, beta_dbl, y_dbl,
./stan/math/prim/scal/prob/neg_binomial_lccdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_lccdf.hpp:          grad_reg_inc_beta(g1, g2, alpha_dbl,
./stan/math/prim/scal/prob/beta_lcdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/beta_lcdf.hpp:          grad_reg_inc_beta(g1, g2, alpha_dbl, beta_dbl, y_dbl,
./stan/math/prim/scal/prob/neg_binomial_2_log_rng.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_lcdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_lcdf.hpp:          grad_reg_inc_beta(g1, g2, alpha_dbl,
./stan/math/prim/scal/prob/neg_binomial_2_lccdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/student_t_lcdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/student_t_lcdf.hpp:            grad_reg_inc_beta(g1, g2, 0.5 * nu_dbl,
./stan/math/prim/scal/prob/student_t_lcdf.hpp:            grad_reg_inc_beta(g1, g2, (T_partials_return)0.5,
./stan/math/prim/scal/prob/beta_lpdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/beta_rng.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_2_log_lpmf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/student_t_lccdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/student_t_lccdf.hpp:            grad_reg_inc_beta(g1, g2, 0.5 * nu_dbl,
./stan/math/prim/scal/prob/student_t_lccdf.hpp:            grad_reg_inc_beta(g1, g2, (T_partials_return)0.5,
./stan/math/prim/scal/prob/student_t_cdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/student_t_cdf.hpp:            grad_reg_inc_beta(g1, g2, 0.5 * nu_dbl,
./stan/math/prim/scal/prob/student_t_cdf.hpp:            grad_reg_inc_beta(g1, g2, (T_partials_return)0.5,
./stan/math/prim/scal/prob/student_t_lpdf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_2_lpmf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_2_rng.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/student_t_rng.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_rng.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal/prob/neg_binomial_lpmf.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>
./stan/math/prim/scal.hpp:#include <stan/math/prim/scal/fun/grad_reg_inc_beta.hpp>

I’m planning to review this tonight.

Awesome, I’ll fix whatever you find tomorrow.