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., 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.