Inverse trigamma function in Stan

I’m wondering if someone in the Stan community (or the Stan developers themselves) has worked on a function to compute the inverse of the trigamma function. This is a particularly important function to use when you are working with a multivariate log-gamma distribution and have re-parameterized the density in terms of the mean and covariance. The reference where this is introduced is here:

Bradley, J.R., Holan, S.H. and Wikle, C.K., 2018. Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion). Bayesian Analysis, 13 (1), pp.253-310.

On page 258, you see that the moments of the multivariate log-gamma are a function of the digamma and trigamma functions, and that when you re-parameterize in terms of the mean and covariance (page 289), the inverse trigamma function shows up.

It seems like there exist implementations of this function outside of Stan (e.g., http://web.mit.edu/~r/current/arch/i386_linux26/lib/R/library/limma/html/trigammainverse.html). Thanks for providing some insight on this issue!

Hi @berchuck, welcome to the Stan forums. I know we have an implementation of trigamma but I don’t think we have the inverse. I don’t know if anyone has worked on it. Are there known numerical issues when computing the trigamma inverse? Often that’s why a function is missing from Stan even though its inverse is available.

If you can find an implementation in C++ then that would be the easiest way to use it in Stan, either by calling the external C++ code from the Stan program or by reimplementing the C++ in Stan code (if possible). We already depend on the boost c++ library so I took a quick look there and I didn’t see an implementation (I saw trigamma but not the inverse, though I may have missed it). I haven’t looked at the source code for the R version that you linked to, but you could also try copying that implementation in Stan code if it’s possible.