Exception: normal_lpdf: Location parameter is nan, but must be finite! in Cmdstan

Modeling
1

Hi all,

I have been running the following model with Cmdstan. I keep on having the following warning:

Chain 1 Exception: normal_lpdf: Location parameter is -inf, but must be finite!

OR:

Chain 1 Exception: normal_lpdf: Location parameter is -inf, but must be finite!

Initially, I wanted to account for a lognormal distribution in the likelihood, since the error message belongs to that line (112). I have tried too different initial values or more constrained prior distributions to avoid the -inf or rare values. However, I keep on having this warning continuously.
Do you have any suggestions on how to avoid this?
My initial values look like this:

rAB0 = rep(0, participants)
rmuab = rep(0, participants)
r = cbind(rAB0, rmuab)
initebl5 = function(){
list(logmus = 1.03, logmul = 1.04 , logphis = 0.3,
logphil = 0.3,logmuab = 0.02,
logAB0 = 5, logsigmaAB0 = 0.5, logsigmamuab = 0.5,
sigma = 1, r = r)
}

functions {
  real[] sldecay(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[3];
    dydt[1] = -theta[1]*y[1];
    dydt[2] = -theta[2]*y[2];
    dydt[3] = theta[3]*y[1] + theta[4]*y[2] - theta[5]*y[3];
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  real ts[N];
  real t0;
}
transformed data {
  real x_r[0];
  int x_i[0]; 
}
parameters {
  real logmus;
  real logmul;
  real logphis;
  real logphil;
  real logmuab;
  
  real logAB0;
  real <lower=0> logsigmaAB0;
  real <lower=0> logsigmamuab;
  real <lower= 0> sigma;
  //real rho_rhobit;
  vector[2] r[nsubjects];
}
model {
  real new_mu[N,3];
  real mu[1,3];
  real eval_time[1];
  real theta[5];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  
  real inits[3];
  
  real mus = exp(logmus);
  real mul = exp(logmul);
  real phis = exp(logphis);
  real phil = exp(logphil);
  real muab = exp(logmuab);
  real AB0 = exp(logAB0);
  real sigmaAB0 = exp(logsigmaAB0);
  real sigmamuab = exp(logsigmamuab);
  
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  //real rho = (exp(rho_rhobit) - 1)/(exp(rho_rhobit) + 1);
  
  //prior distributions
  logmus ~ normal(0, 1);
  logmus ~ normal(0, 1);
  logphis ~ normal(0, 1);
  logphil ~ normal(0, 1);
  logmuab ~ normal(0, 1);
  logAB0 ~ lognormal(5, 10);
  sigma ~ normal(1, 10);
  logsigmaAB0 ~ normal(0, 1);
  logsigmamuab ~ normal(0, 1);
  //rho_rhobit ~ uniform(-10,10);
  
  mean_vec[1] = -pow(sigmaAB0, 2.0)/2;
  mean_vec[2] = -pow(sigmamuab, 2.0)/2;
  
  Sigma[1,1] = pow(sigmaAB0, 2.0);
  Sigma[2,2] = pow(sigmamuab, 2.0);
  Sigma[1,2] = 0;
  Sigma[2,1] = 0;
  
  theta[1] = mus;
  theta[2] = mul;
  theta[3] = phis;
  theta[4] = phil;
  
  
  for (subj_id in 1:nsubjects) {
    r[subj_id] ~ multi_normal(mean_vec, Sigma);
    rAB0[subj_id] = r[subj_id,1];
    rmuab[subj_id] = r[subj_id,2];
  }
  
  for (obs_id in 1:N){
    inits[3] = AB0*exp(rAB0[subj_ids[obs_id]]);
    inits[1] = inits[3]/mus;
    inits[2] = inits[3]/mul;
    
    theta[5] = muab*exp(rmuab[subj_ids[obs_id]]); 
    
    eval_time[1] = ts[obs_id];
    if (eval_time[1] == 0){
      new_mu[obs_id,1] = inits[1];
      new_mu[obs_id,2] = inits[2];
      new_mu[obs_id,3] = inits[3] - pow(sigma, 2.0)/2;
    }
    else{
      mu = integrate_ode_rk45(sldecay, inits, t0, eval_time, theta, x_r, x_i);
      new_mu[obs_id,1] = mu[1,1];
      new_mu[obs_id,2] = mu[1,2];
      new_mu[obs_id,3] = mu[1,3] - pow(sigma, 2.0)/2;
    }  
    //likelihood function for AB levels
    y[obs_id] ~ normal(new_mu[obs_id,3], sigma); 
  }
}
generated quantities {
  real z_pred[N];
  real sigma2 = sigma;
  for (t in 1:N){
    z_pred[t] = normal_rng(log(y[t] + 1e-5), sigma2); 
  }
}