nmfs-ost
[Bug]: missing values for forecast derived quantities in MCMC
Describe the bug
Running the Simple model https://github.com/nmfs-ost/ss3-test-models/tree/main/models/Simple with MCMC produced either 0 or nan values for some forecast quantities had values there for the same model. The error was introduced between 3.30.23.2 and 3.30.24.0.
To Reproduce
Run the Simple model from ss3-test-models via
ss3 -mcmc 10 -mcsave 1 ss3 -mceval
ss_summary.sso includes the following lines:
SSB_2002 5632.59 1447.68
SSB_2003 6615.59 1559.8
...
ForeCatch_2002 167.969 366.748
ForeCatch_2003 403.585 379.522
whereas in derived_posteriors.sso, the columns
| SSB_2002 | SSB_2003 | ... | ForeCatch_2002 | ForeCatch_2003 |
|---|---|---|---|---|
| 5632.59 | 0 | ... | nan | nan |
| 5632.59 | 0 | ... | nan | nan |
| 5632.59 | 0 | ... | nan | nan |
Expected behavior
Values in the forecast years similar to the MLE quantities.
Screenshots
No response
Which OS are you seeing the problem on?
No response
Which version of SS3 are you seeing the problem on?
3.30.24.0 and 3.30.24.1
Additional Context
In discussion around the ongoing widow rockfish assessment, the PFMC made a request to
The STAT include a Bayesian projection to evaluate the probability that the stock would remain above B40% in 2029, if such a method is feasible to achieve within existing workload capacity and model structure.
However, the "feasible" clause makes it easy to just not do this, or explore using an earlier version of SS3. The widow update also used recdev option 1 so additional changes to the model are required to explore MCMC.
I will look into this. Note that you can get the MLE version of the requested probability by setting up the model to give the se on the ratio of B/B40
narrowing in. the numbers-at-age are getting zeroed out in second stage of forecast in mceval
@Rick-Methot-NOAA, Thanks for digging into this. Would be good to get it fixed for other users, but your suggestion of the MLE approximation to B/B40 is a great short-term solution for this particular case (which isn't due for a few weeks anyway).
Problem was with the calculation of the new quantity, HCR_anchor The calc was inside of a if(Show_MSY==1) conditional. in mceval, show_MSY = 0, so hcr_anchor was not being calculated. The first forecast loop was OK because it is doing OFL, but the ABC_buffer needed for loop 2 was dependent on HCR_anchor, so failed its calculation. Moved the HCR_anchor calc to be above the show_MSY==1 conditional in do_forecast.
@iantaylor-NOAA Thank you for catching this and @Rick-Methot-NOAA for fixing it so quickly. @aaronberger-nwfsc and I will test this fix out for Pacific hake as well.
