  /**
   * Apply inverse link function to linear predictor
   * see help(binom) in R
   *
   * @param eta Linear predictor vector
   * @param link An integer indicating the link function
   * @return A vector, i.e. inverse-link(eta)
   */
  vector linkinv_bern(vector eta, int link) {
    if (link == 1)      return(inv_logit(eta)); // logit
    else if (link == 2) return(Phi(eta)); // probit
    else if (link == 3) return(atan(eta) / pi() + 0.5);  // cauchit
    else if (link == 4) return(exp(eta)); // log
    else if (link == 5) return(inv_cloglog(eta)); // cloglog
    else reject("Invalid link");
    return eta; // never reached
  }

  /**
   * Increment with the unweighted log-likelihood
   * @param link An integer indicating the link function
   * @param eta0 A vector of linear predictors | y = 0
   * @param eta1 A vector of linear predictors | y = 1
   * @param N An integer array of length 2 giving the number of
   *   observations where y = 0 and y = 1 respectively
   * @return lp__
   */
  real bern_lpdf(vector eta0, vector eta1, int link, array[] int N) {
    real lp = 0;
    if (link == 1) { // logit
      lp += logistic_lccdf(eta0 | 0, 1);
      lp += logistic_lcdf( eta1 | 0, 1);
    }
    else if (link == 2) {  // probit
      lp += normal_lccdf(eta0 | 0, 1);
      lp += normal_lcdf( eta1 | 0, 1);
    }
    else if (link == 3) {  // cauchit
      lp += cauchy_lccdf(eta0 | 0, 1);
      lp += cauchy_lcdf( eta1 | 0, 1);
    }
    else if(link == 4) {  // log
      lp += sum(log1m_exp(eta0));
      lp += sum(eta1);  // already in log form
    }
    else if(link == 5) {  // cloglog
      lp += sum(log1m_exp(-exp(eta1)));
      lp += sum(-exp(eta0));
    }
    else reject("Invalid link");
    return lp;
  }

  /**
   * Pointwise (pw) log-likelihood vector
   *
   * @param y The integer outcome variable. Note that function is
   *  called separately with y = 0 and y = 1
   * @param eta Vector of linear predictions
   * @param link An integer indicating the link function
   * @return A vector
   */
  vector pw_bern(int y, vector eta, int link) {
    int N = rows(eta);
    vector[N] ll;
    if (link == 1) {  // logit
      for (n in 1:N) ll[n] = bernoulli_logit_lpmf(y | eta[n]);
    }
    else if (link <= 5) {  // link = probit, cauchit, log, or cloglog
      vector[N] pi = linkinv_bern(eta, link); // may not be stable
      for (n in 1:N) ll[n] = bernoulli_lpmf(y | pi[n]);
    }
    else reject("Invalid link");
    return ll;
  }

  /**
   * Log-normalizing constant in the clogit case
   *
   * @param N_j Integer number of observations in the j-th group
   * @param D_j Integer number of successes in the j-th group
   * @param eta_j Vector of linear predictions in the j-th group
   * @return A scalar that normalizes the probabilities on the log-scale
   */
  real log_clogit_denom(int N_j, int D_j, vector eta_j) {
    if (D_j == 1 && N_j == rows(eta_j)) return log_sum_exp(eta_j);
    if (D_j == 0) return 0;
    if (N_j == D_j) {
      if (D_j == 1) return eta_j[N_j];
      return sum(segment(eta_j, N_j - 1, 2));
    }
    else {
      int N_jm1 = N_j - 1;
      return log_sum_exp(log_clogit_denom(N_jm1, D_j, eta_j),
                         log_clogit_denom(N_jm1, D_j - 1, eta_j) + eta_j[N_j]);
    }
    return not_a_number();  // never reaches
  }

  /**
   * Log-likelihood for a clogit model
   * @param eta0 Linear predictors when y == 0
   * @param eta1 Linear predictors when y == 1
   * @param successes Integer array with the number of successes in group j
   * @param failures Integer array with the number of failures in group j
   * @param observations Integer array with the number of observations in group j
   * @return lp__
   */
  real clogit_lpdf(vector eta0, vector eta1,
                    array[] int successes, array[] int failures, array[] int observations) {
    int J = num_elements(observations);
    int pos0 = 1;
    int pos1 = 1;
    vector[J] summands;
    for (j in 1:J) {
      int D_g = successes[j];
      int N_g = observations[j];
      int F_g = failures[j];
      vector[N_g] eta_g = append_row(segment(eta1, pos1, D_g),
                                     segment(eta0, pos0, F_g));
      summands[j] = log_clogit_denom(N_g, D_g, eta_g);
      pos0 += F_g;
      pos1 += D_g;
    }
    return sum(eta1) - sum(summands);
  }
