primary_censored_dist.stan

Go to the documentation of this file.

 
real dist_lcdf(real delay, array[] real params, int dist_id) {
  if (delay <= 0) return negative_infinity();
 
  // Use if-else statements to handle different distribution types
  if (dist_id == 1) return lognormal_lcdf(delay | params[1], params[2]);
  else if (dist_id == 2) return gamma_lcdf(delay | params[1], params[2]);
  else if (dist_id == 3) return normal_lcdf(delay | params[1], params[2]);
  else if (dist_id == 4) return exponential_lcdf(delay | params[1]);
  else if (dist_id == 5) return weibull_lcdf(delay | params[1], params[2]);
  else if (dist_id == 6) return beta_lcdf(delay | params[1], params[2]);
  else if (dist_id == 7) return cauchy_lcdf(delay | params[1], params[2]);
  else if (dist_id == 8) return chi_square_lcdf(delay | params[1]);
  else if (dist_id == 9) return inv_chi_square_lcdf(delay | params[1]);
  else if (dist_id == 10) return double_exponential_lcdf(delay | params[1], params[2]);
  else if (dist_id == 11) return inv_gamma_lcdf(delay | params[1], params[2]);
  else if (dist_id == 12) return logistic_lcdf(delay | params[1], params[2]);
  else if (dist_id == 13) return pareto_lcdf(delay | params[1], params[2]);
  else if (dist_id == 14) return scaled_inv_chi_square_lcdf(delay | params[1], params[2]);
  else if (dist_id == 15) return student_t_lcdf(delay | params[1], params[2], params[3]);
  else if (dist_id == 16) return uniform_lcdf(delay | params[1], params[2]);
  else if (dist_id == 17) return von_mises_lcdf(delay | params[1], params[2]);
  else reject("Invalid distribution identifier");
}
 
real primary_dist_lpdf(real x, int primary_dist_id, array[] real params, real min, real max) {
  // Implement switch for different primary distributions
  if (primary_dist_id == 1) return uniform_lpdf(x | min, max);
  if (primary_dist_id == 2) return expgrowth_lpdf(x | min, max, params[1]);
  // Add more primary distributions as needed
  reject("Invalid primary distribution identifier");
}
 
vector primary_censored_ode(real t, vector y, array[] real theta,
                            array[] real x_r, array[] int x_i) {
  real d = x_r[1];
  int dist_id = x_i[1];
  int primary_dist_id = x_i[2];
  real pwindow = x_r[2];
  int dist_params_len = x_i[3];
  int primary_params_len = x_i[4];
 
  // Extract distribution parameters
  array[dist_params_len] real params;
  if (dist_params_len) {
    params = theta[1:dist_params_len];
  }
  array[primary_params_len] real primary_params;
  if (primary_params_len) {
    int primary_loc = size(theta);
    primary_params = theta[primary_loc - primary_params_len + 1:primary_loc];
  }
 
  real log_cdf = dist_lcdf(t | params, dist_id);
  real log_primary_pdf = primary_dist_lpdf(d - t | primary_dist_id, primary_params, 0, pwindow);
 
  return rep_vector(exp(log_cdf + log_primary_pdf), 1);
}
 
real primary_censored_dist_cdf(data real d, int dist_id, array[] real params,
                               data real pwindow, data real D,
                               int primary_dist_id,
                               array[] real primary_params) {
  real result;
  if (d <= 0) {
    return 0;
  }
 
  if (d >= D) {
    return 1;
  }
 
  // Check if an analytical solution exists
  if (check_for_analytical(dist_id, primary_dist_id)) {
    // Use analytical solution
    result = primary_censored_dist_analytical_cdf(
      d | dist_id, params, pwindow, D, primary_dist_id, primary_params
    );
  } else {
    // Use numerical integration for other cases
    real lower_bound = max({d - pwindow, 1e-6});
    array[size(params) + size(primary_params)] real theta = append_array(params, primary_params);
    array[4] int ids = {dist_id, primary_dist_id, size(params), size(primary_params)};
 
    vector[1] y0 = rep_vector(0.0, 1);
    result = ode_rk45(primary_censored_ode, y0, lower_bound, {d}, theta, {d, pwindow}, ids)[1, 1];
 
    if (!is_inf(D)) {
      real log_cdf_D = primary_censored_dist_lcdf(
        D | dist_id, params, pwindow, positive_infinity(), primary_dist_id,primary_params
      );
      result = exp(log(result) - log_cdf_D);
    }
  }
 
  return result;
}
 
