Skip to content

MV tweedie and D-M for sizefreq #357

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rick-Methot-NOAA merged 18 commits into from
Aug 24, 2022
Merged

MV tweedie and D-M for sizefreq #357

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
28d100d
first commit for MV-Tweedie-ver2
Rick-Methot-NOAA Jul 19, 2022
e74694b
working version of revised D-M logic and parm_names
Rick-Methot-NOAA Jul 21, 2022
c41960a
Merge branch 'main' into MV-Tweedie-ver2
Rick-Methot-NOAA Jul 21, 2022
3fc92e7
create comp_logL_fxns; add D-M to sizefreq; fix diri2 error in bootstrap
Rick-Methot-NOAA Aug 2, 2022
0b07834
working; still need prelim offset for DM-sizefreq
Rick-Methot-NOAA Aug 3, 2022
552c996
D-M working for SizeFreq
Rick-Methot-NOAA Aug 4, 2022
b5fcd1c
just one small bug left in fit_size_summary
Rick-Methot-NOAA Aug 5, 2022
e6abe95
additional cleanup to Fit_Size_Summary table
Rick-Methot-NOAA Aug 5, 2022
628db9b
responding to Ian's review
Rick-Methot-NOAA Aug 5, 2022
ea7a8d8
fix to sizefreq in ssnew; augment report.sso
Rick-Methot-NOAA Aug 8, 2022
3faf8df
Merge branch 'main' into MV-Tweedie-ver2
Rick-Methot-NOAA Aug 11, 2022
d1d776a
Create SS_readcontrol_330.tpl.orig
Rick-Methot-NOAA Aug 11, 2022
8f08cd8
update warning text
Rick-Methot-NOAA Aug 11, 2022
59960a0
convert 2 write_warnings
Rick-Methot-NOAA Aug 11, 2022
f075a8f
begin streamline code for compreport.sso
Rick-Methot-NOAA Aug 12, 2022
a5ceb77
finish compreport update; verify results
Rick-Methot-NOAA Aug 12, 2022
61d2cba
modify associated text in ss_new
Rick-Methot-NOAA Aug 19, 2022
70a0d2e
minor wording change to trigger github actions
iantaylor-NOAA Aug 23, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/call-run-ss3-no-est.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# Build SS3 and run models without estimation and no hessian
# Build SS3 and run models without estimation (from .par) and no hessian
name: call-run-ss3-no-est
on:
workflow_dispatch:
Expand Down
17 changes: 17 additions & 0 deletions SS_miscfxn.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,3 +103,20 @@ FUNCTION void get_catch_mult(int y, int catch_mult_pointer)
return;
}

//********************************************************************
/* SS_Label_FUNCTION 4XX Comp_logL */
FUNCTION dvariable Comp_logL_multinomial(const double& Nsamp, const dvector& obs_comp, const dvar_vector& exp_comp)
{
dvariable logL;
// logL = - Nsamp * obs_comp(tail_L, tail_H) * log(exp_comp(tail_L, tail_H));
// the call to this function does the subsetting to tail_L and tail_H, so this function can operate cleanly on the entirety of the passed vector
logL = - Nsamp * obs_comp * log(exp_comp);
return (logL);
}

