Error with gamma_q but not (1-gamma_p)

Algorithms
1

Hello everyone,

I was trying to fit, using synthetic data, a complicated model (whose code I’ll post soon) in which some variables are related by the regularized upper incomplete gamma function Q. The Stan function library contains this function, gamma_q.

I encounter the following error:

[1] "Error in sampler$call_sampler(args_list[[i]]) : "
[2] " grad_reg_inc_gamma: is not converging "
error occurred during calling the sampler; sampling not done

However, when I substitute in (1 - gamma_p(…)) where i used gamma_q, then the model successfully** starts and completes sampling. (**Successful in running, but recovers samples that are clearly biased…)

Mathematically, the regularized lower incomplete gamma function P is related to Q by the equation P + Q = 1, and there should be no difference when substituting Q and 1-P.

I conclude that there is some sort of algorithmic difference in the way grad_reg_inc_gamma is calculated for gamma_q, compared to the Stan math function for gamma_p’s gradient.

Am I using Stan’s gamma_q and gamma_p functions incorrectly, or did I run into some numerical stability bug?

2

Very likely an issue with Stan’s numerics. Can you make an issue in the Math library: Sign in to GitHub · GitHub ? If possible, can you use print to print out arguments that cause this calculation to blow up? Is there a way we could know which values of arguments to switch from one implementation to the other?

If you’re a C++er and you want to work on a fix that would be appreciated as well :D.

3

Sounds like something Philip could easily sort out. But I can’t find his handle on Discourse.

4

I forwarded him a message!