Gamma distribution quantile and inverse quantile function

General
, , , ,
1

My model involves solving an equation of the F(d, b e^{k}) = \frac{q}{c - 1} for b, where F(.) is the gamma CDF.
Obviously,

b = \frac{F^{-1}(\frac{q}{c - 1}, 1, d)}{e^{k}} ,

where F^{-1} is the inverse of the incomplete gamma quantile function. Is there a way to implement this in stan?

The function I need exists in the boost library of C++ as gamm_q_inv(). If there is no way to implement this in stan, how can I call the existing function in C++ to use within stan? I have come across some tutorials on using C++ libraries in stan, I must confess they are not the friendliest. I will also appreciate an easy to understand tutorial on using C++ libraries in stan if that is the only way out.

Thanks in anticipation of the replies.

1 Like
2

Hi,

Can you clarify exactly how you are parametrising the Gamma here? What does F(d, be^k) mean ?
A Gamma distribution is indexed by two parameters. If we call them \alpha and \beta, I’d write

\operatorname{Pr}(X \leq x) = F(x ; \alpha, \beta) = p,

for X \sim \operatorname{Gamma}(\alpha, \beta).

3

Hello,

I cleared it a bit in the expression for b. F(d, be^k) is actually F(be^k, 1, d). Apologies for the mix-up though.

4

But if the shape parameter is 1, then you’re dealing with an exponential random variable, for which one can write the CDF in closed-form. That’s why it’s important to clarify the exact parametrisation.

5

No, the rate parameter is 1. The shape parameter reamins d. Think of \frac{\gamma(d, be^k)}{\Gamma(d)} = \frac{q}{c-1}, where \gamma(.) is the lower incomplete gamma function. The LHS is nothing but the gamma CDF with shape parameter d and rate parameter b I think. However, we want to solve for b.

6

Right. I understand now. It’s somewhat uncommon to write the rate parameter before the shape parameter. As you note, one route would be to implement your own C++ that calls this function. I’m not the best to advise on that, but @andrjohns and @bbbales2 can probably help.

Another route would be to implement a somewhat quick and dirty approximation. Something like the asymptotic approximations given in this paper.

7

@fomotis which interface are you using to call Stan? The method for including external c++ differs a little between them

8

I am using rstan.

9

I started writing this up before I realised you were calling a function in Boost. This adds a layer of complexity, unfortunately. Stan’s sampler depends on calculating the gradients for each function. Generally, Stan’s autodiff functionality can calculate these automatically, but this doesn’t tend to work with external libraries like Boost. Consequently, you need to manually specify the derivatives wrt to each input for Boost functions.

I’m not familiar with the inverse of the incomplete gamma, but do you have a reference for the gradients? (or @maxbiostat by any chance?), then I can specify those as well

1 Like
10

That can be arranged. First, recall that, for X \sim \operatorname{Gamma}(\alpha, \beta) with shape and rate parametrisation, we have

\operatorname{Pr}(X \leq x) = F_X(x) = \frac{\gamma(\alpha, \beta x)}{\Gamma(\alpha)}.

Now, what we really want is the derivative of

F^{-1}(p) = x, 0 < p < 1.

This will be, by the inverse function theorem,

\frac{d}{dp}F_X^{-1}(p) = \frac{1}{F_X^\prime\left(F_X^{-1}(p)\right)} = \frac{1}{F_X^\prime\left(x\right)}.

However, F_X^\prime(x) is nothing but the probability density function, thus

\begin{align} \frac{d}{dp}F_X^{-1}(p) &= \left(\frac{\beta^\alpha}{\Gamma(\alpha)} \left(F_X^{-1}(p) \right)^{\alpha-1}\exp\left\{-\beta F_X^{-1}(p) \right\}\right)^{-1},\\ &= \frac{\Gamma(\alpha)}{\beta^\alpha} \left(F_X^{-1}(p) \right)^{-\alpha+1}\exp\left\{\beta F_X^{-1}(p) \right\}. \end{align}

WARNING: it is way past my bedtime, so please triple-check.

1 Like
11

Realised I made you do a bunch of work without checking myself first, sorry! Turns out I can get the formula through WolframAlpha embarassingly easily.

When I run it through WolframAlpha I get the following derivative:

\frac{\partial }{\partial x}\left( Q^{-1}(\alpha,x) \right) = \Gamma{(\alpha)}Q^{-1}(\alpha,x)^{1-\alpha}\left( -\exp{\{Q^{-1}(\alpha,x)}\} \right)

Given that according to Boost’s code, the inv_gamma_q(a,x) function is equivalent to the inverse CDF of the gamma distribution with a rate parameter of 1:

Q^{-1}(\alpha,x) = F^{-1}(x;\alpha,1)

The derivative is the same as your derivation (no pun intended) where the rate \beta is 1, except that the exponential should be negated.

Sorry again for the unnecessary work!

1 Like
12

@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

4 Likes
13

Also, I’ll work on getting this function added for a future release, so it will be easier to use (and more optimised)

4 Likes
14

Hello @andrjohns, @maxbiostat

