Gamma distribution quantile and inverse quantile function

General
, , , ,
28

Alright, so for the function:

F^{-1}\left(x; \alpha, \beta \right) = \frac{Q^{-1}(\alpha,1-x)}{\beta}

Which is given in Wolfram Alpha by InverseGammaRegularized[a,1-x]/b, and can be verified by comparing against qgamma(x,shape=a,rate=b).

The partial derivatives for each input are given by:

The c++ 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_icdf(const stan::math::var& x,
                                  const stan::math::var& alpha,
                                  const stan::math::var& beta,
                                  std::ostream* pstream__) {
  // Extract values (not derivatives) of inputs
  double x_val = x.val();
  double a_val = alpha.val();
  double b_val = beta.val();

  // Calculate inverse of the incomplete gamma function
  double gamma_p_inv_val = boost::math::gamma_p_inv(a_val, x_val);
  // Use to derive quantile of Gamma distribution
  double gamma_icdf = gamma_p_inv_val / b_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 - a_val));
  double gamma_a_val = stan::math::tgamma(a_val);

  // Calculate derivative wrt 'x' parameter
  double x_d = (1 / b_val) * gamma_a_val * exp_gamma_p_inv_val * pow_gamma_p_inv_val;

  // Calculate derivative wrt to 'alpha' parameter
  double a_d = (1 / b_val) * exp_gamma_p_inv_val
                * pow_gamma_p_inv_val
                * (stan::math::square(gamma_a_val)
                   * std::pow(gamma_p_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_p_inv_val)
                        / stan::math::square(stan::math::tgamma(a_val+1)))
                   - x_val
                   * gamma_a_val
                   * std::log(gamma_p_inv_val)
                   + boost::math::polygamma(0, a_val)
                   * (gamma_a_val - boost::math::tgamma(a_val, gamma_p_inv_val)));

  // Calculate derivative wrt to 'beta' parameter
  double b_d = -gamma_p_inv_val / stan::math::square(b_val);

  // Return new parameter with values of the function and gradients
  return stan::math::var(new stan::math::precomp_vvv_vari(gamma_icdf, x.vi_, alpha.vi_, beta.vi_, x_d, a_d, b_d));
}

With example model:

functions {
  real gamma_icdf(real x, real alpha, real beta);
}

data {
  real y_mean;
}

parameters {
  real<lower=0> y;
}

model {
  y ~ normal(y_mean, 1);
}
2 Likes