Gamma distribution quantile and inverse quantile function

General
, , , ,
21

From my understanding, I think the derivative you did is for F^{-1}(x, \alpha, \beta). If this is the case, then your derivation should be correct.

22

You can have a look at this https://www.boost.org/doc/libs/1_46_1/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/gamma_dist.html

23

I was misunderstanding the gamma_q_inv function, it gives:

P\left( X > x \right)

Whereas the gamma_p_inv function gives:

P\left( X \leq x \right)

Which is what you were after. This is why my derivatives were different from @maxbiostat (Luiz was right all along!)

So if we denote the inverse of the regularized lower incomplete gamma function (gamma_q_inv) as Q^{-1}(\alpha,x), the inverse of the regularized upper incomplete gamma function (gamma_p_inv) is given by Q^{-1}(\alpha,1-x).

Then, the quantile function is given by (as you noted):

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

The partial derivatives for this are slightly different than those I posted earlier (plus thereā€™s now a partial for \beta), so I need to update the code a little

1 Like
24

Hmm interesting, I had it running with cmdstanr on Windows as well. What are the errors?

I think Iā€™ve got an OSX vm I can spin up for testing somewhere

25

Ah nvm, there was some issue with previous code. Checked again and it works :)

1 Like
26
\frac{\delta F^{-1} (x; \alpha, \beta)}{\delta \beta} = -\frac{ Q^{-1} (x, \alpha, \beta)^{1-\alpha} (Q^{-1} (x, \alpha, \beta)^\alpha + \beta \Gamma(\alpha) exp(Q^{-1} (x, \alpha, \beta)) ) }{\beta^2}

link to derivative

The partial derivative with respect to \alpha is the same as the previous divided by \beta. I will try to make the changes if you agree with the derivatives, and report back if I still get the same errors.

EDIT: @maxbiostat edited the LaTeX to display correctly.

1 Like
27

That wolfram code isnā€™t what youā€™ve specified in the LaTeX, unfortunately.

InverseGammaRegularized[a,q] is the same as gamma_q_inv(a,q), which isnā€™t want you want.
Additionally, InverseGammaRegularized[a,q]/q is the inverse divided by x, not by \beta.

The correct derivative is given by: D[InverseGammaRegularized[a,1-q]/b,q] - Wolfram|Alpha

But as I mentioned before your changes to the code arenā€™t valid for handling the three inputs, which is also why you experienced compilation errors. Just let me code up the partials and I can give you the example for handling three inputs

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
29

Same error as before. I made no changes to the codes.

Error in compileCode(f, code, language = language, verbose = verbose) : 
                   from file6a5c79a83cd3.cpp:14:C:/Users/lucp9544/Documents/R/win-library/4.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:C:/Users/lucp9544/Documents/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: file6a5c79a83cd3.o] Error 1
Error in sink(type = "output") : invalid connection
30

Can you post the stan_model call that youā€™re using and the result of calling cat on the Stan and C++ files?

31

So after some more testing, it looks like including the Boost header (hypergeometric_pFq.hpp) is causing issues with rstanā€™s usage of other boost functions.

The best way forward here is to use cmdstanr: R Interface to CmdStan ā€¢ cmdstanr, which includes external c++ in a different way.

Once you have cmdstanr installed, the following script should work in including external c++:

library(cmdstanr)

hpp = "#include <stan/math/prim.hpp>
#include <stan/math/rev/core.hpp>
#include <boost/math/special_functions/hypergeometric_pFq.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));
}
"

writeLines(hpp,con="gamma_icdf.hpp")

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);
}"

writeLines(model,con="model_gamma.stan")

mod = cmdstan_model("model_gamma.stan", include_paths = getwd(),
                     cpp_options = list(USER_HEADER = file.path(getwd(), "gamma_icdf.hpp")),
                     stanc_options = list("allow-undefined"))

mod$sample(data=list(y_mean=0))

Note that you donā€™t need to specify your model/hpp as strings first, Iā€™ve just done it for reproducibility

EDITED: First half of script didnā€™t copy

