seacode / gmacs

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

CPT issues checklist #170

Closed quantifish closed 8 years ago

quantifish commented 8 years ago
quantifish commented 8 years ago

Jies objective functions look rather different to ours

FUNCTION calculate_obj_function
 int j, i, ii;
 LogLike.initialize();

 // Loglikelihoods (less additive constants)

 // 1a. Retained catch number of "legals"

// LogLike(1) = -0.5*norm2(log(x_ret + 0.001) - log(X_ret + 0.001));

// 1b. Retained catch biomass LogLike(1) = -0.5*norm2(log(x_ret_b + 0.001) - log(X_ret_b + 0.001));

// 2a. Trawl suvey abundance lognormally distributed about predicted value //LogLike(2) = -0.5*norm2(elem_div(log(x_ts)-log(X_ts),sig_ts));

// 2b. Trawl survey biomass lognormally distributed about predicted value LogLike(2) = -0.5*norm2(elem_div(log(b_ts)-log(B_ts),sig_ts_b));

// 3. Pot survey abundance lognormally distributed about predicted value LogLike(3) = -0.5*norm2(elem_div(log(x_ps)-log(X_ps),sig_ps));

// 4. Trawl survey proportions are multinomial wrt predicted proportions // LogLike(4) = effn_ts_rowsum(elem_prod(p_ts,log(P_ts+0.01))); for (j=1;j<=nyrs_ts;j++) { for(i=1;i<=3;i++) { t1 = pts(j,i)(1.0-p_ts(j,i))+0.1/3.0; t0 = log(mfexp(-1._square(p_ts(j,i)-P_ts(j,i))_effn_ts(j)/(2.0_t1))+0.01); // Like_ts(j,i) = (p_ts(j,i)-P_ts(j,i))_sqrt(effn_ts(j)/(2.0_t1)); LogLike(4) += (-0.5)_log(6.29*t1)-log(sqrt(1.0/effn_ts(j)))+t0; } }

// 5. Pot survey proportions are multinomial wrt predicted proportions // LogLike(5) = effn_ps_rowsum(elem_prod(p_ps,log(P_ps+0.01))); for (j=1;j<=nyrs_ps;j++) { for(i=1;i<=3;i++) { t1 = pps(j,i)(1.0-p_ps(j,i))+0.1/3.0; t0 = log(mfexp(-1._square(p_ps(j,i)-P_ps(j,i))_effn_ps(j)/(2.0_t1))+0.01); // Like_ps(j,i) = (p_ps(j,i)-P_ps(j,i))_sqrt(effn_ps(j)/(2.0_t1)); LogLike(5) += (-0.5)_log(6.29*t1)-log(sqrt(1.0/effn_ps(j)))+t0; } }

// 6. Observer proportions are multinomial wrt predicted proportions //LogLike(6) = effn_ob_rowsum(elem_prod(p_ob,log(P_ob+0.01))); for (j=1;j<=nyrs_ob;j++) { for(i=1;i<=3;i++) { t1 = pob(j,i)(1.0-p_ob(j,i))+0.1/3.0; t0 = log(mfexp(-1._square(p_ob(j,i)-P_ob(j,i))_effn_ob(j)/(2.0_t1))+0.01); // Like_ob(j,i) = (p_ob(j,i)-P_ob(j,i))_sqrt(effn_ob(j)/(2.0_t1)); LogLike(6) += (-0.5)_log(6.29_t1)-log(sqrt(1.0/effn_ob(j)))+t0; } } //7. Pot discard male biomass LogLike(7) = -0.5_norm2(log(x_ob1 + 0.001) - log(X_ob1 + 0.001)); LogLike(8) = -0.5*norm2(log(x_ob2 + 0.001) - log(X_ob2 + 0.001));

 // 9. + 10. Groundfish trawl and fixed-gear mortality biomass
 LogLike(9) = 0.0; LogLike(10) = 0.0;
 for(j=1;j<=nyrs_gf;j++)
 {
   LogLike(9) += -0.5*square(log(gft_mort(j)+0.01) - log(B_gft(yid_gf(j))+0.01));
   LogLike(10) += -0.5*square(log(gff_mort(j)+0.01) - log(B_gff(yid_gf(j))+0.01));
 }

 // Quadratic (normal) penalties

 // 1. Model recruit deviations
 Pen(1) = 0.5*norm2(log_New_dev);

 // 2. Directed fishery log fishing mortality deviations
 Pen(2) = 0.5*norm2(log_Fpf_dev);

 // 3. + 4. Gft and Gff log fishing mortality deviations
 Pen(3) = 0.5*norm2(log_Fgft_dev);
 Pen(4) = 0.5*norm2(log_Fgff_dev);

 // Objective function
 f = Pw*Pen - Lw*LogLike;