Add Francis re-weighting

[quantifish]

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]

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
#180

[wStockhausen]

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

---

• 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 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
#180

[quantifish]

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.