stan-dev
Add hypergeometric function of a matrix argument
Overview
Provide an implementation for the hypergeometric function of a matrix argument, a generalization of the hypergeometric series.
Description
The function is commonplace in the densities of matrix variates, e.g. several distributions in the Matrix Variate Distributions text by Gupta and Nagar. The function is also described on Wikipedia.
Expected Output
Implement both the function and its logarithm as functions that can be called within a Stan program:
- [ ]
real hypergeometric_matrix() - [ ]
real log_hypergeometric_matrix()
Additional Resources
- Efficient algorithm for the function: Koev, P., & Edelman, A. (2006). The efficient evaluation of the hypergeometric function of a matrix argument. Mathematics of Computation, 75(254), 833–846. https://doi.org/10.1090/S0025-5718-06-01824-2
- With C code for MATLAB: https://math.mit.edu/~plamen/software/mhgref.html
- R implementation using C++: https://cran.r-project.org/package=HypergeoMat
- GitHub source code: https://github.com/stla/HypergeoMat/blob/master/src/HypergeoMat.cpp
Current Version:
v4.4.0
Hello, I figure there's a lot the team has on its plate, but thought to add this in case the team gets to it someday. I hope I provided enough information in the original post. The function has lots of applications in multivariate stats, random matrix theory. My particular application is covariance structure analysis (psychometrics). Thank you.
The difficulty with implementing hypergeometric series and their extensions comes from the need to also be able to implement their derivatives/gradients, which as I've experienced in this PR can get pretty gnarly.
However, you've opened this issue at the perfect time! Once the above hypergeometric_pFq PR is merged, we can use it to evaluate matrix arguments through one of the methods discussed over in this Mathematica post: https://mathematica.stackexchange.com/questions/58298/hypergeometric-function-with-a-matrix-argument
Although I would note that the stability/speed of these functions will likely be somewhat fragile - since an infinite series is needed to calculate the value, and then a double infinite series needed for the gradients
