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)