DATA_SECTION

  init_int n
  init_int pgr
  init_int ptr
  init_vector sumh(1,n)
  init_vector tau(1,n)
  init_ivector k(1,n)
  init_matrix Zgr(1,n,1,pgr)
  init_matrix Ztr(1,n,1,ptr)
  init_imatrix Y(1,n,1,k)

PARAMETER_SECTION

  init_bounded_vector ctr(1,ptr,-1,1)
  init_bounded_vector cgr(1,pgr,-100,100) 
  vector pop_vec(1,n)
  vector betavec(1,n)
  vector survivorvec(1,n)
  vector incidence_vec(1,n)
  objective_function_value f

PROCEDURE_SECTION

  dvariable fpen = 0.0; //penalty variable
  incidence_vec = Zgr * cgr;
  betavec = Ztr * ctr;
  pop_vec = elem_div(incidence_vec, betavec);
  for (int i=1; i <=n; i++)
  {
    survivorvec(i) = pop_vec(i) - k(i);
    survivorvec(i) = posfun(survivorvec(i),0.001,fpen);
    betavec(i) = posfun(betavec(i),0.000001,fpen);
  }
  // negative log-likelihood
  f = 0.0;
  f += 1000 * fpen;
  for (int i=1; i<=n; i++)
  {
    f -= gammln(survivorvec(i) + k(i) + 1) - gammln(survivorvec(i) + 1);
    if (k(i) >= 1)
    {
      f -= -betavec(i) * sumh(i);
      for (int j=1; j <= k(i); j++)
      {
        f -= log(betavec(i) * Y(i,j));
      }
    }
    f -= survivorvec(i) * -betavec(i) * tau(i);
  }
