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?