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);