Dear Stan community,
I want to check my model with simulated data in R. To this end, first I generated the data in R and then run my model in the cmdstanr with 50 iterations (init=0, n=310,iterations=50,iter_warmup=4000,iter_sampling=4000, adapt_delta=0.99,max_treedepth=15,chains=1). But, I got the following error during the execution (after iterations=30). It should be noted that the model runs without any problems with a low number of iterations, for example, 5 iterations.
Can @rok_cesnovar and anyone else help me to solve this issue? I will appreciate any recommendations.
"
Model executable is up to date!
Running MCMC with 1 chain…
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 Initialization failed.
Warning: Chain 1 finished unexpectedly!
Error: No chains finished successfully. Unable to retrieve the draws.
In addition: Warning message:
Error: No chains finished successfully. Unable to retrieve the draws.
".
My stan codes are here:
functions{
// defines the log spdf
vector a (vector p,matrix kernel, vector exp_x){
vector[num_elements(p)] a_; // `a` saved to output
//{ // `z` not part of output
vector[num_elements(p)] landa;
for (k in 1 : num_elements(p)) {
landa[k] = kernel[k,] * exp_x * p[k];
}
a_ = landa / sum(landa); // assign `a`
for (i in 1:num_elements(p)) {
if (is_nan(a_[i]) || (a_[i] < 0.0)) {
reject("a_[",i,"] is undefined");
}
}
return a_;
}
// defines the log survival
vector log_S(vector t, vector p, matrix kernel, vector exp_x,int[] county, real shape, vector scale) {
vector[num_elements(t)] log_S_;
vector[num_elements(p)] spdf;
spdf = a (p,kernel, exp_x);
for (i in 1:num_elements(t)){
//log_S[i] = log(spdf[county[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)
log_S_[i] = log(spdf[county[i]]) - log1p(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,vector p, matrix kernel, vector exp_x, int[] county, real shape,vector scale){
vector[num_elements(t)] log_f_;
vector[num_elements(p)] spdf;
spdf = a (p,kernel, exp_x);
for (i in 1:num_elements(t)){
//log_f[i] = log(spdf[county[i]]) + log((scale[i]*shape*t[i]^(shape-1))/((1+(scale[i]*t[i]^shape))^2));
log_f_[i] = log(spdf[county[i]]) + log(scale[i]*shape*t[i]^(shape-1)) - 2* log1p((scale[i]*(t[i])^shape));
}
return log_f_;
}
//defines the log-likelihood
real surv_loglogistic_lpdf(vector t,vector d,vector p, matrix kernel, vector exp_x, int[] county, real shape,vector scale){
vector[num_elements(t)] log_lik;
real prob;
vector[num_elements(t)] ls_s;
vector[num_elements(t)] ls_f;
for (i in 1 : num_elements(t)) {
if (t[i] <= 0.0) {
reject("t[", i, "] is ", t[i]);
}
}
ls_s = log_S( t, p, kernel, exp_x, county, shape, scale);
ls_f = log_f( t, p, kernel, exp_x, county, shape, scale);
log_lik = d .* log_f(t, p, kernel, exp_x, county, shape, scale)+ (1 - d) .* log_S(t, p, kernel, exp_x, county, shape, scale);
for (i in 1 : num_elements(log_lik)) {
if (is_inf(log_lik[i]) || is_nan(log_lik[i])) {
reject("log_lik[", i, "] is ", log_lik[i]);
}
}
prob = sum(log_lik);
return prob;
}
}
//data block
data{
int <lower=0> m;
int <lower=0> n_grid;
int <lower=0> n;
int <lower=0> num;
int <lower=0> number;
int <lower =0 , upper = m> county[n];
vector <lower=0> [m] p;
matrix [m,n_grid] kernel;
matrix [n,number] data3;
}
//parameters block
parameters{
vector [num] beta;
real <lower=0> sigma;//scale parameter sigma=1/shape
vector [n_grid] x;
vector <lower=0> [m] r_county;
}
// transformed parameters block
transformed parameters{
vector [num] time_ratio;
vector[n] lambdaa;
real <lower=0> alpha = 1 / sigma;
vector [n_grid] exp_x;
time_ratio = exp(beta);
exp_x = exp(x);
//linpred = x*beta;//linear predictor
lambdaa = exp (((-1 * data3[,4:7] * beta) / sigma)+log(r_county [county]));
}
// model block
model{
target += normal_lpdf(x|0,5);
target += normal_lpdf(beta | 0, 1000);
target += inv_gamma_lpdf(sigma | 0.001, 0.001);
target += normal_lpdf(r_county| 0, 1);
for (i in 1:m) {
if (is_nan(r_county[i]) || (r_county[i] < 0.0)) {
reject("r_county[",i,"] is undefined");
}
}
target += surv_loglogistic_lpdf (data3[,11] | data3[,2],p, kernel, exp_x, county, alpha, lambdaa);
}
.```
Thank you!