delay_multinomial_lpmf.stan

Go to the documentation of this file.

 
real delay_multinomial_lpmf(array[] int obs, vector log_exp_obs) {
  return multinomial_logit_lpmf(obs | log_softmax(log_exp_obs));
}
 
real delay_multinomial_snaps(int start, int end, array[] int obs,
                             vector log_exp_obs, int exp_offset,
                             array[] int total, array[] int obs_lookup,
                             array[] int lsl, array[] int clsl,
                             array[] int nsl, array[] int cnsl) {
  real tar = 0;
  for (i in start:end) {
    if (lsl[i]) {
      // Cutoff block of expected cells; renormalise over all slots up to the
      // cutoff (log_softmax) so before-cutoff interior cells keep weight.
      array[3] int c = filt_obs_indexes(i, i, clsl, lsl);
      vector[lsl[i]] lprob = log_softmax(
        segment(log_exp_obs, c[1] - exp_offset + 1, lsl[i])
      );
      array[3] int o = filt_obs_indexes(i, i, cnsl, nsl);
      if (nsl[i] == lsl[i]) {
        // Fully observed up to the cutoff: plain (truncated) multinomial.
        tar += multinomial_lpmf(segment(obs, o[1], nsl[i]) | exp(lprob));
      } else {
        // Some before-cutoff cells unobserved: marginalise them into one
        // residual category (count = total - observed, prob = 1 - observed).
        vector[nsl[i] + 1] cat_lprob;
        array[nsl[i] + 1] int cat_obs;
        int obs_sum = 0;
        for (j in 1:nsl[i]) {
          int local = obs_lookup[o[1] + j - 1] - c[1] + 1;
          cat_lprob[j] = lprob[local];
          cat_obs[j] = obs[o[1] + j - 1];
          obs_sum += cat_obs[j];
        }
        // Residual probability = 1 - sum(observed probs), on the log scale.
        cat_lprob[nsl[i] + 1] = log1m_exp(log_sum_exp(cat_lprob[1:nsl[i]]));
        cat_obs[nsl[i] + 1] = max(total[i] - obs_sum, 0);
        tar += multinomial_lpmf(cat_obs | exp(cat_lprob));
      }
    }
  }
  return tar;
}