Boost function in rstan

Interfaces RStan
10

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:

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.

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)
1 Like