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.
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
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
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
Ah nvm, there was some issue with previous code. Checked again and it works :)
1 Like
\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}
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
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
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:
- x: D[InverseGammaRegularized[a,1-x]/b,x] - Wolfram|Alpha
- \alpha: D[InverseGammaRegularized[a,1-x]/b,a] - Wolfram|Alpha
- \beta: D[InverseGammaRegularized[a,1-x]/b,b] - Wolfram|Alpha
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
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
Can you post the stan_model call that youāre using and the result of calling cat on the Stan and C++ files?
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
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.
What happens if you run the script I posted above with no changes?
Same error as I posted earlier
mingw32-make.exe: *** [make/program:47: C:\Users\lucp9544\AppData\Local\Temp\RtmpKCGLcy\model-62479601194.hpp] Error 127
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
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?
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
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
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
Hereās a gamma_icdf that has overloads for all combinations of data and parameter:
gamma_icdf.hpp (9.3 KB)
3 Likes