[ADMB Users] big(ger) models

[dave fournier]

I noticed that someone had posted  on the SAS list about a moderate 
sized model
requiring something like 80 hours to fit using GLIMMIX.  I thought this 
might be
an interesting application of some of the sparse matrix stuff.  The 
poster had
posted a simple example of the kind of problem he was having trouble with.

http://listserv.uga.edu/cgi-bin/wa?A2=ind1405b&L=sas-l&P=4001

He also later said that his real problem has an ar(1) structure which we 
know how to deal with.
(I think).  I wrote the attached TPL file to test the ideas.

I ran the model with the following.

time ./ryan -mno 1000000 -ams 10000000 -ilmn 5 -shess -ndi 100000 
-noinit -l1 100000000 -nl1 100000000  -iprint 1  -crit 1.e-3 -noreders

the -noreders  is to prevent the calculation of the estimates of the std 
devs of the random effects via the delta method.

This saves a lot of time and space since there about 75,000 random 
effects.  Fitting the model required less than 4 minutes.
  I don't know how many ADMB people are interested in this kind of 
larger model.

It appears that the -noreders option has never been included in the 
code.  It is only about 8 or so lines of code to add.
Maybe someone is interested.  Per usual I will hold my breath.
-------------- next part --------------
DATA_SECTION
  int nloc
  int nsub
  int nobs
  int iseed
 LOC_CALCS
  nloc=250;
  nsub=300;
  nobs=4;
 END_CALCS
  3darray y(1,nloc,1,nsub,1,nobs)
  matrix x(1,nloc,1,nsub)
 LOC_CALCS
  iseed=12345;
  random_number_generator rng(iseed);

  double lre;
  double sre;
  double prob;
  double eta;

  for (int i=1;i<=nloc;i++)
  {
    lre=sqrt(.3)*randn(rng);
    for (int j=1;j<=nsub;j++)
    {
      sre=sqrt(.5)*randn(rng);
      x(i,j)=randn(rng);
      eta=1.6+2.5*x(i,j)+lre+sre;
      prob=exp(eta)/(1.0+exp(eta));
      for (int k=1;k<=nobs;k++)
      {
        double z=randu(rng);
        if (z<=prob)
          y(i,j,k)=1.0;
        else
          y(i,j,k)=0.0;
      }
    }
  }
        

PARAMETER_SECTION
  objective_function_value f
  init_number a
  init_number b
  init_bounded_number lsd(0.0,5.0,2)
  init_bounded_number ssd(0.0,5.0,2)
 !! lsd.set_scalefactor(100.);
 !! ssd.set_scalefactor(100.);
  random_effects_vector eloc(1,nloc,2)
  random_effects_matrix esub(1,nloc,1,nsub,2)

PROCEDURE_SECTION

  for (int i=1;i<=nloc;i++)
  {
    for (int j=1;j<=nsub;j++)
    {
      f1(i,j,a,b,lsd,ssd,eloc(i),esub(i,j));
    }
  }
       
SEPARABLE_FUNCTION void f1(int i,int j,const prevariable& a,const prevariable& b,const prevariable& lsd,const prevariable& ssd,const prevariable& el,const prevariable& es)

  dvariable eta=a+b*x(i,j)+lsd*el+ssd*es;
  dvariable prob;
  if (value(prob)<20.)
    prob=exp(eta)/(1.0+exp(eta));
  else
    prob=1.0/(1.0+exp(-eta));

  for (int k=1;k<=nobs;k++)
  {
    if (y(i,j,k)>0)
      f-=log(1.e-10+prob);
    else
      f-=log(1.e-10+(1.0-prob));
  }
  if (j==1)
  {
    f+=0.5*square(el);
  }
  f+=0.5*square(es);