// ADMB-RE implementation of from Zeh and Punt (2004) SC/56/BRG1 DATA_SECTION init_int n // Number of years with data init_matrix D(1,n,1,7) // year, N4, SE(N4), indicator, P4, SE(N4), indicator vector cv_N4(1,n) // CV of N4 vector cv_P4(1,n) // CV of P4 ivector indN4(1,n) // Indicator for whether there are N4 data in a particular year ivector indP4(1,n) // Indicator for whether there are P4 data in a particular year PARAMETER_SECTION init_bounded_number log_s_eta(-10,10) // log(sigma_eta) vector log_p(1,n) // log(p) sdreport_vector Py(1,n) // Abundance random_effects_vector log_pi(1,1) // log(pi) random_effects_vector a(0,1) // Intercept and Rate of increase (alternative model) random_effects_vector eta(1,n) objective_function_value g PRELIMINARY_CALCS_SECTION cout << setprecision(4); // Punt's cv's cv_N4 = sqrt(log(square(elem_div(column(D,3),column(D,2)))+1.0)); cv_P4 = sqrt(log(square(elem_div(column(D,6),column(D,5)))+1.0)); indN4 = (ivector) column(D,4); indP4 = (ivector) column(D,7); PROCEDURE_SECTION int i; // Contribution from random effects distribution g = -n*log_s_eta - .5*norm2(eta)/mfexp(2*log_s_eta); log_p = log_pi(1) + eta; for(i=1;i<=n;i++) { if(indN4(i)) g += -0.5*square((log(D(i,2)) - log_p(i) - a(0) - a(1)*(D(i,1)-1980))/cv_N4(i)); // Rate of increase verison if(indP4(i)) g += -0.5*square((log(D(i,5))-log_p(i))/cv_P4(i)); } g *= -1.0; // ADMB does minimization !!! Py = exp(a(0) + a(1)*(column(D,1)-1980.0)); REPORT_SECTION report << "s_eta = " << exp(log_s_eta) << endl; report << "pi = " << exp(log_pi) << endl; report << "Py(1978-1983,1985-1988,1993,2001) = " << Py << endl; TOP_OF_MAIN_SECTION arrmblsize = 4000000L; gradient_structure::set_GRADSTACK_BUFFER_SIZE(300000); gradient_structure::set_CMPDIF_BUFFER_SIZE(20000); gradient_structure::set_MAX_NVAR_OFFSET(50850);