search

seacode / gmacs

A generic size-structured stock assessment model
https://seacode.github.io/gmacs
18 stars 14 forks source link

Add Francis re-weighting #180

Closed quantifish closed 8 years ago

quantifish commented 8 years ago

hardwired for number of indices but for output I have:

FUNCTION Write_R
  // Development--just start to get some output into R
  dvar_vector sdnr_idx(1,5);
  dvar_vector Francis_wts(1,5);
  Francis_wts(1) = franwt(oac_fsh,eac_fsh,sam_fsh); // fishery age compositions
  Francis_wts(2) = franwt(oac_bts,eac_bts,sam_bts); // bts age compositions

FUNCTION dvariable franwt(const dmatrix& obs,const dvar_matrix& pred,const dvector& Nsam)
  dvar_vector ftmp(obs.indexmin(),obs.indexmax());
  for (i=obs.indexmin();i<=obs.indexmax();i++) 
    ftmp(i) = (mn_age(obs(i)) - mn_age(pred(i))) / sqrt(square(Sd_age(obs(i),pred(i))) / Nsam(i));
  return 1./var(ftmp);

From rock lobster assessment:

FUNCTION dvector get_Francis_lf_wts()
// this code based on equation TA1.8 in Francis(2011) should be changed so separate weights if by sex
  int ireg,ii,isx,nobs,ij;
  dvector lfwt(1,nregions);
  double O_yj,v_yj;
  double P_yj;

  for (ireg=1; ireg<=nregions; ireg++)
  {
    nobs=nlfrqs(ireg);
    dvector resid(1,nobs);
    ij=1;
    resid=0;
    ofstream ofsxx(region_name(ireg)+adstring("FrancisLFs.out"));
    ofsxx<<"Year Season Type Obs Pred Sigma Wt Resid N effN"<<endl;
    for (ii=1; ii<=nlfrqs(ireg); ii++)
    {
    //     Lfsigma_risl(ireg,ii)=sigmatilde/(lfrqwt_ri(ireg,ii)*Like_wt(1));
    //     N is 1/Lfsigma
    //     Note that Lfsigma is created for each bin, but LFsigma is the same for each bin

       O_yj=0;
       P_yj=0;
       v_yj=0;
       for (isx=1; isx<=3; isx++)
       {
         if(sum(Olfrq_risl(ireg,ii,isx))>0)
         {
           O_yj+=sum(elem_prod(Olfrq_risl(ireg,ii,isx)(bins(ireg,isx,1),bins(ireg,isx,2)), Size(bins(ireg,isx,1),bins(ireg,isx,2))));
           P_yj+=sum(elem_prod(value(Plfrq_risl(ireg,ii,isx)(bins(ireg,isx,1),bins(ireg,isx,2))), Size(bins(ireg,isx,1),bins(ireg,isx,2))));
           v_yj+=sum(elem_prod(value(Plfrq_risl(ireg,ii,isx)(bins(ireg,isx,1),bins(ireg,isx,2))), square(Size(bins(ireg,isx,1),bins(ireg,isx,2)))));
          }
        }
        v_yj-=square(P_yj);

        resid(ij++)=(O_yj-P_yj)/sqrt(v_yj*value(Lfsigma_risl(ireg,ii,1,1)));
        ofsxx<<lfrqyr_ri(ireg,ii)<<" "<<lfrqseas_ri(ireg,ii)<<" "<<lfrqtype_ri(ireg,ii)<<" "<<O_yj<<" "<<P_yj<<" "<<v_yj<<" "<<sqrt(v_yj*value(Lfsigma_risl(ireg,ii,1,1)))<<" "<<resid(ij-1)<<" "<<lfrqN_ri(ireg,ii)<<" "<<1./Lfsigma_risl(ireg,ii,isx,1)<<" "<<endl;
      }
     lfwt(ireg)=1./(square(std_dev(resid))*((nobs-1.)/nobs*1.));  // nobs adjustment cause Fournier std is pop. est.
     cout<<"LF re-weight "<<lfwt(ireg)<<" ("<<lfwt(ireg)*Like_wt(1)<<")"<< endl;
  }
  return lfwt;
wStockhausen commented 8 years ago

Darcy,

You might consider doing the iterative re-weighting in the ADMB code itself by adding a couple of extra minimization phases (command line switch "maxphs", maybe?) after all parameters are turned on. Once all parameters are turned on, at the end of each subsequent phase you could re-weight the size comps using the Francis weights in the BETWEEN_PHASES_SECTION (or report section). In the next phase you'd be using the new Francis weights to calculate the objective function. Then just check that the Francis weights have converged when you get to the final phase (if not, you could then exit to R to continue re-weighting or you could re-run the model and increase what you set "maxphs" to). Seems easier than exiting to R and re-running the model with re-weighted sample sizes (unless you have to).

