gamma_p and gamma_q gradient problems

[bbbales2]

Description

gamma_p doesn't work great with vars and gamma_q doesn't work great with fvar<double>

Here is test code:

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

TEST(mathMixScalFun, gamma_q) {
  using stan::math::var;
  using stan::math::fvar;
  using stan::math::gamma_q;
  using stan::math::gamma_p;
  using stan::math::gamma_lccdf;
  using stan::math::log;
  using stan::math::log1m;

  var yv = 12;
  var alphav = 3;
  var betav = 3;

  fvar<double> yf(12.0, 1.0);
  fvar<double> alphaf(3.0, 1.0);
  fvar<double> betaf(3.0, 1.0);

  var lp1 = log1m(gamma_p(alphav, betav * yv));
  fvar<double> lp2 = log(gamma_q(alphaf, betaf * yf));

  lp1.grad();

  std::cout << "Computed: " << std::endl;
  std::cout << "lp1.val(): " << lp1.val() << std::endl;
  std::cout << "lp2.val(): " << lp2.val() << std::endl;
  std::cout << "y.adj(): " << yv.adj() << std::endl;
  std::cout << "alpha.adj(): " << alphav.adj() << std::endl;
  std::cout << "beta.adj(): " << betav.adj() << std::endl;
  std::cout << "lp2.d_: " << lp2.d_ << std::endl;
}

TEST(mathMixScalFun, gamma_q_ref) {
  using stan::math::var;
  using stan::math::fvar;
  using stan::math::gamma_q;
  using stan::math::gamma_p;
  using stan::math::gamma_lccdf;
  using stan::math::log;
  using stan::math::log1m;

  var yv = 12;
  var alphav = 3;
  var betav = 3;

  fvar<double> yf(12.0, 1.0);
  fvar<double> alphaf(3.0, 1.0);
  fvar<double> betaf(3.0, 1.0);

  var lp1 = gamma_lccdf(yv, alphav, betav);
  fvar<double> lp2 = gamma_lccdf(yf, alphaf, betaf);

  lp1.grad();

  std::cout << "Reference: " << std::endl;
  std::cout << "lp1.val(): " << lp1.val() << std::endl;
  std::cout << "lp2.val(): " << lp2.val() << std::endl;
  std::cout << "y.adj(): " << yv.adj() << std::endl;
  std::cout << "alpha.adj(): " << alphav.adj() << std::endl;
  std::cout << "beta.adj(): " << betav.adj() << std::endl;
  std::cout << "lp2.d_: " << lp2.d_ << std::endl;
}

The output for gamma_p and gamma_q are (computing the same function, lp1 is with gamma_p and lp2 is with gamma_q):

Computed: 
lp1.val(): -29.4707
lp2.val(): -29.4706
y.adj(): 0
alpha.adj(): 0
beta.adj(): 0
lp2.d_: -6.70499e+12

And the reference values are here:

Reference: 
lp1.val(): -29.4706
lp2.val(): -29.4706
y.adj(): -2.83796
alpha.adj(): 2.68924
beta.adj(): -11.3518
lp2.d_: -11.5005

(edit: updated output)

Current Version:

v3.3.0

[syclik]

Thanks for posting code. I was trying to verify it and I got something different:

Computed: 
lp1.val(): -29.4707
lp2.val(): -29.4706
y.adj(): 0
alpha.adj(): 0
beta.adj(): 0
lp2.d_: -6.70499e+12

Reference: 
lp1.val(): -29.4706
lp2.val(): -29.4706
y.adj(): -2.83796
alpha.adj(): 2.68924
beta.adj(): -11.3518
lp2.d_: -11.5005

Is there some reason my output is different than yours (for even the .val())?

[bbbales2]

@syclik aaah, yeah you're right I updated the numbers in the reference code without regenerating the output. I updated the post.

[bbbales2]

I think the issue is that the fwd gamma_q is an interative algorithm with what might be a coarse tolerance: https://github.com/stan-dev/math/blob/develop/stan/math/fwd/fun/gamma_q.hpp#L32

And for rev gamma_p the gradient gets thresholded to zero: https://github.com/stan-dev/math/blob/develop/stan/math/rev/fun/gamma_p.hpp#L39

[wds15]

Should we try what boost offers here:

https://www.boost.org/doc/libs/1_72_0/libs/math/doc/html/math_toolkit/sf_gamma/gamma_derivatives.html

Maybe this stuff is new? I don't know the history, but this is worth considering.

[bbbales2]

Maybe. It looks like that gives one of the two gradients in each case. I'm not sure which is causing the problems. I just made the issue so it's known (didn't really plan on fixing it)