[ADMB Users] crappy software

[dave fournier]

Prabably a good idea to make this a bit more bulletproof
Bound rho in the choleski decomp a bit

   init_bounded_number rho(-3.,3.,3)


and add 1.e-12 to avoid log of zero
   {
     if (Y(i,t)==1)
       ll-=log(1.e-12+mu(t));
     else
       ll-=log(1.e-12+(1.0-mu(t)));
   }

which gives.



DATA_SECTION
   init_int n
   init_matrix x(1,n,1,3)
   init_matrix Y(1,n,1,3)
PARAMETER_SECTION
   init_number log_tau
   init_bounded_number rho(-3.,3.,3)
   matrix ch(1,2,1,2)
   init_vector beta(1,2)
   matrix b(1,n,1,2)
   sdreport_matrix sigma(1,2,1,2)
   sdreport_number tau
   random_effects_matrix bb(1,n,1,2,2)
   objective_function_value ll;
PROCEDURE_SECTION
   for (int i=1;i<=n;i++)
   {
     f1(i,log_tau,rho,bb(i),beta);
   }
   if (sd_phase())
   {
     tau=exp(log_tau);
     dvar_matrix ch(1,2,1,2);
     ch(1,1)=1.0;
     ch(1,2)=0.0;
     ch(2,1)=rho;
     ch(2,2)=1.0;
     ch(2)/=norm(ch(2));
     ch*=tau;
     sigma=ch*trans(ch);
   }

SEPARABLE_FUNCTION void f1(int& i,const prevariable & log_tau,const 
prevariable & rho,const dvar_vector & bbi,const dvar_vector & beta)
   dvariable tau=exp(log_tau);
   dvar_matrix ch(1,2,1,2);
   ch(1,1)=1.0;
   ch(1,2)=0.0;
   ch(2,1)=rho;
   ch(2,2)=1.0;
   ch(2)/=norm(ch(2));
   ch*=tau;
   dvar_vector b=ch*bbi;
   dvar_vector tmp=beta(1)+b(1)
     +(beta(2)+b(2))*x(i);
   dvar_vector mu=1.0/(1.0+exp(tmp));
   for (int t=1;t<=3;t++)
   {
     if (Y(i,t)==1)
       ll-=log(1.e-12+mu(t));
     else
       ll-=log(1.e-12+(1.0-mu(t)));
   }
   ll+=0.5*norm2(bbi);
                       
                                                               
19,1          54%