Using integrate1d in the function block

Modeling
1

Hey friends. I’m trying to get the following model to compile:

functions {
  real deriv_gamma_density(real x,          // Function argument
                    real xc,         // Complement of function argument
                    real[] theta,    // parameters
                    real[] x_r,      // data (real)
                    int[] x_i) {     // data (integer)
    real alpha = theta[1];
    return (1 / tgamma(alpha)) * x * exp(alpha*x) * exp(-exp(x));
  }
  real lgammaJacLinf(vector x, real alpha, real lambda, real[] x_r, int[] x_i){
    real term_1 = 0; 
    real term_3 = 0;
    real term_4 = 0;
    real vector_length = num_elements(x);
    for ( i in 1:num_elements(x) ){
      real gamma_term = log(gamma_cdf(lambda*x[i], alpha, 1));
      real integral_term = integrate_1d(deriv_gamma_density,
                             negative_infinity(),
                             gamma_term,
                             { alpha }, x_r, x_i);
      term_1 += pow((1/lambda) * ( (digamma(alpha)/tgamma(alpha)) + integral_term) / exp(gamma_lpdf(gamma_cdf(lambda*x[i], alpha, 1) | alpha, 1)),2);
      term_4 += pow(x[i], 2) / pow(lambda, 2);
      term_3 += ((1/lambda) * ( (digamma(alpha)/tgamma(alpha)) + integral_term) / exp(gamma_lpdf(gamma_cdf(lambda*x[i], alpha, 1) | alpha, 1)) ) * (-x[i] / lambda);
    }
    return pow( (1/vector_length) * ( term_1 * term_4 - pow(term_3,2) ), 1/2);
  }
}
data {
  int<lower=1> N;
  vector[N] y;
}
transformed data {
  real x_r[0];
  int x_i[0];
}
parameters {
  real<lower=0>	alpha;
  real<lower=0>	lambda;
}
model {
  y ~ gamma_lpdf(alpha,1/lambda);
  target+=lgammaJacLinf(y,alpha,lambda, x_r, x_i);
}

but I’m getting the error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
 error in 'model10b4c1a410d08_gammaFinal' at line 21, column 32
  -------------------------------------------------
    19:                              alpha, 
    20:                              x_r, 
    21:                              x_i);
                                       ^
    22:       term_1 += pow((1/lambda) * ( (digamma(alpha)/tgamma(alpha)) + integral_term) / exp(gamma_lpdf(gamma_cdf(lambda*x[i], alpha, 1) | alpha, 1)),2);
  -------------------------------------------------

PARSER EXPECTED: ","
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'gammaFinal' due to the above error.

I find this puzzling, since I’m fairly sure I’ve provided all necessary parts of that function and am thus having trouble finding the syntax error.
My questions:

  1. Is it okay to use integrate1d in the function{} block? The example here only does so in the model{} block.
  2. I’m having difficulties understanding the need for x_r and x_i. Is it okay that I’ve just sort of defined them only for use in the integrate1d call?
  3. Is there a syntax error I’m just not seeing here?

Thanks so much!

2

The integrate_1d function in the Stan language accidentally requires that you specify the relative tolerance as the last argument, which is contrary to the documentation. The default value was supposed to be sqrt(epsilon()).

2 Likes
3

Thanks! It worked when I added sqrt(machine_precision()) – I’m not sure epsilon() is a built in function…

4

That’s what I meant