How to do different function integrations for Likelihoods with Rstan

Interfaces RStan
1

Hi, This is the first time I am going to use Stan with my programming. In my case I want to do posterior sampling, but you know that first thing is to get the results for log likelihood function of my problems, but, my difficulty is how to define different integrations with likelihood function and how is it fix for the function block in Stan? I tried several ways and now my code appears like that giving error message
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable “real” does not exist.
error in ‘modelcc0302e11b1_af4dcb4b0b81933d35bb7a79daf485e8’ at line 60, column 10

58:       }
59:       
60:       real WQ_integral_t2(real x, 
             ^
61:                           real xc,

My code looks like:

real likelihood_log(vector param, matrix DataTable, int N){
    int K;
    real Q;
    real W;
    real x;
    real beta;
    real tstar;
    real b0; 
    real b1; 
    real mu; 
    real sigma2; 
    real lambda; 
    real alpha;
    //real[] t;
    vector[6] beta_t;
    int delta;

    b0=param[1];
    b1=param[2];
    mu=param[3];
    sigma2=param[4];
    lambda=param[5];
    alpha=param[6];
.
.
.
.
.delta = 1;
    for (i in 1:rows(DataTable)){
    
      vector[K] D = rep_vector(0,K);
      vector[K] I = rep_vector(0,K);
      vector[K] singlelikelihood = rep_vector(0,K);
      
      vector[6] t = rep_vector(0,K+2); 
      t[1]=0 ;
      t[2]=DataTable[i,1] ; 
      for (j in 3:6) {    
        t[j]=t[j-1]+delta;
      }

      for (j in 1:6){
       beta_t[j]=1/(1+exp(-b0-b1*(t[j]-tstar)));
      }
      
      real WQ_integral_t2(real x, 
                          real xc, 
                          real[] theta, 
                          real[] x_r, 
                          int[] x_i){
                          
           real mu = theta[1];
           real sigma2 = theta[2];
           real lambda = theta[3];
           real alpha = theta[4];
           real t1 = x_r[2];
           return (0.3*exp(-((log(x)-mu)^2)/(2*sigma2))/(sqrt(2*pi()*sigma2)*x))*(-lambda*((t1-x)^(alpha)));
      }

    real WQ_integral_t3(real x, 
                        real xc, 
                        real[] theta, 
                        real[] x_r, 
                        int[] x_i){
                        
          real mu = theta[1];
          real sigma2 = theta[2];
          real lambda = theta[3];
          real alpha = theta[4];
          real t1 = x_r[3];
          return (0.3*exp(-((log(x)-mu)^2)/(2*sigma2))/(sqrt(2*pi()*sigma2)*x))*(-lambda*((t1-x)^(alpha)));
   }

   real W_integral(real x1, 
                   real x1c, 
                   real[] theta, 
                   real[] x1_r, 
                   int[] x1_i){
                   
             real mu = theta[1];
             real sigma2 = theta[2];
             return (0.3*exp(-((log(x1)-mu)^2)/(2*sigma2))/(sqrt(2*pi()*sigma2)*x1));
  }
      
      
      D[1]=beta_t[2]*integrate_1d(WQ_integral_t2,t[1], t[2], { mu, sigma2, alpha, lambda }, {00.0,55.0,51.0,52.0,53.0,54.0} , x_i, 51e-8);
      
      I[1]=(1-beta_t[2])*(integrate_1d(WQ_integral_t2,t[1], t[2],{ mu, sigma2, alpha, lambda }, x_r=t , x_i, 1e-8)-
                          integrate_1d(WQ_integral_t3,t[1], t[2],{ mu, sigma2, alpha, lambda }, x_r=t, x_i, 1e-8))+
                          integrate_1d(W_integral,t[2], t[3],{ mu, sigma2}, x_r, x_i, 1e-8)-
                          integrate_1d(WQ_integral_t3,t[2], t[3],{ mu, sigma2, alpha, lambda }, x_r=t, x_i, 1e-8);
   
.
.

.
   log_sum_singlelikelihoods = sum(log(singlelikelihood));
    }  
    final_log_likelihood = sum(unlist(log_sum_singlelikelihoods));
   return(final_log_likelihood);
}
}
2

That looks like a very complicated model.

Copy an example from the manual and change things one at a time until you’re at the model you want to be at.

That might sound like tedious, way to program, but you kinda have to work in small pieces like this anyway or the statistics will get confused anyway (models usually work worse than you expect :D).

1 Like
3

Your partial code isn’t clear to me but it seems you are trying to define one function inside another? Also, that error tends to be thrown when you try to define a variable after performing calculations. For example, defining j in the function below will throw that error:

functions {
  real fun1() {
    int i = 0;
    i = i + 1;
    real j;
    
    return 0;
  }
}

because I’m trying to define a variable too late.

1 Like
4

Yes, It seems very complex but anyhow I wanna do it. Thanks for the guidance and your ideas in this regard. Now I am modifying it in small pieces using an example.

Thanks :)

1 Like
5

Hello there, As it is a very complex and long code, I copied here some parts of my original code. :)
It is very important me to have your ideas how to define variables within functions.

Thanks for the guidance. I am just trying to improve it in small pieces at once.

6

You can also do good debugging functions like this using expose_stan_functions in Rstan. That exposes functions you wrote in Stan as R functions, and so you can compare them against reference implementations.

7

Wow, It’s a good idea and will follow.
Hope to update my outcomes.

Many thanks
Erandi