Hi Ben @bbbales2,
Thanks for understanding the mess in my code. This is the function block i have used there.
functions{
real Q_integral(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 tn = theta[5];
real Q;
Q= (1/x)*(exp(((-(log(x)-log(tn))^2)/(2*sigma2))+((mu-log(tn))*log(x)/sigma2)-(lambda*(tn-x)^alpha)));
return(Q);
}
real W_integral(real x,
real xc,
real[] theta,
real[] x_r,
int[] x_i){
real mu = theta[1];
real sigma2 = theta[2];
real W;
W = ((0.3 * exp(-square(log(x) - mu) / (2 * sigma2)))/(sqrt(2*pi()*sigma2)*x));
return(W);
}
}
Then, I want to do several integrals through out the my code as follows:
V[1] = beta_t[2]((0.3/sqrt(2pi()sigma2))(exp(((-(mu-log(t[2]))^2)/(2*sigma2))-(((mu-log(t[2]))log(t[2]))/(sigma2)))))(integrate_1d(Q_integral,t[1],t[2],{mu, sigma2,lambda, alpha, t[2]}, {0.0}, x_i, 1e-8));
I[1] = (1-beta_t[2])(((0.3/sqrt(2pi()sigma2))(exp(((-(mu-log(t[2]))^2)/(2*sigma2))-(((mu-log(t[2]))log(t[2]))/(sigma2)))))integrate_1d(Q_integral,t[1], t[2],{ mu, sigma2, alpha, lambda, t[2]}, {0.0}, {0}, 1e-8)-
((0.3/sqrt(2pi()sigma2))(exp(((-(mu-log(t[3]))^2)/(2sigma2))-(((mu-log(t[3]))log(t[3]))/(sigma2)))))integrate_1d(Q_integral,t[1], t[2],{ mu, sigma2, alpha, lambda, t[3]}, {0.0}, {0}, 1e-8))+
integrate_1d(W_integral,t[2], t[3],{ mu, sigma2}, {0.0}, {0}, 1e-8)-
((0.3/sqrt(2pi()sigma2))(exp(((-(mu-log(t[3]))^2)/(2sigma2))-(((mu-log(t[3]))*log(t[3]))/(sigma2)))))*integrate_1d(Q_integral,t[2], t[3],{ mu, sigma2, alpha, lambda, t[3]}, {0.0}, {0}, 1e-8);
.
.
.
.
init_func <- function(){list(b0=runif(1,0,4),mu=runif(1,4,4.5),sigma2=runif(1,0.01,0.05),lambda=runif(1,0.01,0.5),alpha=runif(1,1.5,4))}
fit<-stan(model_code = cancer_code,
data = cancer_data,
iter = 20,
warmup = 10,
init = init_func,
init_r = 0.5,
chains = 2,
control = list(max_treedepth=10 , adapt_delta = 0.95)
)
I think I could understand what you suggested me to do. So, I let the code to print parameter values and limits of integration in following order
print(lower_limit, “,” , upper_limit , “,” , mu , “,” , sigma2 , “,” , alpha , “,” , lambda);
Then stan gives me:
SAMPLING FOR MODEL ‘d5fc77a7c85e8fcd0e4880ff25844065’ NOW (CHAIN 1).
Chain 1:
15,50,4.21296,0.0231476,1.67584,0.309843
15,51,4.21296,0.0231476,1.67584,0.309843
15,52,4.21296,0.0231476,1.67584,0.309843
15,53,4.21296,0.0231476,1.67584,0.309843
15,54,4.21296,0.0231476,1.67584,0.309843
15,55,4.21296,0.0231476,1.67584,0.309843
15,56,4.21296,0.0231476,1.67584,0.309843
15,57,4.21296,0.0231476,1.67584,0.309843
15,58,4.21296,0.0231476,1.67584,0.309843
15,59,4.21296,0.0231476,1.67584,0.309843
15,60,4.21296,0.0231476,1.67584,0.309843
15,61,4.21296,0.0231476,1.67584,0.309843
15,62,4.21296,0.0231476,1.67584,0.309843
15,63,4.21296,0.0231476,1.67584,0.309843
15,64,4.21296,0.0231476,1.67584,0.309843
15,65,4.21296,0.0231476,1.67584,0.309843
15,66,4.21296,0.0231476,1.67584,0.309843
15,67,4.21296,0.0231476,1.67584,0.309843
15,68,4.21296,0.0231476,1.67584,0.309843
15,69,4.21296,0.0231476,1.67584,0.309843
15,70,4.21296,0.0231476,1.67584,0.309843
15,71,4.21296,0.0231476,1.67584,0.309843
15,72,4.21296,0.0231476,1.67584,0.309843
15,73,4.21296,0.0231476,1.67584,0.309843
15,74,4.21296,0.0231476,1.67584,0.309843
Chain 1: 15,50,4.21296,0.0231476,1.67584,0.309843
Chain 1: Exception: integrate: error estimate of integral 2.84703e+13 exceeds the given relative tolerance times norm of integral (in ‘modela682a643eba_d5fc77a7c85e8fcd0e4880ff25844065’ at line 90)
[1] “Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: integrate: error estimate of integral 2.84703e+13 exceeds the given relative tolerance times norm of integral (in ‘modela682a643eba_d5fc77a7c85e8fcd0e4880ff25844065’ at line 90)”
[1] “error occurred during calling the sampler; sampling not done”
SAMPLING FOR MODEL ‘d5fc77a7c85e8fcd0e4880ff25844065’ NOW (CHAIN 2).
Chain 2:
15,50,4.47482,0.0305653,1.52133,0.321217
15,51,4.47482,0.0305653,1.52133,0.321217
15,52,4.47482,0.0305653,1.52133,0.321217
15,53,4.47482,0.0305653,1.52133,0.321217
15,54,4.47482,0.0305653,1.52133,0.321217
15,55,4.47482,0.0305653,1.52133,0.321217
15,56,4.47482,0.0305653,1.52133,0.321217
15,57,4.47482,0.0305653,1.52133,0.321217
15,58,4.47482,0.0305653,1.52133,0.321217
15,59,4.47482,0.0305653,1.52133,0.321217
15,60,4.47482,0.0305653,1.52133,0.321217
15,61,4.47482,0.0305653,1.52133,0.321217
15,62,4.47482,0.0305653,1.52133,0.321217
15,63,4.47482,0.0305653,1.52133,0.321217
15,64,4.47482,0.0305653,1.52133,0.321217
15,65,4.47482,0.0305653,1.52133,0.321217
15,66,4.47482,0.0305653,1.52133,0.321217
15,67,4.47482,0.0305653,1.52133,0.321217
15,68,4.47482,0.0305653,1.52133,0.321217
15,69,4.47482,0.0305653,1.52133,0.321217
15,70,4.47482,0.0305653,1.52133,0.321217
15,71,4.47482,0.0305653,1.52133,0.321217
15,72,4.47482,0.0305653,1.52133,0.321217
15,73,4.47482,0.0305653,1.52133,0.321217
15,74,4.47482,0.0305653,1.52133,0.321217
Chain 2: 15,50,4.47482,0.0305653,1.52133,0.321217
Chain 2: Exception: integrate: error estimate of integral 7.21561e+22 exceeds the given relative tolerance times norm of integral (in ‘modela682a643eba_d5fc77a7c85e8fcd0e4880ff25844065’ at line 90)
[1] “Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: integrate: error estimate of integral 7.21561e+22 exceeds the given relative tolerance times norm of integral (in ‘modela682a643eba_d5fc77a7c85e8fcd0e4880ff25844065’ at line 90)”
[1] “error occurred during calling the sampler; sampling not done”
I have noticed that parameter values didn’t change over the data table which contains 25 rows with different ages.
I would be grateful if anyone can comments here.
Thank you all