FUNCTION dvariable Comp_logL_Dirichlet(const double& Nsamp, const dvariable& dirichlet_Parm, const dvector& obs_comp, const dvar_vector& exp_comp)
{
dvariable logL;
logL = sum(gammln(Nsamp * obs_comp + dirichlet_Parm * exp_comp)) - sum(gammln(dirichlet_Parm * exp_comp));
return (logL);
}
125 changes: 85 additions & 40 deletions SS_objfunc.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -341,44 +341,61 @@ FUNCTION void evaluate_the_objective_function()
if (Comp_Err_L(f) == 0) // multinomial
{
// get female or combined sex logL
// logL functions are at end of file SS_miscfxn.tpl
if (gen_l(f, i) != 2)
length_like(f, i) -= nsamp_l(f, i) *
obs_l(f, i)(tails_w(1), tails_w(2)) * log(exp_l(f, i)(tails_w(1), tails_w(2)));
length_like(f, i) += Comp_logL_multinomial( nsamp_l(f, i), obs_l(f, i)(tails_w(1), tails_w(2)), exp_l(f, i)(tails_w(1), tails_w(2)) );
// add male logL
if (gen_l(f, i) >= 2 && gender == 2)
length_like(f, i) -= nsamp_l(f, i) *
obs_l(f, i)(tails_w(3), tails_w(4)) * log(exp_l(f, i)(tails_w(3), tails_w(4)));
length_like(f, i) += Comp_logL_multinomial( nsamp_l(f, i), obs_l(f, i)(tails_w(3), tails_w(4)), exp_l(f, i)(tails_w(3), tails_w(4)) );
}
else // dirichlet
else if( (Comp_Err_L(f)==1) || (Comp_Err_L(f)==2) ) // dirichlet
{
/* from Thorson: NLL -= gammln(A) - gammln(ninput_t(t)+A) + sum(gammln(ninput_t(t)*extract_row(pobs_ta,t) + A*extract_row(pexp_ta,t))) - sum(lgamma(A*extract_row(pexp_ta,t))) \
dirichlet_Parm=mfexp(selparm(Comp_Err_Parm_Start+Comp_Err_L2(f)))*nsamp_l(f,i);
in option 1, dirichlet_Parm = Theta*n from equation (10) of Thorson et al. 2016
in option 2, dirichlet_Parm = Beta from equation (4) of Thorson et al. 2016 */
dirichlet_Parm = mfexp(selparm(Comp_Err_parmloc(Comp_Err_L2(f),1)));
if (Comp_Err_L(f) == 1)
dirichlet_Parm = mfexp(selparm(Comp_Err_Parm_Start + Comp_Err_L2(f))) * nsamp_l(f, i);
if (Comp_Err_L(f) == 2)
dirichlet_Parm = mfexp(selparm(Comp_Err_Parm_Start + Comp_Err_L2(f)));
// dirichlet_Parm=mfexp(selparm(Comp_Err_Parm_Start+Comp_Err_L2(f)));
dirichlet_Parm *= nsamp_l(f, i);

// note: first term in equations (4) and (10) is calculated
// as offset_l in SS_prelim.tpl and already included in length_like
// now add second term which is only dependent on parameters and sample size
temp = gammln(dirichlet_Parm) - gammln(nsamp_l(f, i) + dirichlet_Parm);
// get female or combined sex logL
// third and final term in equations (4) and (10)
if (gen_l(f, i) != 2) // so not male only
{
temp += sum(gammln(nsamp_l(f, i) * obs_l(f, i)(tails_w(1), tails_w(2)) + dirichlet_Parm * exp_l(f, i)(tails_w(1), tails_w(2))));
temp -= sum(gammln(dirichlet_Parm * exp_l(f, i)(tails_w(1), tails_w(2))));
if (gen_l(f, i) != 2) { // so not male only
temp += Comp_logL_Dirichlet( nsamp_l(f, i), dirichlet_Parm, obs_l(f, i)(tails_w(1), tails_w(2)), exp_l(f, i)(tails_w(1), tails_w(2)) );
}
// add male logL
if (gen_l(f, i) >= 2 && gender == 2)
{
temp += sum(gammln(nsamp_l(f, i) * obs_l(f, i)(tails_w(3), tails_w(4)) + dirichlet_Parm * exp_l(f, i)(tails_w(3), tails_w(4))));
temp -= sum(gammln(dirichlet_Parm * exp_l(f, i)(tails_w(3), tails_w(4))));
if (gen_l(f, i) >= 2 && gender == 2) {
temp += Comp_logL_Dirichlet( nsamp_l(f, i), dirichlet_Parm, obs_l(f, i)(tails_w(3), tails_w(4)), exp_l(f, i)(tails_w(3), tails_w(4)) );
}
length_like(f, i) -= temp;
} else // multivariate-Tweedie
{
dvariable tweedie_Phi;
dvariable tweedie_power;
// Exponentiate [PARAMETER_1]
int k1 = Comp_Err_parmloc(Comp_Err_L2(f),1);
tweedie_Phi = mfexp(selparm(k1));
// One plus logistic-transform [PARAMETER_1]
tweedie_power = 1.0 + mfexp(selparm(k1+1)) / (1.0+mfexp(selparm(k1+1)));
if(gen_l(f,i) !=2) // so not male only
{
// dtweedie( Type y, Type mu, Type phi, Type p, int give_log=0 )
for (int tail_index=tails_w(1); tail_index<=tails_w(2); tail_index++){
temp += 1.; // dtweedie( nsamp_l(f,i)*obs_l(f,i)(tail_index), nsamp_l(f,i)*exp_l(f,i)(tail_index), tweedie_Phi, tweedie_power, true );
}
}
// add male logL
if(gen_l(f,i) >=2 && gender==2)
{
// dtweedie( Type y, Type mu, Type phi, Type p, int give_log=0 )
for (int tail_index=tails_w(3); tail_index<=tails_w(4); tail_index++){
temp += 1.; // dtweedie( nsamp_l(f,i)*obs_l(f,i)(tail_index), nsamp_l(f,i)*exp_l(f,i)(tail_index), tweedie_Phi, tweedie_power, true );
}
}
length_like(f,i)-=temp;
}
if (header_l(f, i, 3) > 0)
length_like_tot(f) += length_like(f, i);
Expand Down Expand Up @@ -480,46 +497,43 @@ FUNCTION void evaluate_the_objective_function()
if (Comp_Err_A(f) == 0) // multinomial
{
if (gen_a(f, i) != 2)
age_like(f, i) -= nsamp_a(f, i) *
obs_a(f, i)(tails_w(1), tails_w(2)) * log(exp_a(f, i)(tails_w(1), tails_w(2)));
age_like(f, i) += Comp_logL_multinomial( nsamp_a(f, i), obs_a(f, i)(tails_w(1), tails_w(2)), exp_a(f, i)(tails_w(1), tails_w(2)) );
// age_like(f, i) -= nsamp_a(f, i) *
// obs_a(f, i)(tails_w(1), tails_w(2)) * log(exp_a(f, i)(tails_w(1), tails_w(2)));
if (gen_a(f, i) >= 2 && gender == 2)
age_like(f, i) -= nsamp_a(f, i) *
obs_a(f, i)(tails_w(3), tails_w(4)) * log(exp_a(f, i)(tails_w(3), tails_w(4)));
age_like(f, i) += Comp_logL_multinomial( nsamp_a(f, i), obs_a(f, i)(tails_w(3), tails_w(4)), exp_a(f, i)(tails_w(3), tails_w(4)) );
// age_like(f, i) -= nsamp_a(f, i) *
// obs_a(f, i)(tails_w(3), tails_w(4)) * log(exp_a(f, i)(tails_w(3), tails_w(4)));
}
else // dirichlet
else if( (Comp_Err_A(f)==1) || (Comp_Err_A(f)==2) ) // dirichlet
{
/* from Thorson: NLL -= gammln(A) - gammln(ninput_t(t)+A) + sum(gammln(ninput_t(t)*extract_row(pobs_ta,t) + A*extract_row(pexp_ta,t))) - sum(lgamma(A*extract_row(pexp_ta,t))) \
dirichlet_Parm=mfexp(selparm(Comp_Err_Parm_Start+Comp_Err_A2(f)))*nsamp_a(f,i);
in option 1, dirichlet_Parm = Theta*n from equation (10) of Thorson et al. 2016
in option 2, dirichlet_Parm = Beta from equation (4) of Thorson et al. 2016
*/
dirichlet_Parm = mfexp(selparm(Comp_Err_parmloc(Comp_Err_A2(f),1)));
if (Comp_Err_A(f) == 1)
dirichlet_Parm = mfexp(selparm(Comp_Err_Parm_Start + Comp_Err_A2(f))) * nsamp_a(f, i);
if (Comp_Err_A(f) == 2)
dirichlet_Parm = mfexp(selparm(Comp_Err_Parm_Start + Comp_Err_A2(f)));
// dirichlet_Parm=mfexp(selparm(Comp_Err_Parm_Start+Comp_Err_A2(f)));

dirichlet_Parm *= nsamp_a(f, i);
// note: first term in equations (4) and (10) is calculated
// as offset_a in SS_prelim.tpl and already included in age_like
// now add second term which is only dependent on parameters and sample size
// second term in equations (4) and (10) which is only dependent on parameters and sample size
temp = gammln(dirichlet_Parm) - gammln(nsamp_a(f, i) + dirichlet_Parm);
// get female or combined sex logL
// final term in equations (4) and (10)
if (gen_a(f, i) != 2) // so not male only
{
temp += sum(gammln(nsamp_a(f, i) * obs_a(f, i)(tails_w(1), tails_w(2)) + dirichlet_Parm * exp_a(f, i)(tails_w(1), tails_w(2))));
temp -= sum(gammln(dirichlet_Parm * exp_a(f, i)(tails_w(1), tails_w(2))));
if (gen_a(f, i) != 2) { // so not male only
temp += Comp_logL_Dirichlet( nsamp_a(f, i), dirichlet_Parm, obs_a(f, i)(tails_w(1), tails_w(2)), exp_a(f, i)(tails_w(1), tails_w(2)) );
}
// add male logL
if (gen_a(f, i) >= 2 && gender == 2)
{
temp += sum(gammln(nsamp_a(f, i) * obs_a(f, i)(tails_w(3), tails_w(4)) + dirichlet_Parm * exp_a(f, i)(tails_w(3), tails_w(4))));
temp -= sum(gammln(dirichlet_Parm * exp_a(f, i)(tails_w(3), tails_w(4))));
if (gen_a(f, i) >= 2 && gender == 2) {
temp += Comp_logL_Dirichlet( nsamp_a(f, i), dirichlet_Parm, obs_a(f, i)(tails_w(3), tails_w(4)), exp_a(f, i)(tails_w(3), tails_w(4)) );
}
age_like(f, i) -= temp;
}
}
else // MV_Tweedie
{
}
if (header_a(f, i, 3) > 0)
age_like_tot(f) += age_like(f, i);
}
Expand Down Expand Up @@ -564,6 +578,7 @@ FUNCTION void evaluate_the_objective_function()
}

// SS_Label_Info_25.7 #Fit to generalized Size composition
dvariable temp_logL;
if (SzFreq_Nmeth > 0) // have some sizefreq data
{
// create super-period expected values
Expand All @@ -581,7 +596,7 @@ FUNCTION void evaluate_the_objective_function()
} // assign back to all obs
}

