[ADMB Users] non-convergence and "inner maxg = 0"

[Ian Fiske]

Thanks Hans, Lennart, and Dave for your suggestions and advice.  After
poking around, I found that the problem only occurs if I initialize
the value of the objective function within the PROCEDURE_SECTION.
Specifically, if g is my objective function, mu_r is random, and y is
data, then

  g = -nroute*log_sigma_mu - 0.5*norm2(mu_r)/mfexp(2*log_sigma_mu);
  g -= -nobs*log_sigma_y - 0.5*norm2(y - Ey)/mfexp(2*log_sigma_y);

causes the inner maxg = 0 infinite loop, whereas

  g -= -nroute*log_sigma_mu - 0.5*norm2(mu_r)/mfexp(2*log_sigma_mu);
  g -= -nobs*log_sigma_y - 0.5*norm2(y - Ey)/mfexp(2*log_sigma_y);

appears to work (notice g is only decremented, not initialized).  Some
examples initialize the objective function (logistic in admb-re
examples) and others do not (e.g., buscycle in admb examples).  Which
is correct?

Dave, I've posted my tpl file is below as you generously requested.
Be warned, it's ugly right now -- I'm working on replacing the hideous
indexing with ragged 3d arrays.

Thanks,
Ian

DATA_SECTION

   init_int nobs                               // Number of data points
   init_int nroute                             // Number of routes
   init_int nrouteyr                           // Number of unique
route*year combinations
   init_int ncovs                              // Number of predictors
for each parameters
   init_int nparams                            // Number of parameters
in the seasonal model
   init_vector y(1,nobs)                       // Response vector
   init_vector jd(1,nobs)                      // Julian day
   init_vector route(1,nobs)                   // Vector of routes
   init_vector routeyr(1,nobs)
   init_vector routeyr_route(1,nrouteyr)
   init_matrix X(1,nrouteyr,1,ncovs)           // Design matrix of covariates

PARAMETER_SECTION

   init_bounded_vector mu_pars(1,ncovs,-10,10,1)
   init_bounded_vector beta_pars(1,ncovs,-10,10,1)
   init_bounded_vector phi_pars(1,ncovs,-10,10,1)
   init_bounded_number log_sigma_y(-5,5,1)
   init_bounded_number log_sigma_mu(-5,5,2)

   random_effects_vector mu_r(1,nroute,2)

   vector Emu_siteyr(1,nrouteyr)
   vector Ebeta_siteyr(1,nrouteyr)
   vector Ephi_siteyr(1,nrouteyr)
   vector mu(1,nrouteyr)
   vector beta(1,nrouteyr)
   vector phi(1,nrouteyr)
   vector Ey(1,nobs)

   objective_function_value g

GLOBALS_SECTION

  #include <df1b2fun.h>

PROCEDURE_SECTION

   Emu_siteyr = X * mu_pars;
   Ebeta_siteyr = X * beta_pars;
   Ephi_siteyr = X * phi_pars;

   // Create route-yr level seasonal curve parameters
   for (int i = 1; i <= nrouteyr; i++) {
     mu(i) = Emu_siteyr(i) + mu_r(routeyr_route(i));
     beta(i) = mfexp(Ebeta_siteyr(i));
     phi(i) = mfexp(Ephi_siteyr(i));
   }

   // Compute expected value for an observation
   for (int i = 1; i <= nobs; i++) {
     dvariable z = (jd(i) - mu(routeyr(i)))/beta(routeyr(i));
     Ey(i) = phi(routeyr(i)) * mfexp(-z + 1) * mfexp(-mfexp(-z));
   }

   g = -nroute*log_sigma_mu - 0.5*norm2(mu_r)/mfexp(2*log_sigma_mu);
   g -= -nobs*log_sigma_y - 0.5*norm2(y - Ey)/mfexp(2*log_sigma_y);

REPORT_SECTION

  report << mu_pars << endl;
  report << beta_pars << endl;
  report << phi_pars << endl;

TOP_OF_MAIN_SECTION

  arrmblsize = 40000000L;
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(3000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(200000);
  gradient_structure::set_MAX_NVAR_OFFSET(60000);

RUNTIME_SECTION

  maximum_function_evaluations 500
  convergence_criteria 1.e-3