Error evaluating the log probability at the initial value

Modeling
, , , ,
1

Hi all,

I have some issues in evaluating the log probability at the initial value.

rt[]    nt[]      rc[]      nc[]
10	595.2	21	640.2
2	762.0	0	756.0
54	5635.0	70	5600.0
47	5135.0	63	4960.0
53	3760.0	62	4210.0
10	2233.0	9	2084.5
25	7056.1	35	6824.0
47	8099.0	31	8267.0
43	5810.0	39	5922.0
25	5397.0	45	5173.0
157	22162.7	182	22172.5
92	20885.0	72	20645.0
")
names(dat) <- gsub(".", "", names(dat), fixed = TRUE)
dataList = list(
  Ntotal = 12,
  nt=log(dat$nt/1000),
  rt=dat$rt,
  nc=log(dat$nc/1000),
  rc=dat$rc
)

################Stan###############

modelString1 = "
data{
  int<lower=0>    Ntotal;
  real  nt[Ntotal];
  int<lower=0>    rt[Ntotal];
  real   nc[Ntotal];
  int<lower=0>    rc[Ntotal];
}



parameters {
    real<lower=-10,upper=10> mu;
    real<lower=0,upper=10> tau2;
    real<lower=0,upper=10>  tauphi;
    real<lower=-10,upper=10> muphi;
    real<lower=-10,upper=10> beta;
    real phi[Ntotal];
    real thetaadjpred[10];
}
transformed parameters {
    real  mt[Ntotal];
    real  mc[Ntotal];
    real <lower=0> mtt[Ntotal];
    real <lower=0> mcc[Ntotal];
    real theta[Ntotal];
    real thetaadj[Ntotal];

for(i in 1:12){
    theta[i] = thetaadj[i] + beta*(phi[i]-mean(phi[1:12]));
    
    mt[i] = nt[i] + phi[i] + theta[i];
   
    mc[i] = nc[i] + phi[i];
    mtt[i] = exp(mt[i]);//+ (1e-10)
    mcc[i] = exp(mc[i]);//+ (1e-10)
   
}
}

model
{

for(j in 1:10){
  thetaadjpred[j]  ~ normal(mu,tau2);
}
  
  for(i in 1:12){
    target += normal_lpdf(thetaadj[i] | mu,tau2);
    phi[i] ~ normal(muphi, tauphi);
    rt[i]~poisson(mtt[i]);
    rc[i]~poisson(mcc[i]);
 
  }
  

}
generated quantities {
    
    
    real phi0;
    real <lower=0> invph0;
    real thetapred[10];
    real <lower=0> rr[Ntotal];
    real <lower=0> rbase[Ntotal];
    
  phi0 = -mu/beta + mean(phi[1:12]);
  invph0=exp(phi0);

  for(j in 1:10){
    thetapred[j] = thetaadjpred[j] + beta*( (j-3)/2.0  -mean(phi[1:12]));

  }
    for(i in 1:10){
      rr[i]= exp( theta[i]);
    rbase[i]=exp(phi[i]);
    }
}

"
writeLines( modelString1 , con="TEMPmodel.stan" )
stanModel <- stan(
  data = dataList ,
  file = "TEMPmodel.stan",  
  chains = 4,            
  warmup = 5000,
  iter = 25000,
  control = list(adapt_delta = 0.99)
  
)

stanModel

When I try to debug, I try to swith the order of the code, but it never work.

here are the issues I have:
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: mtt[i_0__] is nan, but must be greater than or equal to 0 (in ‘model27e02255220e_TEMPmodel’ at line 24)

Thanks for helping to look into this.

2

The problem probably is that you run

for(i in 1:12){
   ....
  mtt[i] = exp(mt[i]);
  ...
}

i.e. you always initialize 12 elements of mtt, but mtt has Ntotal elements, so if Ntotal > 12 there will be uninitialized (and thus nan) elements in mtt.

Does that make sense?

3

Another issue is that you never initialize thetaadj. Maybe it should be declared in parameters instead of transformed parameters?

2 Likes
4

I think I assign the Ntotal as 12, why this could be an issue??

5

If you always initialize it to 12, then yes it is not an issue (although it is bad practice, writing it as for(i in 1:Ntotal) is less likely to cause issues as you change your code). It was just something that jumped at me at first glance. @nhuurre noticed another issue that is not affected by this - generally you need to take care to initialize everything except the variables in the data and parameters block.

Best of luck with your model.