Error message: Initialization between (-2, 2) failed after 100 attempt

Modeling
, ,
1

Dear Stan developers,

Hi! I am fitting four parameters with Matt trick.
The four parameters are 1) s_dm_early, 2) s_dm_late, 3) s_em_early, 4) s_em_late.
When I ran the following Stan code in R, errors appeared as:

Chain 3:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 3:   Stan can't start sampling from this initial value
Chain 3: Initialization between (-2, 2) failed after 100 attempts. 
Chain 3:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."

The Stan code is as follows:

data {
  int<lower=1>                 nSub;
  int<lower=1>                 maxT;
  int<lower=0>                 SigmaBound;
  int<lower=1>                 T[nSub];            
  int<lower=0,upper=1>         Missing_dm[nSub,maxT];
  int<lower=0,upper=1>         Missing_em[nSub,maxT];
  real<lower=0,upper=1000>     st[nSub,maxT];    
  int<lower=-21,upper=1000>    rf[nSub,maxT];   
  real<lower=0,upper=1000>     dur_dm[nSub,maxT];    
  real<lower=0,upper=1000>     dur_em[nSub,maxT];   
  int<lower=0,upper=1000>      dm[nSub,maxT];    
  real<lower=-180,upper=1000>  em[nSub,maxT];    
}

transformed data {
}

parameters {
  vector[4] mu_p;
  vector<lower=0>[4] sigma;
  
  vector[nSub] s_dm_early_pr;
  vector[nSub] s_dm_late_pr;
  vector[nSub] s_em_early_pr;
  vector[nSub] s_em_late_pr;
}

transformed parameters {
  vector<lower=0,upper=SigmaBound>[nSub] s_dm_early;
  vector<lower=0,upper=SigmaBound>[nSub] s_dm_late;
  vector<lower=0,upper=SigmaBound>[nSub] s_em_early;
  vector<lower=0,upper=SigmaBound>[nSub] s_em_late;
  
  for (iSub in 1:nSub) {
    s_dm_early[iSub] = Phi_approx( mu_p[1] + sigma[1] * s_dm_early_pr[iSub]) * SigmaBound;
    s_dm_late[iSub] = Phi_approx( mu_p[2] + sigma[2] * s_dm_late_pr[iSub]) * SigmaBound;
    s_em_early[iSub] = Phi_approx( mu_p[3] + sigma[3] * s_em_early_pr[iSub]) * SigmaBound;
    s_em_late[iSub] = Phi_approx( mu_p[4] + sigma[4] * s_em_late_pr[iSub]) * SigmaBound;
  }
}

model {
  mu_p  ~ normal(0,1);
  sigma ~ normal(0,1);
  
  s_dm_early_pr ~ normal(0,1);
  s_dm_late_pr ~ normal(0,1);
  s_em_early_pr ~ normal(0,1);
  s_em_late_pr ~ normal(0,1);
  
  for (iSub in 1:nSub) {
    int inT = T[iSub];
    for (t in 1:inT) {
      if (Missing_dm[iSub,t]==0){
        
        real ip;
        real is_dm;
        real is_em;
        
        // DM 
        if (dur_dm[iSub,t] == 4.5){
          is_dm = s_dm_early[iSub];
        }
        else {
          is_dm = s_dm_late[iSub];
        }
        ip = normal_cdf(rf[iSub,t], 0, is_dm);
        dm[iSub,t] ~ bernoulli(ip);
        
        // EM
        if (Missing_em[iSub,t]==0){
          if (dur_em[iSub,t] == 4.5){
            is_em = s_em_early[iSub];
          } 
          else {
            is_em = s_em_late[iSub];
          }
          
          if (dm[iSub,t]==0){
            target += normal_lpdf(em[iSub,t] | 0, sqrt(is_dm^2 + is_em^2))
                    + normal_lcdf(em[iSub,t] - (is_dm^2 + is_em^2) * rf[iSub,t] / (is_dm^2)
                                 | 0, is_em);
          }
          if (dm[iSub,t]==1){
            target += normal_lpdf(em[iSub,t] | 0, sqrt(is_dm^2 + is_em^2))
                    + normal_lcdf(-(em[iSub,t] - (is_dm^2 + is_em^2) * rf[iSub,t] / (is_dm^2))
                                   | 0, is_em);
         } 
        }
      }
    }
  }
}

generated quantities {
}

You help will be really grateful to me!

2

Probably nSub is a large number and there are large regions of the parameter space with vanishingly low probability (hence the log(0), i.e. negative infinity). If you know this model can fit the data or got it working at least once, it indicates a problem with finding an initial value, so that’s the (relatively easy) problem to solve.

If you can try making nSub small for testing purposes and see if it manages to initialize the chain, or alternatively, if you can choose initial parameters init (maybe inits in rstan, I’m not sure) arguments passed on to the sampling function that will produce a nonzero likelihood (even if it’s a very bad fit, as long as it’s nonzero HMC will take it from there and move towards higher probability regions of parameter space).

If this doesn’t work it may indicate a problem or bug with the model, but that will require more investigation (and maybe some more detail about what the model is and what the functions in your code are supposed to do).