flr / FLa4a

The repository for the JRC initiative on stock assessment
https://fishreg.jrc.ec.europa.eu/web/a4a
12 stars 6 forks source link

a4a.tpl "fixes" in admb_dev branch #72

Closed jimianelli closed 9 years ago

jimianelli commented 9 years ago

Worked through a couple of the problem examples...specifically the shketest and ple4test which had different but similar issues. This can be seen at 9a533685c13779b50ffb16502e11158f957db9c4 basically, the key code mod is shown below

    // init_vector qpar(1,noQpar,1)           
    init_bounded_vector qpar(1,noQpar,-10,10,1)        
    // init_vector vpar(1,noVpar,2)        
    init_bounded_vector vpar(1,noVpar,-10,10,2)        

resulted in converged models in both pathological examples--unfortunately results and interpretation becomes invalid since they end up at their bounds (some elements do).

So the question is how to resolve. Well the reason q or variance terms would go large or small on their own is likely a scaling problem so ideally, something in the input data should check that scales are reasonable numerically for computational purposes (incidentally, I changed all instances of exp() to mfexp() to give some robustness when say, the optimization tests something ridiculous (e.g., exp(58)). Perhaps nudging these arbitrary bounds to be -15,15 or a little higher would result in enough latitude for the estimation not to fall over, but then I'm less certain about how the vpar vector in particular should be scaled...normally I would put things in the standard deviation scale for what seem to be log-scale parameters[?]...

In the shketest figure below you can see the 2nd vpar parameter settled on the bound (and has a stupidly tiny std error--because of being on the bound).

image

Another route would be to put a reasonable penalty on these types of parameters until the last phase (which might be one more than the maximum phase).

Haven't thought much more about "optimum" phasing here but it did strike me as odd to have the vpar in the 1st phase, seems having a base "ballpark" sigma of say 0.2 to start out would be reasonable, then estimate the values based on the data in later phases (easy to implement).

Anyway, all for today. I put the test directories in the inst/admb folder of the admb_dev branch (need to delete later; just put them for examining). The kludgy steps I followed to run a test example were:

  1. navigate to the FLa4a/inst/admb directory in a terminal
  2. compile a4a.tpl: admb a4a
  3. copy the exe into a folder to test, e.g.: copy admb.exe shketest
  4. navigate to shketest cd shketest
  5. run the model a4a -nox (-nox option distills what's most essential in checking convergence progress, imho).
ejardim commented 9 years ago

Hi Jim,

Sorry for the late reply. I'm on a sequence of meetings ...

The input data is centered before being passed to ADMB. Need to check if and how the scaling is done. The issue with bounds has been a problem since the beggining. What we're fitting are the the coefficients of the different submodels, not directly the values like q etc. So the bounds are not so clear.

IMO a smart phasing would be a better option, I did my best but I'm far from being an ADMB expert. So if you have any recommendations it would be very easy to implement.

Something else that is puzzling me is that when I fit a model to simulated data from a previous fit most of the times the fit doesn't converge. See the example

library(FLa4a) data(ple4) data(ple4.indices) fit <- sca(ple4, ple4.indices, fit="assessment") fit <- sca(ple4 + simulate(fit, 1), ple4.indices, fit="assessment") Warning message: In a4aInternal(fmodel = fmodel, qmodel = qmodel, srmodel = srmodel, : Hessian was not positive definite.

The simulation is quite simple. It uses a normal distribution with the covariance matrix and the vector of parameters from ADMB, to simulate a new set of parameters for the submodels. The predict method does the computation of the new catch-at-age matrix.

Best

EJ

On 11/06/15 02:42, Jim Ianelli wrote:

Worked through a couple of the problem examples...specifically the shketest and ple4test which had different but similar issues. modifying the code below

// init_vector qpar(1,noQpar,1) init_bounded_vector qpar(1,noQpar,-10,10,1) // init_vector vpar(1,noVpar,2) init_bounded_vector vpar(1,noVpar,-10,10,2)

resulted in converged models in both pathological examples--unfortunately results and interpretation becomes invalid since they end up at their bounds (some elements do).

So the question is how to resolve. Well the reason q or variance terms would go large or small on their own is likely a scaling problem so ideally, something in the input data should check that scales are reasonable numerically for computational purposes (incidentally, I changed all instances of |exp()| to |mfexp()| to give some robustness when say, the optimization tests something ridiculous (e.g., exp(58)). Perhaps nudging these arbitrary bounds to be -15,15 or a little higher would result in enough latitude for the estimation not to fall over, but then I'm less certain about how the vpar vector in particular should be scaled...normally I would put things in the standard deviation scale for what seem to be log-scale parameters[?]...

In the shketest figure below you can see the 2nd vpar parameter settled on the bound (and has a stupidly tiny std error--because of being on the bound).

image https://cloud.githubusercontent.com/assets/2715618/8097448/31922980-0f97-11e5-8ff4-d1b7aec2ef4e.png

Another route would be to put a reasonable penalty on these types of parameters until the last phase (which might be one more than the maximum phase).

Haven't thought much more about "optimum" phasing here but it did strike me as odd to have the vpar in the 1st phase, seems having a base "ballpark" sigma of say 0.2 to start out would be reasonable, then estimate the values based on the data in later phases (easy to implement).

Anyway, all for today. I put the test directories in the inst/admb folder of the admb_dev branch (need to delete later; just put them for examining). The kludgy steps I followed to run a test example were:

  1. navigate to the FLa4a/inst/admb directory in a terminal
  2. compile a4a.tpl: |admb a4a|
  3. copy the exe into a folder to test, e.g.: |copy admb.exe shketest|
  4. navigate to shketest |cd shketest|
  5. run the model |a4a -nox| (-nox option distills what's most essential in checking convergence progress, imho).

— Reply to this email directly or view it on GitHub https://github.com/flr/FLa4a/issues/72.

ejardim commented 9 years ago

Coming back to this again. Scaling seems o be the problem here, although is a difficult one to sort out. Note that we're dealing with coefficients of the linear model, not the values themselves. The scale of the parameters will depend on the submodel being used. If it's a numeric "year" the beta is in a completely different scale from a factor(year) or a smoothing parameters.

For now I added the possibility to choose which "fleet" (piece of information) to be centered, which improved the fit in some cases, where scaling the surveys was doing worst.