@fomotis
The first step is to specify the c++ function. The c++ for the gamma_q_inv function and its gradients is:
#include <boost/math/special_functions/hypergeometric_pFq.hpp>
#include <stan/math/prim/scal.hpp>
#include <stan/math/rev/core.hpp>
// Function to handle Stan parameters
inline stan::math::var gamma_q_inv_stan(const stan::math::var& a,
const stan::math::var& q, std::ostream* pstream__) {
// Extract values (not derivatives of inputs)
double a_val = a.val();
double q_val = q.val();
// Calculate inverse of the incomplete gamma function
double gamma_q_inv_val = boost::math::gamma_q_inv(a_val, q_val);
// Pre-compute some quantities for the gradients
double exp_gamma_q_inv_val = std::exp(gamma_q_inv_val);
double pow_gamma_q_inv_val = std::pow(gamma_q_inv_val,(1 - a_val));
double gamma_a_val = stan::math::tgamma(a_val);
// Calculate derivative wrt to 'a' parameter
double a_d = exp_gamma_q_inv_val
* pow_gamma_q_inv_val
* (stan::math::square(gamma_a_val)
* std::pow(gamma_q_inv_val,a_val)
// Boost doesn't have a normalised Hypergeometric pFq, implement our own
* (boost::math::hypergeometric_pFq({a_val,a_val},{a_val+1,a_val+1},-gamma_q_inv_val)
/ stan::math::square(stan::math::tgamma(a_val+1)))
+ (q_val - 1)
* gamma_a_val
* std::log(gamma_q_inv_val)
+ boost::math::polygamma(0, a_val)
* (gamma_a_val - boost::math::tgamma(a_val, gamma_q_inv_val)));
// Calculate derivative wrt 'q' parameter
double q_d = gamma_a_val * -exp_gamma_q_inv_val * pow_gamma_q_inv_val;
// Return new parameter with values of the function and gradients
return stan::math::var(new stan::math::precomp_vv_vari(gamma_q_inv_val, a.vi_, q.vi_, a_d, q_d));
}
I’ve checked the return values and gradients for a couple of inputs against those returned by Wolfram Alpha.
Save this in the same folder as your stan code, with the name gamma_q_inv.hpp. The name of the file isn’t important, but just so you can replicate my results.
Next, you need to include a functions definition in your Stan model that has the required signature, for example:
functions {
real gamma_q_inv_stan(real a, real q);
}
data {
real y_mean;
}
parameters {
real<lower=0> y;
}
model {
y ~ normal(y_mean, 1);
}
Next in R, navigate to the same directory as your Stan model and gamma_q_inv.hpp file. It’s not a requirement that you be in the same directory, you just have to adjust your file paths accordingly.
First compile the Stan model, including the external .hpp file. using the stan_model function:
mod = stan_model("model_gamma.stan", allow_undefined = TRUE,
includes = paste0('\n#include "',file.path(getwd(), 'gamma_q_inv.hpp'), '"\n'))
Now you can sample your model using the sampling function:
samp = sampling(mod, data=list(y_mean=0))
It’s a bit convoluted, but I’ll add updating/revamping the instructions to my to-do list