Exception: poisson_lpmf: Rate parameter is nan, but must not be nan!

Dear all,

I have a time series model (see below), with parameter k and f[5]. The stan model checking was is syntactically correct. I started with giving a fixed initial value, just to try whether model works. Somehow when I run the model, I have the error below and couldnt figure out why the likelohood does not get any values. Can someone help me ?

thanks

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Rate parameter is nan, but must not be nan! (in ‘model3700927c4_eel_stan_model_01’ at line 137)

data {
  int<lower=0> Nyear;
  int<lower=0> Nage;
  // projection matrix from natural mortality and maturity
  matrix[Nage, Nage] FemMat;
  matrix[Nage, Nage] MaleMat;
  // recruitment indices at age 0
  matrix[1, Nyear] FemRec;  
  matrix[1, Nyear] MaleRec;
  // fisheries selectivity 
  real FishSelectFem[Nage];
  real FishSelectMale[Nage];
  // observed survey indices
  int Ydat[Nage,Nyear];
}

// parameters: f, continuous variable (5 periods) and k (ratio, >0)
parameters {
  vector[5] f;
  real<lower=0> k;
}

transformed parameters {
  // extend from f to full time series
  real f1[22];
  real f2[10];
  real f3[6];
  real f4[3];
  real f5[10];
  real FishMort[Nyear];
  real FemFish[Nyear,Nage,Nage];  
  real MaleFish[Nyear,Nage,Nage];
  // projected number 
  matrix[Nage, Nyear] FemNumbers;
  matrix[Nage, Nyear] MaleNumbers;
  matrix[Nage, Nyear] Numbers;
  
  f1 = rep_array(f[1],22);
  f2 = rep_array(f[2],10);
  f3 = rep_array(f[2],6);
  f4 = rep_array(f[2],3);
  f5 = rep_array(f[2],10);
  FishMort = append_array(f1, append_array(f2, append_array(f3, append_array(f4, f5))));
  
  // calculate FemFish from FishMort
  FemFish = rep_array(0,Nyear,Nage,Nage);
  
  for (i in 1:Nyear) {
    for (j in 1:Nage) {
      FemFish[i,j,j]= exp(-FishMort[i]*FishSelectFem[j]);
    }
  }
  
  //
    MaleFish = rep_array(0,Nyear,Nage,Nage);
  
  for (i in 1:Nyear) {
    for (j in 1:Nage) {
      MaleFish[i,j,j]= exp(-FishMort[i]*FishSelectMale[j]);
    }
  }
  
  ///////// intermediate calculations /////////////////////////
    
    // first year age 0 number
  FemNumbers[1, 1] = FemRec[1,1]*k;
  MaleNumbers[1, 1] = MaleRec[1,1]*k;
  
  // calculate projection for female 
  for (i in 2:Nyear) {
    // convert FemFish (array) of i year into a matrix
    matrix[Nage,Nage] iFemFish_mat;
    for (j in 1: Nage) {
      for (l in 1:Nage) {
        iFemFish_mat[j,l] = FemFish[i-1,j,l];
      }
    }
    FemNumbers[,i]= FemMat*iFemFish_mat*FemNumbers[,i-1]; // projection from i-1 year to i year
    FemNumbers[1,i]=FemRec[1,i]*k; //ith year, age 0 from recruitment
  }
  // calculate projection for male 
  for (i in 2:Nyear) {
    // convert MaleFish (array) of i year into a matrix
    matrix[Nage,Nage] iMaleFish_mat;
    for (j in 1: Nage) {
      for (l in 1:Nage) {
        iMaleFish_mat[j,l] = MaleFish[i-1,j,l];
      }
    }
    MaleNumbers[,i]= MaleMat*iMaleFish_mat*MaleNumbers[,i-1];
    MaleNumbers[1, i] = MaleRec[1,i]*k;
  }
  
  // calculate total number
  Numbers = MaleNumbers + FemNumbers;
  
}


model {
  
  //////////////////////////////////
    // prior
  f ~ uniform(0, 100);
  k ~ uniform(0, 100);
  
  // likelihood
  for (i in 1:Nage) {
    for (j in 1:Nyear) {
      Ydat[i,j] ~ poisson(Numbers[i,j]);
    }
  }
  
}

init = function() list(
f=c(0.1, 0.1, 0.1, 0.1, 0.1),
#f=rnorm(0.1,0.5,0.5),
#k=rnorm(0.1,0.5,0.5),
k=0.1
)

fit <- stan(file=“mystan.stan”, data=mydata, iter=10, chains = 1)

I’m not sure if it’s the problem but looks like you need a lower bound for f.

parameters {
  vector<lower=0>[5] f; // lower=0
  real<lower=0> k;
}
2 Likes

Thanks @nhuurre! Yes, you spotted one issue that I need to set a lower bound for f. The main issue i think is that I did not assign initial zero values for the 2 matrix:

FemNumbers = rep_matrix(0, Nage, Nyear);

MaleNumbers = rep_matrix(0, Nage, Nyear);

After I included these statements, the model is running.

Thanks

3 Likes