SzFreq_like = -SzFreq_like_base; // initializes
SzFreq_like = -offset_Sz_tot; // initializes for each Sz_Method
for (iobs = 1; iobs <= SzFreq_totobs; iobs++)
{
if (SzFreq_obs_hdr(iobs, 3) > 0)
Expand All @@ -590,12 +605,42 @@ FUNCTION void evaluate_the_objective_function()
f = abs(SzFreq_obs_hdr(iobs, 3));
z1 = SzFreq_obs_hdr(iobs, 7);
z2 = SzFreq_obs_hdr(iobs, 8);
SzFreq_like(SzFreq_LikeComponent(f, k)) -= SzFreq_sampleN(iobs) * SzFreq_obs(iobs)(z1, z2) * log(SzFreq_exp(iobs)(z1, z2));
int Sz_method = SzFreq_obs1(iobs, 1); // sizecomp method
int logL_method = Comp_Err_Sz(Sz_method);
temp_logL = 0.0;
switch (logL_method)
{
case 0: // multinomial
{
temp_logL += Comp_logL_multinomial( SzFreq_sampleN(iobs), SzFreq_obs(iobs)(z1, z2), SzFreq_exp(iobs)(z1, z2));
break;
}
case 1: // dirichlet with theta*n
{
dirichlet_Parm = mfexp(selparm(Comp_Err_parmloc(Comp_Err_Sz2(Sz_method),1))) * SzFreq_sampleN(iobs); // theta * n
temp_logL -= gammln(dirichlet_Parm) - gammln( SzFreq_sampleN(iobs) + dirichlet_Parm );
temp_logL -= Comp_logL_Dirichlet( SzFreq_sampleN(iobs), dirichlet_Parm, SzFreq_obs(iobs)(z1, z2), SzFreq_exp(iobs)(z1, z2));
break;
}
case 2: // dirichlet with beta
{
dirichlet_Parm = mfexp(selparm(Comp_Err_parmloc(Comp_Err_Sz2(Sz_method),1))); // beta
temp_logL -= gammln(dirichlet_Parm) - gammln( SzFreq_sampleN(iobs) + dirichlet_Parm );
temp_logL -= Comp_logL_Dirichlet( SzFreq_sampleN(iobs), dirichlet_Parm, SzFreq_obs(iobs)(z1, z2), SzFreq_exp(iobs)(z1, z2));
break;
}
case 3: // MV Tweedie
{
break;
}
}
SzFreq_like(SzFreq_LikeComponent(f, k)) += temp_logL;
SzFreq_eachlike(iobs) = value(temp_logL) - SzFreq_each_offset(iobs);
}
}

