primarycensored_pmf.stan

Go to the documentation of this file.

 
int enw_to_pcens_dist_id(int dist) {
  if (dist == 1) return 2;  // exponential -> primarycensored gamma (shape 1)
  if (dist == 2) return 1;  // lognormal   -> primarycensored lognormal
  if (dist == 3) return 2;  // gamma       -> primarycensored gamma
  reject(
    "The primarycensored discretisation does not support distribution id ",
    dist, ". Supported distributions are exponential, lognormal and gamma."
  );
}
 
array[] real enw_to_pcens_params(real mu, real sigma, int dist) {
  array[2] real params;
  if (dist == 1) {
    // Exponential: epinowcast uses exp(-mu) as the rate. Routed through the
    // gamma path as Gamma(shape = 1, rate = exp(-mu)), which is identical to
    // Exponential(rate = exp(-mu)) but hits the analytical gamma-uniform
    // solution.
    params[1] = 1;
    params[2] = exp(-mu);
  } else if (dist == 2) {
    // Lognormal: shared (meanlog, sdlog) parameterisation.
    params[1] = mu;
    params[2] = sigma;
  } else if (dist == 3) {
    // Gamma: epinowcast uses exp(mu) as the shape with sigma as the rate;
    // primarycensored gamma takes (shape, rate).
    params[1] = exp(mu);
    params[2] = sigma;
  } else {
    reject("Unsupported distribution id for primarycensored params: ", dist);
  }
  return params;
}
 
vector log_hazard_to_logit_hazard(vector lhaz) {
  int n = num_elements(lhaz);
  // Logit transform all but last, then append Inf (last hazard h[n] = 1)
  return append_row(
    lhaz[1:(n-1)] - log1m_exp(lhaz[1:(n-1)]),
    positive_infinity()
  );
}
 
vector lprob_to_log_hazard(vector lprob, int u) {
  vector[u] lhaz;
  // h_0 = p_0 (survival before delay 0 is 1).
  lhaz[1] = lprob[1];
  // Running log survival: log S_d = log(1 - sum_{i <= d} p_i).
  real log_surv = log1m_exp(lprob[1]);
  for (d in 2:u) {
    // h_{d-1} = p_{d-1} / S_{d-2}; here lprob[d] is p_{d-1} and log_surv is
    // log S_{d-2} (survival up to but excluding the current delay).
    lhaz[d] = lprob[d] - log_surv;
    // Update survival to exclude the current delay before the next iteration.
    if (d < u) {
      log_surv = log_diff_exp(log_surv, lprob[d]);
    }
  }
  // Pin the terminal hazard to exactly 1 (log(1) = 0) so that
  // log_hazard_to_logit_hazard() yields +inf and all remaining mass is
  // reported at the maximum delay, rather than relying on floating-point
  // cancellation in the final iteration above.
  lhaz[u] = 0;
  return log_hazard_to_logit_hazard(lhaz);
}
 
vector discretised_pcens_logit_hazard(real mu, real sigma, int dmax, int dist,
                                      int ref_as_p) {
  int pcens_dist = enw_to_pcens_dist_id(dist);
  // All supported distributions (exponential routed via gamma, lognormal,
  // gamma) take two parameters under the primarycensored convention.
  array[2] real params = enw_to_pcens_params(mu, sigma, dist);
  array[0] real primary_params;
  vector[dmax] lprob = primarycensored_sone_lpmf_vectorized(
    dmax - 1, 0.0, dmax * 1.0, pcens_dist, params, 1.0, 1, primary_params
  );
  if (ref_as_p == 1) {
    return lprob;
  }
  return lprob_to_log_hazard(lprob, dmax);
}