real primary_censored_dist_lcdf(data real d, int dist_id, array[] real params,
                                data real pwindow, data real D,
                                int primary_dist_id,
                                array[] real primary_params) {
  real result;
 
  if (d <= 0) {
    return negative_infinity();
  }
 
  if (d >= D) {
    return 0;
  }
 
  // Check if an analytical solution exists
  if (check_for_analytical(dist_id, primary_dist_id)) {
    result = primary_censored_dist_analytical_lcdf(
      d | dist_id, params, pwindow, positive_infinity(), primary_dist_id, primary_params
    );
  } else {
    // Use numerical integration
    result = log(primary_censored_dist_cdf(
      d | dist_id, params, pwindow, positive_infinity(), primary_dist_id, primary_params
    ));
  }
 
  // Handle truncation
  if (!is_inf(D)) {
    real log_cdf_D = primary_censored_dist_lcdf(
      D | dist_id, params, pwindow, positive_infinity(), primary_dist_id, primary_params
    );
    result = result - log_cdf_D;
  }
 
  return result;
}
 
real primary_censored_dist_lpmf(data int d, int dist_id, array[] real params,
                                data real pwindow, data real d_upper,
                                data real D, int primary_dist_id,
                                array[] real primary_params) {
  if (d_upper > D) {
    reject("Upper truncation point is greater than D. It is ", d_upper,
           " and D is ", D, ". Resolve this by increasing D to be greater or equal to d + swindow or decreasing swindow.");
  }
  if (d_upper <= d) {
    reject("Upper truncation point is less than or equal to d. It is ", d_upper,
           " and d is ", d, ". Resolve this by increasing d to be less than d_upper.");
  }
  real log_cdf_upper = primary_censored_dist_lcdf(
    d_upper | dist_id, params, pwindow, positive_infinity(), primary_dist_id, primary_params
  );
  real log_cdf_lower = primary_censored_dist_lcdf(
    d | dist_id, params, pwindow, positive_infinity(), primary_dist_id, primary_params
  );
  if (!is_inf(D)) {
    real log_cdf_D;
 
    if (d_upper == D) {
      log_cdf_D = log_cdf_upper;
    } else {
      log_cdf_D = primary_censored_dist_lcdf(
        D | dist_id, params, pwindow, positive_infinity(), primary_dist_id, primary_params
      );
    }
    return log_diff_exp(log_cdf_upper, log_cdf_lower) - log_cdf_D;
  } else {
    return log_diff_exp(log_cdf_upper, log_cdf_lower);
  }
}
 
real primary_censored_dist_pmf(data int d, int dist_id, array[] real params,
                               data real pwindow, data real d_upper,
                               data real D, int primary_dist_id,
                               array[] real primary_params) {
  return exp(
    primary_censored_dist_lpmf(
      d | dist_id, params, pwindow, d_upper, D, primary_dist_id, primary_params
    )
  );
}
 
vector primary_censored_sone_lpmf_vectorized(
  int max_delay, data real D, int dist_id,
  array[] real params, data real pwindow,
  int primary_dist_id, array[] real primary_params
) {
 
  int upper_interval = max_delay + 1;
  vector[upper_interval] log_pmfs;
  vector[upper_interval] log_cdfs;
  real log_normalizer;
 
  // Check if D is at least max_delay + 1
  if (D < upper_interval) {
    reject("D must be at least max_delay + 1");
  }
 
  // Compute log CDFs
  for (d in 1:upper_interval) {
    log_cdfs[d] = primary_censored_dist_lcdf(
      d | dist_id, params, pwindow, positive_infinity(), primary_dist_id,
      primary_params
    );
  }
 
  // Compute log normalizer using upper_interval
  if (D > upper_interval) {
    if (is_inf(D)) {
      log_normalizer = 0; // No normalization needed for infinite D
    } else {
      log_normalizer = primary_censored_dist_lcdf(
        D | dist_id, params, pwindow, positive_infinity(),
        primary_dist_id, primary_params
      );
    }
  } else {
    log_normalizer = log_cdfs[upper_interval];
  }
 
  // Compute log PMFs
  log_pmfs[1] = log_cdfs[1] - log_normalizer;
  for (d in 2:upper_interval) {
    log_pmfs[d] = log_diff_exp(log_cdfs[d], log_cdfs[d-1]) - log_normalizer;
  }
 
  return log_pmfs;
}
 
vector primary_censored_sone_pmf_vectorized(
  int max_delay, data real D, int dist_id,
  array[] real params, data real pwindow,
  int primary_dist_id,
  array[] real primary_params
) {
  return exp(
    primary_censored_sone_lpmf_vectorized(
      max_delay, D, dist_id, params, pwindow, primary_dist_id, primary_params
    )
  );
}