Alright so I’ve found the time to fiddle with the model. First, my results for fitting the provided posterior:
Using 6 parallel chains with 6 cores I get on my local notebook within 5m+15s (warmup+sampling) a minimal N_Eff
of ~147 with iter_sampling=84
. Diagnostics are fine (max(rhat) < 1.01, E-BFMI=.9, no divergences). During sampling the mean number of leapfrog steps is 16.
I believe the main performance gain is due to a custom warm-up procedure which not only generally performs much better than Stan’s regular warm-up, but crucially allows me to continuously adapt the “centeredness” of the parametrization, which leads to very efficient sampling. Sadly, this procedure is not yet ready for public consumption and also currently I don’t have the time to make it so.
Still, I have some other comments:
- Change your model such that positivity of all physically positive parameters is ensured. This includes the fishwise
pAm_i
andv_i
parameters, which currently may become negative. - Similarly, take care that the ODE solution only yields physically reasonable (positive) values.
- Change your ODE such that instead of states {E,V,H,UR} you have states {E,V,H+UR}. H and UR can be recovered from H+UR via H=min(Hp,H+UR) etc.
- Think about whether this fixed threshold Hp really makes sense (I don’t believe it does, even if it were a fishwise unknown parameter).
- Zero-fecundity observations can be handled via something like
0~bernoulli_logit(H+UR-Hp)
. This will encourage simulated fish to stay immature. Personally, I would change the model to get rid of the fixed threshold and incorporate something like the above. - Zero-fecundity simulations can also be handled using
H+UR-Hp
, appropriately mapped from R to R+. - Finally, it looks like you will have to play with the centeredness of your parametrization.
Edit: Updating to add that for 32 fish everything took 12m+1m.