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!