1 Like
32

After several attempts I finally have cmdstan and cmdstanr installed. I tried running the codes below but ended up with Error 127 mingw32-make.exe: *** [make/program:47: C:\Users\lucp9544\AppData\Local\Temp\Rtmpa05r64\model-4714e0f2622.hpp] Error 127. I tried the bernoulli example which works fine.

mod = cmdstan_model(stan_file = "test_gammaNew.stan", 
                    include_paths = getwd(),
                    cpp_options = list(USER_HEADER = 
                                         file.path(getwd(), 
                                                   "gamma_icdf.hpp")),
                    stanc_options = list("allow-undefined"))

test_gammaNew.stan is exactly as model_gamma.stan. The same for gamma_icdf.hpp.
From some check on Google, I figured it has something to do with paths but I canā€™t seem to figure it out.
Also, the code works fine without the external C++ file.

33

What happens if you run the script I posted above with no changes?

34

Same error as I posted earlier

mingw32-make.exe: *** [make/program:47: C:\Users\lucp9544\AppData\Local\Temp\RtmpKCGLcy\model-62479601194.hpp] Error 127
35

Hmm, @spinkney would you mind running the above script to see if it works for you? Just to see whether itā€™s only my installation that it works for

36

Hello @andrjohns ,
Any update? I still tried again with no success. If I run the model without the C++ file, it works just fine. Do you also use a windows system?

37

Good timing! I was just looking into this. It looks like the error has something to do with how cmdstanr handles Windows paths. I can replicate your error, and then fix it if I use the Windows short path name.

Can you try the cmdstan_model call again, but this time wrap the getwd() calls in a call to utils::shortPathName:

mod = cmdstan_model("model_gamma.stan", include_paths = utils::shortPathName(getwd()),
                    cpp_options = list(USER_HEADER = file.path(utils::shortPathName(getwd()), "gamma_icdf.hpp")),
                    stanc_options = list("allow-undefined"))
1 Like
38

Hello @andrjohns , it works with no errors, but if I attempt to use the gamma_icdf function in any stan model, I get the following errors.

undefined reference to `boost::math::tools::promote_args<double, double, double, float, float, float>::type gamma_model1_model_namespace::gamma_icdf<double, double, double>(double const&, double const&, double const&, std::ostream*)'

undefined reference to `boost::math::tools::promote_args<double, double, double, float, float, float>::type gamma_model1_model_namespace::gamma_icdf<double, double, double>(double const&, double const&, double const&, std::ostream*)'

Here is the stan model I used in testing;

functions {
  //real gamma_icdf(real p, real alpha, real beta);
  real gamma_icdf(real x, real alpha, real beta);
}
data {
  int<lower = 1> K; 
  vector[K] p;
  real<lower = 0> alpha;
  real<lower = 0> beta;
}
parameters {
}
generated quantities {
  vector[K] q;
  for(k in 1:K) {
    //q[k] = gamma_icdf(p[k], alpha, beta);
    q[k] = gamma_icdf(p[k], alpha, beta);
  }
}

I also tested with another program where I called the gamma_icdf() function in the generated quantities block, I also ended with the same error. Here is the stan program for that

functions {
  real gamma_icdf(real p, real k, real theta);
}
data {
  int<lower = 1> N;
  int<lower = 1> K; 
  vector[N] x;
  vector[K] p;
  real q;
}
parameters {
  real<lower = 0> alpha;
  real<lower = 0> beta;
}
transformed parameters{
  
  real b;
  b=gamma_icdf(q, alpha, beta);
  
}
model {
  alpha ~ normal(10, 10);
  beta ~ normal(5, 5);
  x ~ gamma(alpha, beta);
}
1 Like
39

Ah, the gamma_icdf implementation I gave you will only work when all inputs are parameters sorry. Iā€™ll modify it to work for all combinations and post an updated version

1 Like
40

Hereā€™s a gamma_icdf that has overloads for all combinations of data and parameter:
gamma_icdf.hpp (9.3 KB)

3 Likes