Numerically Stable gamma_lcdf Gradients

[andrjohns]

Summary

This PR adapts the improved gamma_p gradients from #780 for use with gamma_lcdf, to improve stability of estimation with large inputs.

Tests

New mix/prob tests added for a range of inputs

Side Effects

N/A

Release notes

Improved numerical stability of gamma_lcdf gradients

Checklist

• Math issue change gamma_lcdf to just use gamma_p and the rev-mode derivatives #2763

• Copyright holder: Andrew Johnson

  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)
  • dependencies checks pass, (make test-math-dependencies)
  • 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

[rok-cesnovar]

Thanks Andrew!

[spinkney]

Would it work to just call gamma_p like in the below? Seems easier to understand.

Like

namespace stan {
namespace math {

template <typename T_y, typename T_shape, typename T_inv_scale>
return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
                                                    const T_shape& alpha,
                                                    const T_inv_scale& beta) {
  using T_y_ref = ref_type_t<T_y>;
  using T_alpha_ref = ref_type_t<T_shape>;
  using T_beta_ref = ref_type_t<T_inv_scale>;
  static const char* function = "gamma_lcdf";
  check_consistent_sizes(function, "Random variable", y, "Shape parameter",
                         alpha, "Inverse scale parameter", beta);
  T_y_ref y_ref = y;
  T_alpha_ref alpha_ref = alpha;
  T_beta_ref beta_ref = beta;
  check_positive_finite(function, "Shape parameter", alpha_ref);
  check_positive_finite(function, "Inverse scale parameter", beta_ref);
  check_nonnegative(function, "Random variable", y_ref);

  if (size_zero(y, alpha, beta)) {
    return 0;
  }

  auto const y_ = as_column_vector_or_scalar(y_ref);
  auto const alpha_ = as_column_vector_or_scalar(alpha_ref);
  auto const beta_ = as_column_vector_or_scalar(beta_ref);
  
  return log(gamma_p(alpha, beta_ * y_));
}

}  // namespace math
}  // namespace stan
#endif

[andrjohns]

Yep, but then the gradients for the outer log() and the b * y get auto-diffed through, so the sampling won't be as efficient

[spinkney]

I think only one thing. You can remove lines 66 - 77 as those were only needed in the previous version.
These lines

  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
      gamma_vec(math::size(alpha));
  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
      digamma_vec(math::size(alpha));
  if (!is_constant_all<T_shape>::value) {
    for (size_t i = 0; i < stan::math::size(alpha); i++) {
      const T_partials_return alpha_dbl = alpha_vec.val(i);
      gamma_vec[i] = tgamma(alpha_dbl);
      digamma_vec[i] = digamma(alpha_dbl);
    }
  }

Otherwise, looks great!

[andrjohns]

Good catch, thanks!