[ADMB Users] admb re gammln

[Anders Nielsen]

Hi, 

Sorry to respond this late, but I have been on vacation. 

Aparently the gammln is there, but not for the vector case, so 
a quick fix is to add a GLOBALS_SECTION like below: 

GLOBALS_SECTION
  #include <df1b2fun.h>

  df1b2vector gammln(const df1b2vector& z){
    int from=z.indexmin();
    int to=z.indexmax();
    df1b2vector ret(from,to); 
    for(int i=from; i<=to; ++i){
      ret(i)=gammln(z(i));
    }
    return(ret);
  }

A complete test case is attached. 

Unless there are big objections I will add the function to ADMB. 

Cheers, 

Anders. 




On Thu, 2010-08-05 at 10:40 -0700, Mark Maunder wrote:
> 
> Does anyone know if gammln has been implemented in admb-re?
-------------- next part --------------
GLOBALS_SECTION
  #include <df1b2fun.h>

  df1b2vector gammln(const df1b2vector& z){
    int from=z.indexmin();
    int to=z.indexmax();
    df1b2vector ret(from,to); 
    for(int i=from; i<=to; ++i){
      ret(i)=gammln(z(i));
    }
    return(ret);
  }

DATA_SECTION
  init_vector X(1,15);

PARAMETER_SECTION
  init_number logsize;
  init_bounded_number p(0,1);
  random_effects_vector dummy(1,1);
  sdreport_number size;
  objective_function_value jnll;

PROCEDURE_SECTION
  size=exp(logsize);
  int N=X.indexmax()-X.indexmin()+1;
  jnll=-sum(gammln(X+size))+N*gammln(size)+
       sum(gammln(X+1.0))-N*size*log(p)-sum(X)*log(1.0-p);

  jnll+=norm2(dummy)/2.0;
-------------- next part --------------
13 5 28 28 15 4 13 4 10 17 11 13 12 17 3