[ADMB Users] SEPARABLE_FUNCTION

[Richard Chandler]

Thank you Hans and Dave. The tpl file below now works. Weihai Liu was
also kind enough to send me a suggestion, which is pasted in too.
Since Hans is one of the principal developers, I assume his
implementation is the preferred one.

Richard


// Based on Hans' suggestions

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
  random_effects_vector u(1,nG,2)

  objective_function_value nll

PROCEDURE_SECTION

  for(int id=1; id<=nG; id++) {
      nll_group(id, p0,p1,log_lambda,log_sigma,u(id));
      }


SEPARABLE_FUNCTION void nll_group(int id, const dvariable& p0,const
dvariable& p1,const dvariable& log_lambda, const dvariable& log_sigma,
const dvariable& u)
  dvar_vector p(1,R);
  dvar_vector f(1,S);
  dvar_vector g(1,S);
  dvariable sigma = exp(log_sigma);
  dvariable lambda = exp(log_lambda);

  nll += log_sigma + 0.5*square(u/sigma);

  int q=10*(id-1)+1;
  for(int i=q;i<=q+9;i++) {
      p(i) = 1.0/(1.0+exp(-1.0*(p0 + p1*x(i) + u)));
      for(int k=1;k<=S;k++) {
          f(k) = pow(lambda, N(k)) * exp(-1.0*lambda - gammln(N(k)+1));
          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(f*g);
  }






// based on Weihai's example

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
  random_effects_vector u(1,nG,2)

  objective_function_value nll


PROCEDURE_SECTION
  nll_group(p0,p1,log_lambda,log_sigma,u);



SEPARABLE_FUNCTION void nll_group(const dvariable& p0,const dvariable&
p1,const dvariable& log_lambda, const dvariable& log_sigma, const
dvar_vector& u)
  dvar_vector p(1,R);
  dvar_vector f(1,S);
  dvar_vector g(1,S);
  dvariable sigma = exp(log_sigma);
  dvariable lambda = exp(log_lambda);

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


  for(int i=1;i<=R;i++) {
      int 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++) {
          f(k) = pow(lambda, N(k)) * exp(-1.0*lambda -
gammln(N(k)+1)); //poisson
          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));//binomial
              else
                  g(k) =0.;
              }
          }
      nll -= log(f*g + 1.e-20);
  }






On Fri, Mar 11, 2011 at 6:08 AM, H. Skaug <hskaug at gmail.com> wrote:
> Hi,
>
> Yes, it seems that your model is separable, but of type "4.1 The first example"
> in ADMBRE manual. Hence, you do not need the sparse matrix stuff.
> In your code you need to loop over the elements of "u" like in "4.1
> The first example"
> Your PROCEDURE_SECTION should basically consist of:
>
>  for(int id=1;id<=xxxxx;id++)
>   SEPARABLE_FUNCTION(u(id), + the rest of your arguments)
>
>
> Inside the SEPARABLE_FUNCTION you need to do all the calculations
> that involve u(id).
>
> If you do this you will be able to handle large models. Make sure
> that you get the same result on small examples with and with SEPARABLE_FUNCTION.
>
>
> Hans
>
> On Thu, Mar 10, 2011 at 2:54 PM, Richard Chandler
> <richard.chandlers at gmail.com> wrote:
>> 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
>>
>