Dear all,
I’m new to RStan, and I’m writing my spatial survival models as following codes. I find similar method to my codes for vectorizing integrate_1d in RStan(file:///G:/Dr%20Nazar/UBUNTU%20SERVER%20-%20UJI/stan%20files/Vectorizing%20integrate_1d%20-%20General%20-%20The%20Stan%20Forums.html). Now, I would like to speed this process up because it longs about 14 days for N=35000 cases, but I get as error as follows:
"
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: Error in function tanh_sinh::integrate: The tanh_sinh quadrature evaluated your function at a singular point and got nan. Please narrow the bounds of integration or check your function for singularities. (in ‘model38147f0761c4_Onicesco_interval_independentmodel_Log_logistic34744’ at line 129)
[1] “Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: Error in function tanh_sinh::integrate: The tanh_sinh quadrature evaluated your function at a singular point and got nan. Please narrow the bounds of integration or check your function for singularities. (in ‘model38147f0761c4_Onicesco_interval_independentmodel_Log_logistic34744’ at line 129)”
[1] “error occurred during calling the sampler; sampling not done”
"
Thank you in advance!
The code I used here as below:
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;
vector <lower=0>[N] ub;
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters{
vector [num] beta;
real <lower=0> sigma;
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 top level transformed parameter
time_ratio = exp(beta);
exp_x = exp(x);
lambdaa = exp ((-1 * respon[,7:18] * beta) / sigma);
{// `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;
vector[N] method_two;
int x_idx [N]= sort_indices_asc(ub);
vector[N] sx = ub[x_idx];
vector[N] sorted_m2;
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) {
real lb;
real ub1 = sx[i];
real increment;
real prev_val;
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-method_two[i]);
if (i==1) {
lb = 0;
prev_val = 0;
} else {
lb = sx[i-1];
prev_val = sorted_m2[i-1];
}
increment = integrate_1d(integrand, lb,ub1 ,{ alpha, lambdaa[i]}, {0.0} , {0}, 1e-8);
sorted_m2[i] = prev_val + increment;
}
method_two[x_idx] = sorted_m2;
target += respon[,6] .* (log(f)-log(s)) + log(s);
}
generated quantities{
vector[N] log_lik;//log-likelihood
vector[N] method_two;
{for(n in 1:N){
log_lik[n] =respon[n,6] * (log(((exp ((-1 * respon[n,7:18] * beta) * alpha)*alpha*respon[n,4]^(alpha-1))/((1+(exp ((-1 * respon[n,7:18] * beta) *alpha)*respon[n,4]^alpha))^2)))-log(1-method_two[n])) + log(1- method_two[n]);
}
}
.```