nonparametric.stan

Go to the documentation of this file.

 
real pstep_lcdf(real t, vector boundaries, vector pmf) {
  int K = num_elements(pmf);
  if (t < boundaries[2]) return negative_infinity();
  if (t >= boundaries[K + 1]) return 0;
  // Right-continuous CDF with jumps at the right edges
  // boundaries[2], ..., boundaries[K + 1]. F(t) = cum_pmf[k] for
  // t in [boundaries[k + 1], boundaries[k + 2]); equivalently the
  // largest k with boundaries[k + 1] <= t. Boundary-on-jump cases
  // (t == boundaries[k + 1]) advance k, matching R's
  // `findInterval(left.open = FALSE)`.
  int k = 1;
  while (k < K && boundaries[k + 2] <= t) k += 1;
  return log(cumulative_sum(pmf)[k]);
}
 
vector hazards_to_pmf(vector hazards) {
  int K = num_elements(hazards);
  vector[K] log_surv;
  log_surv[1] = 0;
  if (K > 1) {
    log_surv[2:K] = cumulative_sum(log1m(hazards[1:(K - 1)]));
  }
  return hazards .* exp(log_surv);
}
 
real phazard_lcdf(real t, vector boundaries, vector hazards) {
  return pstep_lcdf(t | boundaries, hazards_to_pmf(hazards));
}
 
vector primary_lcdf_vec(vector p, int primary_id,
                        array[] real primary_params, data real pwindow) {
  int N = num_elements(p);
  vector[N] out;
  for (i in 1:N) {
    out[i] = primary_lcdf(p[i] | primary_id, primary_params, pwindow);
  }
  return out;
}
 
real discretestep_lcdf(
  data real d, vector boundaries, vector pmf,
  int primary_id, array[] real primary_params, data real pwindow
) {
  int K = num_elements(pmf);
  // Integration support in u = d - p for p in [0, pwindow]. It is not
  // clipped at 0 so boundaries that start below zero (delays with negative
  // support) are handled; for non-negative boundaries the per-bin clip to
  // [boundaries[k], boundaries[k + 1]] below gives the same result.
  real u_min = d - pwindow;
  real u_max = d;
 
  // Structural-zero short-circuit. Below `boundaries[2]` F_step is zero
  // and the bin-1 contribution carries `cum_before = 0`, so the integral
  // collapses to 0. Returning `negative_infinity()` directly keeps
  // `log(0)` off the autodiff tape so downstream `log_diff_exp(a, -inf)`
  // evaluates cleanly with a zero gradient w.r.t. `pmf`.
  if (u_max <= boundaries[2]) return negative_infinity();
 
  // Sub-interval endpoints in u-space, clipped to [u_min, u_max].
  vector[K] lo = fmax(u_min, head(boundaries, K));
  vector[K] hi = fmin(u_max, tail(boundaries, K));
 
  // F_step is right-continuous and on [b_k, b_{k+1}) takes the value
  // sum_{j < k} pmf[j] (mass before bin k). cumulative_sum(pmf) gives
  // the mass through and including bin k, so we shift right by one.
  vector[K] cum_before;
  cum_before[1] = 0;
  if (K > 1) cum_before[2:K] = head(cumulative_sum(pmf), K - 1);
 
  // 0/1 mask drops bins with `hi <= lo` from the reduction without a
  // branch in the inner expression. Built on `data`-level inputs.
  vector[K] active;
  for (k in 1:K) active[k] = hi[k] > lo[k] ? 1 : 0;
 
  // F_primary at lo/hi via two vectorised calls; one masked subtraction
  // gives the per-bin difference for the dot product.
  vector[K] f_lo = primary_lcdf_vec(d - lo, primary_id, primary_params,
                                    pwindow);
  vector[K] f_hi = primary_lcdf_vec(d - hi, primary_id, primary_params,
                                    pwindow);
  vector[K] f_diff = (exp(f_lo) - exp(f_hi)) .* active;
 
  real integral = dot_product(cum_before, f_diff);
 
  // Tail region [boundaries[K+1], u_max]: F_step = 1, contributing
  // F_primary(d - tail_start) - F_primary(d - u_max).
  real tail_start = fmax(boundaries[K + 1], u_min);
  if (tail_start < u_max) {
    real fp_tail = exp(primary_lcdf(d - tail_start | primary_id,
                                    primary_params, pwindow));
    real fp_end = exp(primary_lcdf(d - u_max | primary_id,
                                   primary_params, pwindow));
    integral += fp_tail - fp_end;
  }
 
  return log(integral);
}
 
real discretehazard_lcdf(
  data real d, vector boundaries, vector hazards,
  int primary_id, array[] real primary_params, data real pwindow
) {
  return discretestep_lcdf(
    d | boundaries, hazards_to_pmf(hazards), primary_id, primary_params,
    pwindow
  );
}