if (do_once == 1)
cout << " did sizefreq obj_fun: " << SzFreq_like << " base: " << SzFreq_like_base << endl;
cout << " did sizefreq obj_fun: " << SzFreq_like << " base: " << offset_Sz_tot << endl;
}

// SS_Label_Info_25.8 #Fit to morph composition
Expand Down
47 changes: 33 additions & 14 deletions SS_prelim.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -265,26 +265,29 @@
obs_l(f, i)(tails_l(f, i, 3), tails_l(f, i, 4)) * log(obs_l(f, i)(tails_l(f, i, 3), tails_l(f, i, 4)));
}
}
else
else if( (Comp_Err_L(f)==1) || (Comp_Err_L(f)==2) ) // dirichlet
{
// Dirichlet-Multinomial (either 1 = linear, 2 = saturating)
// cannot use fxn Comp_Err_Dirichlet for this calc because only need the first part here
offset_l(f, i) = gammln(nsamp_l(f, i) + 1.);
if (gen_l(f, i) != 2)
{
int z1 = tails_l(f, i, 1);
int z2 = tails_l(f, i, 2);
offset_l(f, i) += gammln(nsamp_l(f, i) + 1.) -
offset_l(f, i) -= sum(gammln(1. + nsamp_l(f, i) * obs_l(f, i)(z1, z2)));
// sum(gammln(1. + nsamp_l(f,i)*obs_l(f,i)(tails_l(f,i,3),tails_l(f,i,4))));
sum(gammln(1. + nsamp_l(f, i) * obs_l(f, i)(z1, z2)));
}
if (gen_l(f, i) >= 2 && gender == 2)
{
int z1 = tails_l(f, i, 3);
int z2 = tails_l(f, i, 4);
offset_l(f, i) += gammln(nsamp_l(f, i) + 1.) -
// sum(gammln(1. + nsamp_l(f,i)*obs_l(f,i)(tails_l(f,i,3),tails_l(f,i,4))));
sum(gammln(1. + nsamp_l(f, i) * obs_l(f, i)(z1, z2)));
offset_l(f, i) -= sum(gammln(1. + nsamp_l(f, i) * obs_l(f, i)(z1, z2)));
}
}
else if( (Comp_Err_L(f)==3)) // MV Tweedie
{
// no MV Tweedie offset
}
}
// echoinput<<" length_comp offset: "<<offset_l<<endl;
echoinput << " length comp var adjust has been set-up " << endl;
Expand Down Expand Up @@ -515,24 +518,28 @@
obs_a(f, i)(tails_a(f, i, 3), tails_a(f, i, 4)) * log(obs_a(f, i)(tails_a(f, i, 3), tails_a(f, i, 4)));
}
}
else
else if( (Comp_Err_A(f)==1) || (Comp_Err_A(f)==2) ) // dirichlet
{
// Dirichlet-Multinomial (either 1 = linear, 2 = saturating)
offset_a(f, i) = gammln(nsamp_a(f, i) + 1.);
if (gen_a(f, i) != 2)
{
int z1 = tails_a(f, i, 1);
int z2 = tails_a(f, i, 2);
offset_a(f, i) += gammln(nsamp_a(f, i) + 1.) -
sum(gammln(1. + nsamp_a(f, i) * obs_a(f, i)(z1, z2)));
offset_a(f, i) -= sum(gammln(1. + nsamp_a(f, i) * obs_a(f, i)(z1, z2)));
}
if (gen_a(f, i) >= 2 && gender == 2)
{
int z1 = tails_a(f, i, 3);
int z2 = tails_a(f, i, 4);
offset_a(f, i) += gammln(nsamp_a(f, i) + 1.) -
sum(gammln(1. + nsamp_a(f, i) * obs_a(f, i)(z1, z2)));
offset_a(f, i) -= sum(gammln(1. + nsamp_a(f, i) * obs_a(f, i)(z1, z2)));
}
}
else if( (Comp_Err_A(f)==3) ) // MV Tweedie
{
// MV Tweedie has no offset, at least yet
}

}
// echoinput<<" agecomp offset "<<offset_a<<endl;
echoinput << " age comp var adjust has been set-up " << endl;
Expand Down Expand Up @@ -576,9 +583,21 @@
z1 = SzFreq_obs_hdr(iobs, 7);
z2 = SzFreq_obs_hdr(iobs, 8);
g = SzFreq_LikeComponent(f, k);
SzFreq_like_base(g) -= SzFreq_sampleN(iobs) * SzFreq_obs(iobs)(z1, z2) * log(SzFreq_obs(iobs)(z1, z2));
if (Comp_Err_Sz(k) == 0) // Multinomial
{
offset_Sz_tot(g) -= SzFreq_sampleN(iobs) * SzFreq_obs(iobs)(z1, z2) * log(SzFreq_obs(iobs)(z1, z2));
SzFreq_each_offset(iobs) -= SzFreq_sampleN(iobs) * SzFreq_obs(iobs)(z1, z2) * log(SzFreq_obs(iobs)(z1, z2));
}
else if (Comp_Err_Sz(k) == 1 || Comp_Err_Sz(k) == 2 ) // Dirichlet
{
offset_Sz_tot(g) += gammln(SzFreq_sampleN(iobs) + 1.) - sum(gammln(1. + SzFreq_sampleN(iobs) * SzFreq_obs(iobs)(z1, z2)));
SzFreq_each_offset(iobs) += gammln(SzFreq_sampleN(iobs) + 1.) - sum(gammln(1. + SzFreq_sampleN(iobs) * SzFreq_obs(iobs)(z1, z2)));
}
else if (Comp_Err_Sz(k) == 3) // MV Tweedie
{
// MV Tweedie not available
}
}

// identify super-period starts and stops
if (s < 0) // start/stop a super-period ALL observations must be continguous in the file
{
Expand All @@ -595,7 +614,7 @@
}
}
}
echoinput << " gen size comp var adjust has been set-up " << endl;
echoinput << " Sizefreq comp var adjust has been applied and offset calculated " << endl;

if (N_suprper_SzFreq > 0)
{
Expand Down
Loading