My final goal is to build a hierarchical model of logistic population growth models for multiple fish species. The problem I have now is that the model work “ok” independently, but can not work if I combine them together.

The error I get is as follows:

"Chain 1: Rejecting initial value:

Chain 1: Error evaluating the log probability at the initial value.

Chain 1: Exception: model474863befb1__namespace::log_prob: P[sym1__, sym2__] is 0.0001, but must be greater than or equal to 0.001 (in ‘string’, line 24, column 1 to column 47)

Chain 1:

Chain 1: Initialization between (-2, 2) failed after 100 attempts.

Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”

[1] “error occurred during calling the sampler; sampling not done” "

Notice: The time-series data have missing entries, which are labeled as “0” and filtered out in the model section in stan. The demo data are provided.

Independent stan model is as follows

```
data{
int<lower=1> N;
vector[N] Catch;
vector[N] CPUE;
vector[N] Stocking;
vector[N] Environment;
real area;
}
parameters{
real logK;
real logR;
real logQ;
real logSigma;
real impact_E;
real<lower=0.001,upper=1.5> P0;
}
transformed parameters{
real<lower=0> k;
real<lower=0> r;
real<lower=0> q;
real<lower=0> sigma;
vector<lower=0>[N] P;
vector[N] CPUE_hat;
k = exp(logK);
r = exp(logR);
q = exp(logQ);
sigma = exp(logSigma);
P[1] = P0;
CPUE_hat[1] = k*q*P[1];
for (t in 2:N){
P[t]=P[t-1]+(r*exp(impact_E*Environment[t-1]))*P[t-1]*(1-P[t-1])-Catch[t-1]/(k*area)+Stocking[t-1]/(k*area);
if(P[t]<0.0001){
P[t]=0.0001; // prevent negative P
} else {
P[t]=P[t];
}
CPUE_hat[t] = k*area*q*P[t];
}
}
model{
logR ~ uniform(log(0.1),log(0.5));
logK ~ uniform(log(1),log(1000));
logQ ~ uniform(log(0.00001),log(1));
for (t in 1:N){
if (CPUE[t]>0){
log(CPUE[t]) ~ normal(log(CPUE_hat[t]),sigma); // filter out the missing entries
}
}
}
```

Combined model is as follows.

```
data{
int Nstocks;
int Nyears;
matrix[Nyears,Nstocks] CPUE;
matrix[Nyears,Nstocks] Catch;
matrix[Nyears,Nstocks] Stocking;
matrix[Nyears,Nstocks] Environment;
real Area [Nstocks];
int Nsurvey_year[Nstocks];
}
parameters{
vector[Nstocks] logK;
vector[Nstocks] logR;
vector[Nstocks] logQ;
vector[Nstocks] logSigma;
vector[Nstocks] impact_E;
vector<lower=0.001,upper=2>[Nstocks] P0;
}
transformed parameters{
vector[Nstocks] k;
vector[Nstocks] r;
vector[Nstocks] q;
vector[Nstocks] sigma;
matrix<lower=0.001,upper=2>[Nyears,Nstocks] P;
matrix[Nyears,Nstocks] CPUE_hat;
for (i in 1:Nstocks){
k[i] = exp(logK[i]);
r[i] = exp(logR[i]);
q[i] = exp(logQ[i]);
sigma[i] = exp(logSigma[i]);
P[1,i] = P0[i];
for (o in 2:Nsurvey_year[i]){
P[o,i] = P[o-1,i]+(r[i]*exp(impact_E[i]*Environment[o-1,i]))*P[o-1,i]*(1-P[o-1,i])-Catch[o-1,i]/(k[i])+Stocking[o-1,i]/(k[i]);
if (P[o,i]<0.0001){
P[o,i] = 0.0001;
} else {
P[o,i] = P[o,i];
}
CPUE_hat[o,i] = k[i]*P[o,i]*q[i];
}
}
}
model{
for (i in 1:Nstocks){
logR[i] ~ uniform(log(0.1),log(0.5));
logK[i] ~ uniform(log(1),log(1000));
logQ[i] ~ uniform(log(0.00001),log(1));
for (o in 1:Nsurvey_year[i]){
if (CPUE[o,i]>0){
log(CPUE[o,i]) ~ normal(log(CPUE_hat[o,i]),sigma[i]);
}
}
}
}
```

combine_data.RData (1.3 KB)