[ADMB Users] SEPARABLE_FUNCTION

[Richard Chandler]

Hi,

I implemented a binomial-Poisson mixture model (Royle Biometrics 2004)
with random effects in ADMB (.tpl below). It works very well on small
simulated datasets; however, I run into memory issues when I try
complex variations of the model on large datasets. I have fiddled with
many of the command line options and the TOP_OF_MAIN_SECTION, but I
believe I need to use ADMB's sparse matrix capabilities via the
SEPARABLE_FUNCTION. My approach has been to start with something
simple:

SEPARABLE_FUNCTION void nll_group(const dvariable& log_sigma, const
dvar_vector& u)
  nll = nG*log_sigma + 0.5*norm2(u)/mfexp(2*log_sigma);

but this does not converge. I am new to ADMB and would greatly
appreciate it if someone could suggest a better way to use a
SEPARABLE_FUNCTION here. The RE manual is very helpful, but I must be
missing something. Also, I would appreciate any other suggestions for
improving the implementation of this model.

Thanks for the great software.

Richard


Details: 32-bit ADMB 10.0 on 64-bit Windows 7, MinGW




// Binomial-Poisson mixture model (Royle Biometrics 2004)
// Site-level covariate of p
// Random group effect

DATA_SECTION

  init_int R              // Number of sites
  init_int T              // Number of occasions
  init_int S              // Number of possible values of N
  init_int nG             // Number of groups

  init_vector N(1,S)      // Possible values of N
  init_vector ID(1,R)     // Group ID
  init_vector x(1,R)      // Site-specific covariate
  init_matrix y(1,R,1,T)  // Vector. Count data with R rows and T columns


PARAMETER_SECTION

  init_bounded_number log_lambda(-20.0,20.0,1)
  init_bounded_number p0(-20.0,20.0,1)
  init_bounded_number p1(-20.0,20.0,1)

  init_bounded_number log_sigma(-3.0,3.0,2) // log of random effect SD
  sdreport_number sigma

  number lambda
  vector p(1,R)

  vector f(1,S)
  vector g(1,S)

  random_effects_vector u(1,nG,2)

  objective_function_value nll

GLOBALS_SECTION
  #include <df1b2fun.h>


PROCEDURE_SECTION

  sigma = exp(log_sigma);

  nll = nG*log_sigma + 0.5*norm2(u)/mfexp(2*log_sigma);

  lambda = exp(log_lambda);

  for(int k=1;k<=S;k++) {
    f(k) = pow(lambda, N(k)) * exp(-1.0*lambda - gammln(N(k)+1));
    }

  int id;
  for(int i=1;i<=R;i++) {
      id = ID(i);
      p(i) = 1.0/(1.0+exp(-1.0*(p0 + p1*x(i) + u(id))));
      for(int k=1;k<=S;k++) {
          g(k) = 1.0;
          for(int j=1;j<=T;j++) {
              if(N(k)>=y(i,j))
                  g(k) *= exp(gammln(N(k)+1) - gammln(y(i,j)+1) -
                      gammln(N(k)-y(i,j)+1)) * pow(p(i), y(i,j)) *
                      pow(1.0-p(i), N(k)-y(i,j));
              else
                  g(k) = 0.0;
              }
          }
      nll -= log(sum(elem_prod(f, g)));
      }

REPORT_SECTION

  report << "lambda =  "  << lambda  << endl;
  report << "p0 =  "  << p0  << endl;
  report << "p1 =  "  << p1  << endl;
  report << "sigma = "  << sigma  << endl;
  report << "u = "  << u  << endl;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: nmix.dat
Type: application/octet-stream
Size: 6940 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20110310/e2e23751/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: nmix.tpl
Type: application/octet-stream
Size: 2073 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20110310/e2e23751/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: nmix.pin
Type: application/octet-stream
Size: 351 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20110310/e2e23751/attachment-0002.obj>