Gamma distribution quantile and inverse quantile function

General
, , , ,
41

Hello @andrjohns ,

same error message.

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*)'
42

Ah, something else Iā€™d missed. One of the quirks of cmdstanr's approach to including external c++ is that the function needs to be wrapped in an appropriate namespace.

This namespace is given by the model name followed by _model_namespace, so the for model_gamma.stan example, the functions need to be declared within the model_gamma_model_namespace.

Hereā€™s an updated version of the header with the functions declared in the model_gamma_model_namespace:
gamma_icdf.hpp (10.0 KB)

But because that would only allow you use the gamma_icdf functions with Stan models called gamma_model.stan, hereā€™s an R function which will take a Stan model name and the path to the hpp file to update that namespace automatically:

prep_hpp = function(stan_filename,hpp_file) {
    namespace = paste0(gsub(".stan","",stan_filename),"_model_namespace")
    writeLines(gsub('[^ ]*_model_namespace',namespace,readLines(hpp_file)),
               con=hpp_file)
}

Which you would use like so:

prep_hpp("model_gamma_new.stan",file.path(getwd(),"gamma_icdf.hpp"))
mod = cmdstan_model("model_gamma_new.stan", include_paths = utils::shortPathName(getwd()),
                    cpp_options = list(USER_HEADER = file.path(utils::shortPathName(getwd()), "gamma_icdf.hpp")),
                    stanc_options = list("allow-undefined"))
2 Likes
43

Hello @andrjohns ,

Thank you so much for your support so far. Everything works fine now, but there is a lot of error during the MCMC sampling. The error message I typically get is provided below

Error in function boost::math::hypergeometric_pFq<long double>: Cancellation is so severe that no bits in the reuslt are correct

I am thinking this has to do with the model itself, but I want to be sure it isnā€™t technical.

Also, the function returns NAN if I call the gamma_icdf() function in the generated quantities block. I wanted to use this to test that the function returns the same quantiles as in qgamma() function in R. It returns the error below if I only call the function within the generated quantities block.

Chain 4 Exception: Error in function boost::math::gamma_p_inv<double>(double, double): Overflow Error  in line 30, column 4 to column 42
44

Yeah this looks like the values being passed to gamma_icdf are causing issues with Boostā€™s gamma_p_inv and hypergeometric_pFq.

Can you extract some of the values that are being passed? You can either insert print() statements to print the parameter values, or you can extract the values at each iteration from the stanfit object itself

1 Like
45

posteriorDraws.csv (917.2 KB)

Find attached the posterior draws.

c = exp(par[3])+1;
d = exp(par[4])+0.000000001;
b = gamma_icdf(q/(c-1), d, 1)*exp(-par[2]);

and q = 0.05 is supplied in the data block.

46

Having a look at those draws, it appears that all values of b are valid and correct?

Do these errors occur frequently throughout the sampling or just a few times?

47

I wouldnā€™t say it occurs frequently, but it does occur more than a few times especially between 10%- 20% sampling. A lot of the proposed values gets rejected. I am using the message from the sampler to judge that.

Also, when I call the function within the generated quantities block, it returns the following error

 Error in function boost::math::gamma_p_inv<double>(double, double): Overflow Error in line 30, column 4 to column 42

and the posterior means of the variable is always NaN. You can try the following stan code

functions {
  real gamma_icdf(real x, real alpha, real beta);
}
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);
}
generated quantities {
  vector[K] q2;
  for(k in 1:K) {
    q2[k] = gamma_icdf(p[k], alpha, beta);
  }
}

and I called it using the following R codes

prep_hpp("gamma_model2.stan",
         file.path(getwd(),"gamma_icdf.hpp"))
probs <- c(0, 0.001, 0.01, 0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975, 0.99, 0.999, 1)

mod <- cmdstan_model("gamma_model2.stan", 
                     include_paths = utils::shortPathName(getwd()),
                    cpp_options = list(USER_HEADER = 
                                         file.path(utils::shortPathName(getwd()),
                                                               "gamma_icdf.hpp")),
                    stanc_options = list("allow-undefined"))
x <- rgamma(100, shape = 10, rate = 5)
fit <- mod$sample(
  data = list(
    N = length(x),
    K = length(probs),
    p = probs,
    x = x,
    q = 0.5
  ),
  seed = 123,
  chains = 3,
  parallel_chains = 3,
  refresh = 500
)

All the values eturned for q2 are NaNs.