Problem with application of integrate_1d

General
, ,
1

Dear all,
In my model, I need to compute an integral (using integrate_1d) with the same function with several different upper limits. But, I encountered with following error in rstan:

’
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Expressions for elements of array must have the same or promotable types; found first element type=real; type at position 2=vector error in ‘model390c40e3ef1_Onicesco_interval_independentmodel_Log_logistic34744’ at line 70, column 74

70:     log_integ[i]=log(1-integrate_1d(integrand, lb, ub[i],{ alpha, lambdaa[i]}, {0.0} , {0}));
  71:     }

'.

This is my codes in the RSTAN . Is there any idea how to solve this problem?

functions{
  real integrand(real u, real xc, real[] theta, real[] x_r, int[] x_i)
  {     // data (integer)
 
   real df_dx = ((theta[2] *theta[1]*u^(theta[1]-1))/((1+(theta[2] *u^theta[1]))^2));
return(df_dx);
}
}
data {
  int <lower=0> province;
  int <lower=0> n_grid;
  int <lower=0> N;
  int <lower=0> num;
  int <lower=0> number;
 int <lower =0 , upper = province> county[N];
vector <lower=0> [province] p;
row_vector[n_grid] kernel[province];
matrix [N,number] respon;
real lb;
}
transformed data {
 real x_r[0];
 int x_i[0];
 vector <lower=lb>[N] ub;
 for (i in 1:N){
 ub[i]=respon[i,4];
    }
    }
parameters{
  vector [num] beta;
  real <lower=0> sigma;
  //real <lower=0> sigma1;
  vector [n_grid] x;
}
transformed parameters{
  vector [num] time_ratio;
  vector [N] lambdaa;
  real alpha = 1 / sigma;
  vector [n_grid] exp_x;
  vector[province] a; // `a` saved to output because it's a toplevel transformed parameter
 vector[N] log_integ;
 time_ratio  = exp(beta);
exp_x = exp(x);  
 lambdaa = exp ((-1 * respon[,7:18] * beta) / sigma);
for (i in 1:N) {
 log_integ[i]=log(1-integrate_1d(integrand, lb, ub[i],{ alpha, lambdaa[i]}, {0.0} , {0}));
    }
{ // `z` not part of output
vector[province] landa;
 for (k in 1 : province) {
 landa[k] = kernel[k] * exp_x * p[k];
         }
a = landa / sum(landa); // assign `a`
         }
         }
model{
  vector [N] f;
  vector [N] s;
 target += normal_lpdf(x|0,5);
 target += normal_lpdf(beta | 0, 1000);
 target += gamma_lpdf(alpha | 0.001, 0.001);
for (i in 1:N) {
f[i]= a[county[i]]*((lambdaa [i]*alpha*respon[i,4]^(alpha-1))/((1+(lambdaa [i]*respon[i,4]^alpha))^2));
//s[i]= a[county[i]]*(1- ((lambdaa [i]*respon[i,4]^alpha)/ (1+(lambdaa [i]*respon[i,4]^alpha)))); s[i]= a[county[i]]  * (log_integ[i]);
 target += respon[i,6] *  (log(f[i])-log(s[i])) +  log(s[i]);
 }
}
.```
2

Your model code compiled successfully for me after I explicitly specified the default relative tolerance:

log_integ[i]=log(1-integrate_1d(integrand, lb, ub[i],{ alpha, lambdaa[i]}, {0.0} , {0}, 1e-8));
1 Like
3

Thanks a ton for your reply and help. Actually, I’m new to Stan and I had trouble running my hierarchical survival model (about ~ 8 days for iter=10000, warm up=5000, N=35000, province=31, num=12).

How can I run my model using within-chain parallelization method(25.1 Reduce-sum | Stan User’s Guide)? Can you help me to change my codes in RSTAN?

Thank you in advance!

functions{
   real icar_normal_lpdf(vector phi, int province, int[] node1, int[] node2) {
    return -0.5 * dot_self(phi[node1] - phi[node2]);
  }
  // defines the log survival
  vector log_S (vector t, real shape,vector scale){
    vector[num_elements(t)] log_S;
    for (i in 1:num_elements(t)){
      log_S[i] =  log(1- ((scale [i]*(t[i])^shape)/ (1+(scale [i]*(t[i])^shape)))); // county is a vector that shows the location of patients with dim = num_elements(t) 
    }
    return log_S;
  }
  //defines the log pdf
  vector log_f (vector t,real shape,vector scale){
    vector[num_elements(t)] log_f;
    for (i in 1:num_elements(t)){
      log_f[i] =log((scale[i]*shape*t[i]^(shape-1))/((1+(scale[i]*t[i]^shape))^2));
    }
    return log_f;
  }
  //defines the log-likelihood for right censored data
  real surv_loglogistic_lpdf(vector t,vector d, real shape,vector scale){
    vector[num_elements(t)] log_lik;
    real prob;
    log_lik = d .* log_f(t, shape, scale)+ (1 - d) .* log_S(t, shape, scale);
    prob = sum(log_lik);
    return prob;
  }
   real partial_sum(array[] int y_slice,
                   int start, int end,vector d,
                   real shape, vector scale) {
    return surv_loglogistic_lpdf(y_slice | d, shape, scale);
  }
}
//data block
data{
  int <lower=0> province;
  int <lower=0> N;
  int <lower=0> num;
  int <lower=0> number;
  int <lower =0 , upper = province> county[N];
  matrix [N,number] respon;
   int<lower=0> N_edges;
  int<lower=1, upper=province> node1[N_edges];  // node1[i] adjacent to node2[i]
  int<lower=1, upper=province> node2[N_edges];  // and node1[i] < node2[i]
}
//parameters block
parameters{
  vector [num] beta;
  real <lower=0> sigma;//scale parameter sigma=1/shape
  vector <lower=0> [province] sigma1;
  vector[province] w1;  // random effect (CAR)
  vector[province] theta;       // heterogeneous effects  
}
// transformed parameters block
transformed parameters{
  vector [num] time_ratio;
  //vector[N] linpred;
  vector[N] lambdaa;
  real <lower=0> alpha = 1 / sigma;
time_ratio  = exp(beta);
  lambdaa = exp ((-respon[,7:18] * beta - w1[county] - theta[county]) / sigma);
}
// model block
model{
  target += normal_lpdf(beta | 0, 1000);
  //target += student_t_lpdf(sigma| 3,0,1); 
  target += gamma_lpdf(alpha | 0.001, 0.001);
   w1 ~ icar_normal_lpdf(province, node1, node2);
  sum(w1) ~ normal(0, 0.001 * province);
  target += normal_lpdf(theta| 0, sigma1);
  target += student_t_lpdf(sigma1| 3,0,1); 
  respon[,4] ~ surv_loglogistic//model for the data
}
// generated quantities block
generated quantities{
  vector[N] log_lik;//log-likelihood
  
  { for(n in 1:N){
  log_lik[n] = respon[n,6] * (log((exp ((-respon[n,7:18] * beta - w1[county[n]] - theta[county[n]])*alpha)*alpha*respon[n,4]^(alpha-1)) / ((1+(exp ((-respon[n,7:18] * beta - w1[county[n]] - theta[county[n]])*alpha)*(respon[n,4])^(alpha)))^2))) + (1-respon[n,6]) * ( log(1- ((exp ((-respon[n,7:18] * beta - w1[county[n]] - theta[county[n]])*alpha)*(respon[n,4])^alpha) / (1+(exp ((-respon[n,7:18] * beta - w1[county[n]] - theta[county[n]])*alpha)*(respon[n,4])^alpha)))));}
  }
}
.```