# Boost function in rstan

**URL:** <https://discourse.mc-stan.org/t/boost-function-in-rstan/16290>\
**Category:** RStan\
**Created:** [June 29, 2020, 8:37pm UTC](https://discourse.mc-stan.org/t/boost-function-in-rstan/16290 "2020-06-29T20:37:54Z")\
**Posts on this page:** 1\
**Showing post:** 10

<div class="post-metadata">

**Author:** ![jstag](https://avatars.discourse-cdn.com/v4/letter/j/9dc877/32.png) [@jstag](https://discourse.mc-stan.org/u/jstag)\
**Post date:** [July 1, 2020, 10:13pm UTC](https://discourse.mc-stan.org/t/boost-function-in-rstan/16290/10 "2020-07-01T22:13:09Z")

</div>

I have tried several approaches, but am now thinking it might be easiest if I follow the suggestion here, avoid Boost altogether and use the algebraic solver to calculate:

> [@Quantile function gamma revisited](https://discourse.mc-stan.org/t/quantile-function-gamma-revisited/5468/4):
>
> Formule: Γ^-1(p; λ^-2), with mu 0 and sigma 1 and p the percentage you want to compute the time for and lambda sometimes expressed as Q (third parameter generalized gamma) I hope this helps.

I have coded an extremely silly toy example to test out this approach and it appears to work. I’m sure there are many more efficient ways to do this (vectorization, etc.), but this appears to get the job done for now. Any improvements would be much appreciated.

```no-highlight
require(rstan)
n <- 500
rand_y <- rgamma(n, shape = 0.5, scale = 5)

fit_data <- list(N = n, y = rand_y)

mc <- 
'
functions {
vector qgamma(vector y, // unknowns
              vector theta, // parameters
              real[] x_r, // data (real)
              int[] x_i) { // data (integer)
  vector[1] z;
  z[1] = gamma_cdf(y, theta[1], 1/theta[2]) - x_r[1];
  return z;
}
}
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real<lower=0> shape;
real<lower=0> scale;
}
model {
shape ~ gamma(0.5,0.2);
scale ~ gamma(0.5,0.2);

for (n in 1:N)
  y[n] ~ gamma(shape, 1/scale);
}
generated quantities {
vector[1] qgamma_result;
vector[2] theta = [shape, scale]';
vector[1] y_guess = [0.5]';

int x_i[0];
real x_r[1]; // Return the 40th percentile

x_r[1] = 0.4;

qgamma_result = algebra_solver(qgamma, y_guess, theta, x_r, x_i);
}
'

mod <- stan_model(model_code = mc, data = fit_data, iter = 2000, init =list(list(shape = 0.5, scale = 10)), chains = 1)

print(fit)
plot(fit)

require(dplyr)
est_df <- data.frame(shape = rstan::extract(fit, "shape")$shape, scale = rstan::extract(fit, "scale")$scale, qgamma_draw = rstan::extract(fit, "qgamma_result")$qgamma_result)
est_df <- est_df %>% mutate(qgamma_correct = qgamma(0.4, shape = shape, scale = scale))
head(est_df)

```

---

_[View the full topic](https://discourse.mc-stan.org/t/boost-function-in-rstan/16290)._