adjustments to reflect that max = length - 1

R/estimate_secondary.R

@@ -497,7 +497,7 @@ simulate_secondary <- function(data, type = "incidence", family = "poisson",
       list(i = index, m = meanlog, s = sdlog),
       function(i, m, s) {
         discretised_lognormal_pmf_conv(
-          scaled[max(1, i - delay_max):i],
+          scaled[max(1, i - delay_max - 1):i],
           meanlog = m, sdlog = s
         )
       }

R/utilities.R

@@ -306,17 +306,17 @@ convert_to_logsd <- function(mean, sd) {
 }
 
 discretised_lognormal_pmf <- function(meanlog, sdlog, max_d, reverse = FALSE) {
-  pmf <- plnorm(1:max_d, meanlog, sdlog) -
-    plnorm(0:(max_d - 1), meanlog, sdlog)
-  pmf <- as.vector(pmf) / as.vector(plnorm(max_d, meanlog, sdlog))
+  pmf <- plnorm(1:(max_d + 1), meanlog, sdlog) -
+    plnorm(0:max_d, meanlog, sdlog)
+  pmf <- as.vector(pmf) / as.vector(plnorm(max_d + 1, meanlog, sdlog))
   if (reverse) {
     pmf <- rev(pmf)
   }
   return(pmf)
 }
 
 discretised_lognormal_pmf_conv <- function(x, meanlog, sdlog) {
-  pmf <- discretised_lognormal_pmf(meanlog, sdlog, length(x), reverse = TRUE)
+  pmf <- discretised_lognormal_pmf(meanlog, sdlog, length(x) - 1, reverse = TRUE)
   conv <- sum(x * pmf, na.rm = TRUE)
   return(conv)
 }
@@ -325,9 +325,9 @@ discretised_gamma_pmf <- function(mean, sd, max_d, zero_pad = 0,
                                   reverse = FALSE) {
   alpha <- exp(2 * (log(mean) - log(sd)))
   beta <- exp(log(mean) - 2 * log(sd))
-  pmf <- pgamma(1:max_d, shape = alpha, scale = beta) -
-    pgamma(0:(max_d - 1), shape = alpha, scale = beta)
-  pmf <- as.vector(pmf) / as.vector(pgamma(max_d, shape = alpha, scale = beta))
+  pmf <- pgamma(1:(max_d + 1), shape = alpha, scale = beta) -
+    pgamma(0:max_d, shape = alpha, scale = beta)
+  pmf <- as.vector(pmf) / as.vector(pgamma(max_d + 1, shape = alpha, scale = beta))
   if (zero_pad > 0) {
     pmf <- c(rep(0, zero_pad), pmf)
   }

inst/stan/estimate_infections.stan

@@ -61,15 +61,15 @@ transformed parameters {
   vector[t] infections;                                     // latent infections
   vector[ot_h] reports;                                     // estimated reported cases
   vector[ot] obs_reports;                                   // observed estimated reported cases
-  vector[estimate_r * delay_type_max[gt_id]] gt_rev_pmf;
+  vector[estimate_r * (delay_type_max[gt_id] + 1)] gt_rev_pmf;
   // GP in noise - spectral densities
   if (!fixed) {
     noise = update_gp(PHI, M, L, alpha[1], rho[1], eta, gp_type);
   }
   // Estimate latent infections
   if (estimate_r) {
     gt_rev_pmf = get_delay_rev_pmf(
-      gt_id, delay_type_max[gt_id], delay_types_p, delay_types_id,
+      gt_id, delay_type_max[gt_id] + 1, delay_types_p, delay_types_id,
       delay_types_groups, delay_max, delay_np_pmf,
       delay_np_pmf_groups, delay_mean, delay_sd, delay_dist,
       1, 1, 0
@@ -89,8 +89,8 @@ transformed parameters {
   }
   // convolve from latent infections to mean of observations
   if (delay_id) {
-    vector[delay_type_max[delay_id]] delay_rev_pmf = get_delay_rev_pmf(
-      delay_id, delay_type_max[delay_id], delay_types_p, delay_types_id,
+    vector[delay_type_max[delay_id] + 1] delay_rev_pmf = get_delay_rev_pmf(
+      delay_id, delay_type_max[delay_id] + 1, delay_types_p, delay_types_id,
       delay_types_groups, delay_max, delay_np_pmf,
       delay_np_pmf_groups, delay_mean, delay_sd, delay_dist,
       0, 1, 0
@@ -109,8 +109,8 @@ transformed parameters {
  }
  // truncate near time cases to observed reports
  if (trunc_id) {
-    vector[delay_type_max[trunc_id]] trunc_rev_cmf = get_delay_rev_pmf(
-      trunc_id, delay_type_max[trunc_id], delay_types_p, delay_types_id,
+    vector[delay_type_max[trunc_id] + 1] trunc_rev_cmf = get_delay_rev_pmf(
+      trunc_id, delay_type_max[trunc_id] + 1, delay_types_p, delay_types_id,
       delay_types_groups, delay_max, delay_np_pmf,
       delay_np_pmf_groups, delay_mean, delay_sd, delay_dist,
       0, 1, 1
@@ -171,8 +171,8 @@ generated quantities {
       normal_rng(delay_mean_mean, delay_mean_sd);
     array[delay_n_p] real delay_sd_sample =
       normal_rng(delay_sd_mean, delay_sd_sd);
-    vector[delay_type_max[gt_id]] sampled_gt_rev_pmf = get_delay_rev_pmf(
-      gt_id, delay_type_max[gt_id], delay_types_p, delay_types_id,
+    vector[delay_type_max[gt_id] + 1] sampled_gt_rev_pmf = get_delay_rev_pmf(
+      gt_id, delay_type_max[gt_id] + 1, delay_types_p, delay_types_id,
       delay_types_groups, delay_max, delay_np_pmf,
       delay_np_pmf_groups, delay_mean_sample, delay_sd_sample,
       delay_dist, 1, 1, 0

inst/stan/estimate_secondary.stan

@@ -40,7 +40,7 @@ transformed parameters {
     vector[delay_type_max[delay_id]] delay_rev_pmf;
     if (delay_id) {
       delay_rev_pmf = get_delay_rev_pmf(
-        delay_id, delay_type_max[delay_id], delay_types_p, delay_types_id,
+        delay_id, delay_type_max[delay_id] + 1, delay_types_p, delay_types_id,
         delay_types_groups, delay_max, delay_np_pmf,
         delay_np_pmf_groups, delay_mean, delay_sd, delay_dist,
         0, 1, 0
@@ -61,7 +61,7 @@ transformed parameters {
  // truncate near time cases to observed reports
  if (trunc_id) {
     vector[delay_type_max[trunc_id]] trunc_rev_cmf = get_delay_rev_pmf(
-      trunc_id, delay_type_max[trunc_id], delay_types_p, delay_types_id,
+      trunc_id, delay_type_max[trunc_id] + 1, delay_types_p, delay_types_id,
       delay_types_groups, delay_max, delay_np_pmf,
       delay_np_pmf_groups, delay_mean, delay_sd, delay_dist,
       0, 1, 1

inst/stan/estimate_truncation.stan

@@ -24,7 +24,7 @@ transformed data{
 
   for (i in 1:obs_sets) {
     end_t[i] = t - obs_dist[i];
-    start_t[i] = max(1, end_t[i] - delay_type_max[trunc_id] + 1);
+    start_t[i] = max(1, end_t[i] - delay_type_max[trunc_id]);
   }
 }
 parameters {
@@ -35,11 +35,11 @@ parameters {
 }
 transformed parameters{
   real sqrt_phi = 1 / sqrt(phi);
-  matrix[delay_type_max[trunc_id], obs_sets - 1] trunc_obs = rep_matrix(
-    0, delay_type_max[trunc_id], obs_sets - 1
+  matrix[delay_type_max[trunc_id] + 1, obs_sets - 1] trunc_obs = rep_matrix(
+    0, delay_type_max[trunc_id] + 1, obs_sets - 1
   );
-  vector[delay_type_max[trunc_id]] trunc_rev_cmf = get_delay_rev_pmf(
-    trunc_id, delay_type_max[trunc_id], delay_types_p, delay_types_id,
+  vector[delay_type_max[trunc_id] + 1] trunc_rev_cmf = get_delay_rev_pmf(
+    trunc_id, delay_type_max[trunc_id] + 1, delay_types_p, delay_types_id,
     delay_types_groups, delay_max, delay_np_pmf,
     delay_np_pmf_groups, delay_mean, delay_sd, delay_dist,
     0, 1, 1
@@ -75,10 +75,10 @@ model {
   }
 }
 generated quantities {
-  matrix[delay_type_max[trunc_id], obs_sets] recon_obs = rep_matrix(
-    0, delay_type_max[trunc_id], obs_sets
+  matrix[delay_type_max[trunc_id] + 1, obs_sets] recon_obs = rep_matrix(
+    0, delay_type_max[trunc_id] + 1, obs_sets
   );
-  matrix[delay_type_max[trunc_id], obs_sets - 1] gen_obs;
+  matrix[delay_type_max[trunc_id] + 1, obs_sets - 1] gen_obs;
   // reconstruct all truncated datasets using posterior of the truncation distribution
   for (i in 1:obs_sets) {
     recon_obs[1:(end_t[i] - start_t[i] + 1), i] = truncate(
@@ -87,7 +87,7 @@ generated quantities {
   }
  // generate observations for comparing
   for (i in 1:(obs_sets - 1)) {
-    for (j in 1:delay_type_max[trunc_id]) {
+    for (j in 1:(delay_type_max[trunc_id] + 1)) {
       if (trunc_obs[j, i] == 0) {
         gen_obs[j, i] = 0;
       } else {

inst/stan/functions/delays.stan

@@ -4,10 +4,10 @@ array[] int get_delay_type_max(
 ) {
   array[delay_types] int ret;
   for (i in 1:delay_types) {
-    ret[i] = 1;
+    ret[i] = 0;
     for (j in delay_types_groups[i]:(delay_types_groups[i + 1] - 1)) {
       if (delay_types_p[j]) { // parametric
-        ret[i] += delay_max[delay_types_id[j]] - 1;
+        ret[i] += delay_max[delay_types_id[j]];
       } else { // nonparametric
         ret[i] += delay_np_pmf_groups[delay_types_id[j] + 1] -
           delay_np_pmf_groups[delay_types_id[j]] - 1;
@@ -30,14 +30,14 @@ vector get_delay_rev_pmf(
   int new_len;
   for (i in delay_types_groups[delay_id]:(delay_types_groups[delay_id + 1] - 1)) {
     if (delay_types_p[i]) { // parametric
-      vector[delay_max[delay_types_id[i]]] new_variable_pmf =
+      vector[delay_max[delay_types_id[i]] + 1] new_variable_pmf =
         discretised_pmf(
           delay_mean[delay_types_id[i]],
           delay_sigma[delay_types_id[i]],
-          delay_max[delay_types_id[i]],
+          delay_max[delay_types_id[i]] + 1,
           delay_dist[delay_types_id[i]]
       );
-      new_len = current_len + delay_max[delay_types_id[i]] - 1;
+      new_len = current_len + delay_max[delay_types_id[i]];
       if (current_len == 1) { // first delay
         pmf[1:new_len] = new_variable_pmf;
       } else { // subsequent delay to be convolved

inst/stan/simulate_infections.stan

@@ -40,7 +40,7 @@ generated quantities {
     // generate infections from Rt trace
     vector[delay_type_max[gt_id]] gt_rev_pmf;
     gt_rev_pmf = get_delay_rev_pmf(
-      gt_id, delay_type_max[gt_id], delay_types_p, delay_types_id,
+      gt_id, delay_type_max[gt_id] + 1, delay_types_p, delay_types_id,
       delay_types_groups, delay_max, delay_np_pmf,
       delay_np_pmf_groups, delay_mean[i], delay_sd[i], delay_dist,
       1, 1, 0
@@ -53,7 +53,7 @@ generated quantities {
 
     if (delay_id) {
       vector[delay_type_max[delay_id]] delay_rev_pmf = get_delay_rev_pmf(
-        delay_id, delay_type_max[delay_id], delay_types_p, delay_types_id,
+        delay_id, delay_type_max[delay_id] + 1, delay_types_p, delay_types_id,
         delay_types_groups, delay_max, delay_np_pmf,
         delay_np_pmf_groups, delay_mean[i], delay_sd[i], delay_dist,
         0, 1, 0

inst/stan/simulate_secondary.stan

@@ -33,8 +33,8 @@ generated quantities {
   array[n, all_dates ? t : h] int sim_secondary;
   for (i in 1:n) {
     vector[t] secondary;
-    vector[delay_type_max[delay_id]] delay_rev_pmf = get_delay_rev_pmf(
-        delay_id, delay_type_max[delay_id], delay_types_p, delay_types_id,
+    vector[delay_type_max[delay_id] + 1] delay_rev_pmf = get_delay_rev_pmf(
+        delay_id, delay_type_max[delay_id] + 1, delay_types_p, delay_types_id,
         delay_types_groups, delay_max, delay_np_pmf,
         delay_np_pmf_groups, delay_mean[i], delay_sd[i], delay_dist,
         0, 1, 0

tests/testthat/test-dist.R

@@ -8,18 +8,20 @@ test_that("distributions are the same in R and stan", {
   gamma_params <-
     do.call(gamma_dist_def, (c(args, list(samples = 1))))$params[[1]]
   pmf_r_lognormal <- dist_skel(
-    n = seq_len(args$max_value) - 1, dist = TRUE, cum = FALSE,
+    n = seq_len(args$max_value + 1) - 1, dist = TRUE, cum = FALSE,
     model = "lognormal", params = lognormal_params, max_value = args$max,
     discrete = TRUE
   )
   pmf_r_gamma <- dist_skel(
-    n = seq_len(args$max_value) - 1, dist = TRUE, cum = FALSE,
+    n = seq_len(args$max_value + 1) - 1, dist = TRUE, cum = FALSE,
     model = "gamma", params = gamma_params, max_value = args$max,
     discrete = TRUE
   )
 
-  pmf_stan_lognormal <- discretised_pmf(args$mean, args$sd, args$max_value, 0)
-  pmf_stan_gamma <- discretised_pmf(args$mean, args$sd, args$max_value, 1)
+  pmf_stan_lognormal <- discretised_pmf(
+    args$mean, args$sd, args$max_value + 1, 0
+  )
+  pmf_stan_gamma <- discretised_pmf(args$mean, args$sd, args$max_value + 1, 1)
 
   expect_equal(pmf_r_lognormal, pmf_stan_lognormal)
   expect_equal(pmf_r_gamma, pmf_stan_gamma)