[ADMB Users] something suspicious in my TPL codes (semi-implicit)

[Saang-Yoon]

Dr. Fournier.  Thank you very much for your help and time.   Now the
estimation improved very much thanks to your advice, however, the
estimate of K is stuck to the upper boundary of K, and some estimate
of Ft (e.g., Ft of 1993) is unreasonably large.  I show corrected TPL
codes and data (DAT file) below.     Would you mind suggesting what I
need to do for avoid the problem?   Best Wishes,  Saang-Yoon


DATA_SECTION
  int i;  //for loop;
  int j;  //for loop;

  init_vector years(1,3);
  int nyrs;
  int ncol;  //data column
  !!nyrs=years(2)-years(1)+1;
  !!ncol=size_count(years);

  init_matrix yieldsurvey(1,nyrs,1,ncol);
  vector yrs(1,nyrs);
  vector Yt(1,nyrs);
  vector It(1,nyrs);

  !!yrs=trans(yieldsurvey)(1);
  !!Yt=trans(yieldsurvey)(2);  //yield;
  !!It=trans(yieldsurvey)(3);

  number minK;
  number maxK;
  !!minK=max(Yt)-max(Yt)*15/100;
  !!maxK=minK+minK*3000/100;

  number minY;
  number maxY;
  number meanY;
  !!minY=min(Yt);
  !!maxY=max(Yt);
  !!meanY=mean(Yt);

  !!cout<<"yrs: "<<yrs<<endl;
  !!cout<<"Yt: "<<Yt<<endl;
  !!cout<<"It: "<<It<<endl;
  !!cout<<"K range: "<<minK<<" "<<maxK<<endl;


PARAMETER_SECTION
  init_bounded_number q(0.0,1.0,1);
  init_bounded_number r(0.0,5,2);
  init_bounded_number K(minK,maxK,3);   //aka B_infinitive
  init_bounded_number b1(0.0,1.0,4);     //B1=b1*K;
  init_bounded_dev_vector F_devs(1,nyrs,-5.0,5.0,5);
  //init_number log_avg_F;
  init_bounded_number log_avg_F(-5.0,0.4,6);

  number B1;
  vector Bt(1,nyrs+1);
  vector Ft(1,nyrs);

  vector predYt(1,nyrs); //expected yield in year t
  vector predIt(1,nyrs);

  vector sig2It(1,nyrs);

  vector meanlogIt(1,nyrs);
  number varlogI;

  number varlogY;

  number Bpen;
  number bio_tmp; //biomass temporary
  number y_tmp; //yield temporary

  number Bmsy;
  number msy; //msy;
  number Fmsy;

  objective_function_value f;

PROCEDURE_SECTION
  B1=b1*K;

  f=0.0;
  Bpen=0;
  //cout<<"f1: "<<f<<endl;

  //Project B_t
  Bt(1)=B1;
  Bt(1)=posfun(Bt(1),Yt(1),Bpen); //keeps it above

  //Ft
  for(i=1;i<=nyrs;i++)
    Ft(i)=exp(log_avg_F+F_devs(i));


  //four time steps per year desired
  int nsteps=20;
  double delta=1.0/nsteps;
  //integrate the logistic dynamics over n time steps per year
  for(i=1;i<=nyrs;i++)   {
     bio_tmp=1.e-20+Bt(i);
	 y_tmp=0.0;  //yield
	 for(j=1;j<=nsteps;j++) {
	    bio_tmp=bio_tmp*(1.0+r*delta)/(1.0+(r*bio_tmp/K+Ft(i))*delta);
	    y_tmp=y_tmp+Ft(i)*bio_tmp*delta;
     };
     predYt(i)=y_tmp; //expected yield;
     Bt(i+1)=bio_tmp;
     Bt(i+1)=posfun(Bt(i+1),meanY,Bpen);  //keeps it above
  };
  f+=Bpen;

  //
  //survey as ROUTINE;
  for(i=1;i<=nyrs;i++)
	 predIt(i)=q*Bt(i);
  varlogI=norm2(log(It)-log(1.e-10+predIt))/nyrs;
  f+=0.5*nyrs*log(varlogI)+norm2(log(It)-log(1.e-10+predIt))/
(2.0*varlogI);
  cout<<"f3: "<<f<<endl;

  //
  //commerical yield as ROUTINE;
  varlogY=norm2(log(Yt)-log(1.e-10+predYt))/nyrs;
  f+=0.5*nyrs*log(varlogY)+norm2(log(Yt)-log(1.e-10+predYt))/
(2.0*varlogY);
  cout<<"f4: "<<f<<endl;


  //management references
  Bmsy=K/2;
  msy=r*K/4;  //MSY;
  Fmsy=msy/Bmsy;


REPORT_SECTION
  report<<"yr Yt It Bt predIt varlogI Ft Ft/Fmsy (Bt-Yt) f Bpen
max.grad"<<endl;
  for(i=1;i<=nyrs;i++)   {
    report<<yrs(i)<<" "<<Yt(i)<<" "<<It(i)<<" "<<Bt(i)<<"
"<<predIt(i)<<" ";
    report<<varlogI<<" "<<Ft(i)<<" "<<Ft(i)/Fmsy<<" ";
    report<<Bt(i)-Yt(i)<<" "<<f<<" "<<Bpen<<"
"<<objective_function_value::gmax<<endl;
  };


GLOBALS_SECTION
  #include <admodel.h>
  #include <math.h>
  #include <stdio.h>
  #include <stddef.h>
  #include <stdlib.h>



#==================== DAT (data file below)
#====================
#years
1968 2010 -99
#
#yieldsurvey:	data	on	yield	and	survey values
#commerical	yield	in	mt/1000
#survey	index:	relative	total	biomass	in mt/1000
#yr	yield	survery
1968	18.321	10.227
1969	21.271	9.519
1970	21.41	4.833
1971	15.61	6.178
1972	18.039	6.142
1973	16.953	6.299
1974	17.211	3.561
1975	16.75	2.257
1976	14.988	1.463
1977	10.639	2.699
1978	6.944	2.274
1979	6.935	1.45
1980	7.539	6.412
1981	6.979	2.5
1982	12.52	2.203
1983	11.989	2.068
1984	6.28	0.576
1985	3.267	0.688
1986	3.474	0.796
1987	3.58	0.494
1988	2.759	0.165
1989	1.783	0.948
1990	4.089	0.703
1991	2.564	0.708
1992	5.299	0.559
1993	4.3	0.529
1994	4.158	0.871
1995	1.135	0.344
1996	1.7	1.265
1997	2.464	3.67
1998	3.985	4.22
1999	4.963	7.738
2000	7.341	5.666
2001	7.419	11.213
2002	5.663	3.644
2003	6.562	3.919
2004	6.815	4.966
2005	3.851	2.391
2006	2.109	4.388
2007	1.662	7.912
2008	1.504	6.9
2009	1.806	6.797
2010	1.16	2.242


On Apr 1, 4:55 pm, dave fournier <da... at otter-rsch.com> wrote:
> This looks bad to me.
>
>   //Ft
>    for(i=1;i<=nyrs;i++)
>      Ft(i)=Ft(i)*exp(F_devs(i));
>
> I think you should have something like
>
>     init_number log_avg_F
>
>   //Ft
>    for(i=1;i<=nyrs;i++)
>      Ft(i)=exp(log_avf_F+F_devs(i));
> _______________________________________________
> Users mailing list
> Us... at admb-project.orghttp://lists.admb-project.org/mailman/listinfo/users