Buck

On Mon, May 23, 2016 at 1:43 PM, Darcy Webber notifications@github.com wrote:

hardwired for number of indices but for output I have: FUNCTION Write_R // Development--just start to get some output into R dvar_vector sdnr_idx(1,5); dvar_vector Francis_wts(1,5); Francis_wts(1) = franwt(oac_fsh,eac_fsh,sam_fsh); // fishery age compositions Francis_wts(2) = franwt(oac_bts,eac_bts,sam_bts); // bts age compositions [8:41:11 AM] Jim Ianelli: FUNCTION dvariable franwt(const dmatrix& obs,const dvar_matrix& pred,const dvector& Nsam) dvar_vector ftmp(obs.indexmin(),obs.indexmax()); for (i=obs.indexmin();i<=obs.indexmax();i++) ftmp(i) = (mn_age(obs(i)) - mn_age(pred(i))) / sqrt(square(Sd_age(obs(i),pred(i))) / Nsam(i)); return 1./var(ftmp);

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/seacode/gmacs/issues/180

wStockhausen commented 8 years ago

Just noticed that the function "franwt" returns a dvariable. You definitely want to return a double if you decide to re-weight in the BETWEEN_PHASES_SECTION. You don't want any derivative information from the calculation of "franwt"s affecting the minimization in the subsequent phase. No problems, obviously, if you're just writing to R, though.

Buck

On Wed, May 25, 2016 at 7:40 AM, William Stockhausen - NOAA Federal < william.stockhausen@noaa.gov> wrote:

Darcy,

You might consider doing the iterative re-weighting in the ADMB code itself by adding a couple of extra minimization phases (command line switch "maxphs", maybe?) after all parameters are turned on. Once all parameters are turned on, at the end of each subsequent phase you could re-weight the size comps using the Francis weights in the BETWEEN_PHASES_SECTION (or report section). In the next phase you'd be using the new Francis weights to calculate the objective function. Then just check that the Francis weights have converged when you get to the final phase (if not, you could then exit to R to continue re-weighting or you could re-run the model and increase what you set "maxphs" to). Seems easier than exiting to R and re-running the model with re-weighted sample sizes (unless you have to).

Buck

  • Dr. William T. Stockhausen *
  • Resource Ecology and Fisheries Management *
  • Alaska Fisheries Science Center *
  • National Marine Fisheries Service *
  • National Oceanic and Atmospheric Administration *
  • 7600 Sand Point Way N.E. *
  • Seattle, Washington 98115-6349 *
  • email: William.Stockhausen@noaa.gov *
  • voice: 206-526-4241 fax: 206-526-6723 *
  • web : http://www.afsc.noaa.gov *

    All models are wrong, some are useful.--G.E.P. Box Beware of geeks bearing equations. --W. Buffett

    Disclaimer: The opinions expressed above are personal and do not necessarily reflect official NOAA policy.

On Mon, May 23, 2016 at 1:43 PM, Darcy Webber notifications@github.com wrote:

hardwired for number of indices but for output I have: FUNCTION Write_R // Development--just start to get some output into R dvar_vector sdnr_idx(1,5); dvar_vector Francis_wts(1,5); Francis_wts(1) = franwt(oac_fsh,eac_fsh,sam_fsh); // fishery age compositions Francis_wts(2) = franwt(oac_bts,eac_bts,sam_bts); // bts age compositions [8:41:11 AM] Jim Ianelli: FUNCTION dvariable franwt(const dmatrix& obs,const dvar_matrix& pred,const dvector& Nsam) dvar_vector ftmp(obs.indexmin(),obs.indexmax()); for (i=obs.indexmin();i<=obs.indexmax();i++) ftmp(i) = (mn_age(obs(i)) - mn_age(pred(i))) / sqrt(square(Sd_age(obs(i),pred(i))) / Nsam(i)); return 1./var(ftmp);

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/seacode/gmacs/issues/180

quantifish commented 8 years ago

Done. Also added in sdnr and MAR. These are simply added to the gmacs.rep output file. I didn't add any ADMB code that will do the iterative re-weighting automatically at this stage (although this may become a feature in the future). The reason why I didn't do this is that data-weighting is very subjective, and can differ depending on whether or not the user chooses to use sdnr, MAR or Francis method. There is now also an example of Francis-weighting used in one the the SMBKC runs.