DATA_SECTION

  init_int cgr_phase
  init_int ctr_phase
  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_ivector maxk(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,-10,10,ctr_phase)
  init_bounded_vector cgr(1,pgr,-10,10,cgr_phase) 
  vector survivorvec(1,n)
  vector betavec(1,n)
  vector extras(1,n)  // people not infected beyond the minimum required by the model
  objective_function_value f

PROCEDURE_SECTION

  f = 0;
  extras = mfexp(Zgr * cgr);
  betavec = mfexp(Ztr * ctr);
  for (int i=1; i <=n; i++)
  {
    survivorvec(i) = extras(i) + maxk(i) - k(i);
  }
  // negative log-likelihood
  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);
  }
