[ADMB Users] Beta deviate

[Coilin Minto]

Dear ADMB Help List,

I would like to implement beta-distributed random effects via the
probability integral transform in the "beta_deviate" function, as in
section 3.4 "Non-Gaussian random effects" of the ADMB-RE manual.

adlink fails with the following error message:

  df1b2invcumdbeta_incbet.cpp:58: undefined reference to
`incbet(df3_three_variable const&, df3_three_variable const&, double
const&)'

Line 58 of df1b2invcumdbeta_incbet.cpp needs to implement the
incomplete beta integral for a double x. There is a reference to:

  df3_three_variable incbet(const df3_three_variable& a, const
df3_three_variable& b,const  double& x);

on line 35 but I can't find a defined function for this in the source.

Below is a nonsensical (given estimated values of gamma-distributed
random effects) but hopefully illustrative example of the problem
based on the 'gamma' example (examples/admb-re/liver_gamma.tpl).

I am running admb-11 with gcc 4.4.3 on Ubuntu 10.04.

Many thanks for your time and any assistance.

Best regards,
Coilin

//----------------
// liver_beta.tpl
//----------------

DATA_SECTION
  init_int np
  init_int nh
  init_matrix data(1,np,1,4)
  ivector I(1,np)
  ivector nump(1,nh)
  vector S(1,np)
  vector TRT(1,np)
  vector CARD(1,np)
 LOC_CALCS
  I=ivector(column(data,1));
  S=column(data,2);
  TRT=column(data,3);
  CARD=column(data,4);
  int ic=1;
  for (int i=1;i<np;i++)
  {
    if (I(i+1)>I(i))
    {
      nump(I(i))=ic;
      ic=1;
    }
    else
    {
      ic++;
    }
  }
  nump(nh)=ic;

PARAMETER_SECTION
  init_vector beta(0,2);
  init_bounded_number log_theta1(-5.0,3.0,2)
  init_bounded_number log_theta2(-5.0,3.0,2)

  random_effects_vector u(1,nh,2)
  objective_function_value f

PROCEDURE_SECTION
  int i;
  int j=0;
  for (i=1;i<=nh;i++)
  {
    fun(i,j,u(i),log_theta1,log_theta2,beta);
    //fun(i,j,u(i),log_theta1,beta);
  }

SEPARABLE_FUNCTION void fun( int i,int & j ,const prevariable& ui,
const prevariable& log_theta1, const prevariable& log_theta2, const
dvar_vector& beta)
//SEPARABLE_FUNCTION void fun( int i,int & j ,const prevariable& ui,
const prevariable& log_theta1, const dvar_vector& beta)

 f += 0.9189385 + 0.5*square(ui);  // N(0,1) likelihood contribution from u's

 // Likelihood contribution
 dvariable theta1=exp(log_theta1);
 dvariable theta2=exp(log_theta2);
 int ii;
 //dvariable z=cumd_norm(ui);         // z has uniform (0,1) distribution
 //z = 0.99999999*z + 0.000000005;      // To gain numerical stability
 //dvariable gi = theta1*inv_cumd_gamma(z,1.0/theta1);
 dvariable gi = beta_deviate(ui,theta1,theta2);
 for (ii=1;ii<=nump(i);ii++)
 {
   j++;

   dvariable log_lambda=beta(0)+beta(1)*TRT(j)+beta(2)*CARD(j)+log(gi);
   dvariable lambda=mfexp(log_lambda);
   f += lambda*S(j) - log_lambda;
 }


//----------------
// liver_beta.dat
//----------------

#INST SURVT treat heart
# number of patients  number of hospitals
    21                    7
#
1    23.286    1    1
1    6.429    0    0
1    26.857    1    0
1    11.143    0    1
2    3.857    0    1
2    9.000    1    0
2    8.714    0    0
3    1.143    0    1
3    23.143    1    0
3    2.571    1    0
3    2.857    0    1
4    76.429    1    0
4    35.857    1    1
4    25.857    0    0
4    52.286    0    0
5    25.143    1    1
5    25.143    0    0
6    29.286    1    1
6    28.857    0    0
7    1.857    0    1
7    14.286    1    0