I actually have a script that works only within the generated quantities block using the gamma_p_inv() (Incomplete Gamma Functions - 1.53.0) . I believe it doesn’t work within the transformed parameters block, because of the required derivative. I think we should be differentiating the inverse of the CDF of F, but I think that’s what you both did though.

15

@andrjohns thanks for writing this up! Very minor modification to get this to work using cmdstanr on linux. I had to remove the scal.hpp reference and put in #include <stan/math/prim.hpp> as this was changed from rstan. After doing that and calling with

mod <- cmdstan_model(stan_file, include_paths = getwd(),
                     cpp_options = list(USER_HEADER = file.path(getwd(), hpp_file)),
                     stanc_options = list("allow-undefined")
)

works!

However, trying it on macOS 11.4 with intel I’m getting a bunch of errors related to Undefined symbols for architecture x86_64. Any idea why?

16

I made the following changes based on my last comment. Could you please have a look at this?

** gamma_icdf.hpp** C++ file to compute F^{-1}(x, \alpha, \beta)

#ifndef GAMMA_ICDF_HPP
#define GAMMA_ICDF_HPP

#include <boost/math/special_functions/gamma.hpp>

double gamma_icdf(const double &p, const double &alpha, const double &beta, 
                  std::ostream *pstream__) {
  if(p <= 0) {
    return 0.0; 
  } else if(p >= 1) {
    return positive_infinity();
  } else {
    return boost::math::gamma_p_inv(alpha, p) / beta;
  }
}

#endif

gamma_p_inv.hpp C++ File

#include <boost/math/special_functions/hypergeometric_pFq.hpp>
#include <stan/math/prim/scal.hpp>
#include <stan/math/rev/core.hpp>
#include "gamma_icdf.hpp"

// Function to handle Stan parameters
inline stan::math::var gamma_p_inv_stan(const stan::math::var& p,
                                        const stan::math::var& alpha,
                                        const stan::math::var& beta, 
                                        std::ostream* pstream__) {
  // Extract values (not derivatives of inputs)
  double p_val = p.val();
  double alpha_val = alpha.val();
  double beta_val = beta.val();
  
  // Calculate inverse of the CDF of the gamma function
  //double gamma_p_inv_val = boost::math::gamma_q_inv(a_val, q_val);
  double gamma_p_inv_val = gamma_icdf(p_val, alpha_val, beta_val);
  
  // Pre-compute some quantities for the gradients
  double exp_gamma_p_inv_val = std::exp(gamma_p_inv_val);
  double pow_gamma_p_inv_val = std::pow(gamma_p_inv_val,(1 - alpha_val));
  double gamma_alpha_val = stan::math::tgamma(alpha_val);
  
  // Calculate derivative wrt to 'alpha' parameter
  double alpha_d = exp_gamma_p_inv_val
    * pow_gamma_p_inv_val
    * (stan::math::square(gamma_alpha_val)
    * std::pow(gamma_p_inv_val,alpha_val)
    // Boost doesn't have a normalised Hypergeometric pFq, implement our own
    *  (boost::math::hypergeometric_pFq({alpha_val,alpha_val},
                                        {alpha_val+1,alpha_val+1},
                                        -gamma_p_inv_val)
    / stan::math::square(stan::math::tgamma(alpha_val+1)))
    + (beta_val - 1)
    * gamma_alpha_val
    * std::log(gamma_p_inv_val)
    + boost::math::polygamma(0, alpha_val)
    * (gamma_alpha_val - boost::math::tgamma(alpha_val, 
                                             gamma_p_inv_val)));
    
    // Calculate derivative wrt 'q' parameter
    double beta_d = gamma_alpha_val * -exp_gamma_p_inv_val * pow_gamma_p_inv_val;
    
    // Return new parameter with values of the function and gradients
    return stan::math::var(new stan::math::precomp_vv_vari(gamma_p_inv_val, 
                                                           alpha.vi_, beta.vi_, 
                                                           alpha_d, beta_d));
}
17

I get the following error message when I compile the simple stan code provided by @andrjohns

Error in compileCode(f, code, language = language, verbose = verbose) : 
                   from file6a5c525817f6.cpp:14:~/R/win-library/4.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:~/R/win-library/4.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:13: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function] static void set_zero_all_adjoints() {             ^~~~~~~~~~~~~~~~~~~~~make: *** [C:/PROGRA~1/R/R-devel/etc/x64/Makeconf:244: file6a5c525817f6.o] Error 1
In addition: Warning message:
In readLines(file, warn = TRUE) :
  incomplete final line found on ~\test_gammaNew.stan'
Error in sink(type = "output") : invalid connection
18

Just to double-check what your overall aim is here - you need the CDF and the quantile function (inverse CDF) of the Gamma distribution, is that right?

19

Yes. I only need the the quantile of the gamma distribution. I already have a C++ script for computing the F^{-1}(x, \alpha, \beta), i.e. the quantile. I already pasted that up there. I try to read through your codes and make adjustments where I think is necessary using the function from the script I provided. I provided the whole code to that above as well, but running it I encounter the error message in my last post.

20

Gotcha. Yep the code I posted will only be valid for the gamma_q_inv function, since the derivatives are specific to that and only having two arguments.

Let me look up the gamma_p_inv implementation and the derivatives when using it for the quantile