I would recommend against using interval priors unless things have physically constrained bounds. For instance, it looks like you’re constraining standard deviations between 0- and 1 and effects between upper and lower log scale values. One thing you can check to see if this is a problem is whether probability is bunching up against the boundaries—that will cause divergences, etc. as it drives values on the unconstrained scale to +/- infinity.
For priors on your random effects, you should just be looping over N to make sure all of the variables get priors, e.g.,
log_thetaA ~ normal(log_bar_thetaA, omega_thetaA);
This produces a centered parametrization, so you need to go back and declare the variables with the bounds to get a non-centered parameterization, e.g.,
vector<offset=log_bar_thetaA, multiplier=>omega_thetaA> log_thetaA;
Then loop through the row_indx values.
I was curious about why you’re multiplying y_pred[row_indx] * tau to get the scale of the observations.
You may also be running into numerical instabilities here (and you never need to multiply by 1!):
y_pred[row_indx]
= exp(log_thetaA[i])
* (1 - ages[row_indx]^exp(log_thetaC[I])
/ (exp(log_thetaB[i])^exp(log_thetaC[i])
+ ages[row_indx]^exp(log_thetaC[I])));
I would check that the right hand side here isn’t overflowing or underflowing. Particularly, taking exp(a)^exp(b) is going to overflow or underflow whenever the absolute value of b is greater than 6 unless exp(a) is very close to 1.
The observations for y_obs are constrained to be positive so you have to actually include the bounds for the normal_lpdf to get things to normalize because the arguments are not constant. That is, you want something like:
y_obs[row_indx] ~ normal(y_pred[row_indx], y_pred[row_indx] * tau) T[0, ];
where the T provides the truncation below at 0. Otherwise, the normal’s not going to be normalized with y_obs constrained to be positive.