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:
 "Error in sampler$call_sampler(args_list[[i]]) : "
 " 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?