From c7d48a1ba05df1179749e87916d15865225d07d2 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Fri, 10 Jan 2025 14:18:30 +0000
Subject: [PATCH 01/34] estimate r from initial R

---
 inst/stan/data/simulation_rt.stan   |  1 -
 inst/stan/estimate_infections.stan  |  8 +----
 inst/stan/functions/infections.stan | 11 ++++---
 inst/stan/functions/rt.stan         | 45 +++++++++++++++++++++++++++++
 inst/stan/simulate_infections.stan  |  2 +-
 5 files changed, 52 insertions(+), 15 deletions(-)

diff --git a/inst/stan/data/simulation_rt.stan b/inst/stan/data/simulation_rt.stan
index 3beae3161..3bad2c61a 100644
--- a/inst/stan/data/simulation_rt.stan
+++ b/inst/stan/data/simulation_rt.stan
@@ -1,5 +1,4 @@
   array[n, 1] real initial_infections; // initial logged infections
-  array[n, seeding_time > 1 ? 1 : 0] real initial_growth; //initial growth
 
   matrix[n, t - seeding_time] R; // reproduction number
   int pop;                       // susceptible population
diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
index fe845b8fe..db1746d75 100644
--- a/inst/stan/estimate_infections.stan
+++ b/inst/stan/estimate_infections.stan
@@ -60,12 +60,6 @@ transformed parameters {
   vector[ot_h] reports;                                     // estimated reported cases
   vector[ot] obs_reports;                                   // observed estimated reported cases
   vector[estimate_r * (delay_type_max[gt_id] + 1)] gt_rev_pmf;
-  array[estimate_r && seeding_time > 1 ? 1 : 0] real initial_growth; // seed growth rate
-
-  if (num_elements(initial_growth) > 0) {
-    initial_growth[1] = initial_growth_estimate +
-      initial_infections_estimate - initial_infections[1];
-  }
 
   // GP in noise - spectral densities
   profile("update gp") {
@@ -108,7 +102,7 @@ transformed parameters {
         params
       );
       infections = generate_infections(
-        R, seeding_time, gt_rev_pmf, initial_infections, initial_growth, pop,
+        R, seeding_time, gt_rev_pmf, initial_infections, pop,
         future_time, obs_scale, frac_obs
       );
     }
diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan
index a599be8c2..5a2dc09bf 100644
--- a/inst/stan/functions/infections.stan
+++ b/inst/stan/functions/infections.stan
@@ -18,14 +18,13 @@ real update_infectiousness(vector infections, vector gt_rev_pmf,
   return(new_inf);
 }
 // generate infections by using Rt = Rt-1 * sum(reversed generation time pmf * infections)
-vector generate_infections(vector oR, int uot, vector gt_rev_pmf,
-                           array[] real initial_infections, array[] real initial_growth,
-                           int pop, int ht, int obs_scale, real frac_obs) {
+vector generate_infections(vector R, int uot, vector gt_rev_pmf,
+                           array[] real initial_infections, int pop, int ht,
+                           int obs_scale, real frac_obs) {
   // time indices and storage
-  int ot = num_elements(oR);
+  int ot = num_elements(R);
   int nht = ot - ht;
   int t = ot + uot;
-  vector[ot] R = oR;
   real exp_adj_Rt;
   vector[t] infections = rep_vector(0, t);
   vector[ot] cum_infections;
@@ -36,7 +35,7 @@ vector generate_infections(vector oR, int uot, vector gt_rev_pmf,
     infections[1] = infections[1] / frac_obs;
   }
   if (uot > 1) {
-    real growth = exp(initial_growth[1]);
+    real growth = exp(R_to_r(R[1], gt_rev_pmf));
     for (s in 2:uot) {
       infections[s] = infections[s - 1] * growth;
     }
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 462e4d259..5bfa42eb3 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -65,3 +65,48 @@ void rt_lp(array[] real initial_infections, vector bp_effects,
   initial_infections ~ normal(prior_infections, sqrt(prior_infections));
   
 }
+
+/**
+ * Calculate the negative moment generating function (nMGF) of a probability
+ * distribution with given probability mass function.
+ *
+ * @param r function argument of the nMGF
+ * @param pmf probability mass function as vector (first index: 0)
+ */
+  real neg_MGF(real r, vector pmf) {
+    int len = num_elements(pmf);
+    vector[len] exp_r = exp(-r * linspaced_vector(len, 0, len - 1));
+    return(dot_product(pmf, exp_r));
+  }
+
+
+/**
+ * Helper function used by the [R_to_r()] function to estimate r from R
+ *
+ * @param r A vector of length 1, the growth rate
+ * @param gt_pmf probability mass function of the generation time
+ * @param R reproduction number
+ */
+vector eq(vector r, vector gt_pmf, real R) {
+  return([ neg_MGF(r[1], gt_pmf) - 1 / R ]');
+}
+
+/**
+ * Estiamte the growth rate r from reproduction number R. Used in the model to
+ * estimate the initial growth rate
+ *
+ * @param R reproduction number
+ * @param gt_rev_pmf reverse probability mass function of the generation time
+ */
+real R_to_r(real R, vector gt_rev_pmf) {
+  int gt_len = num_elements(gt_rev_pmf);
+  vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
+  real mean_gt = dot_product(
+    gt_pmf,
+    linspaced_vector(gt_len, 0, gt_len - 1)
+  );
+  vector[1] r_approx = [ (R - 1) / (R * mean_gt) ]';
+  vector[1] r = solve_newton(eq, r_approx, gt_pmf, R);
+
+  return(r[1]);
+}
diff --git a/inst/stan/simulate_infections.stan b/inst/stan/simulate_infections.stan
index 43ea6cba9..8e3c40aa9 100644
--- a/inst/stan/simulate_infections.stan
+++ b/inst/stan/simulate_infections.stan
@@ -62,7 +62,7 @@ generated quantities {
 
       infections[i] = to_row_vector(generate_infections(
         to_vector(R[i]), seeding_time, gt_rev_pmf, initial_infections[i],
-        initial_growth[i], pop, future_time, obs_scale, frac_obs[i]
+        pop, future_time, obs_scale, frac_obs[i]
       ));
 
       if (delay_id) {

From 07b2b51b59b0ce1e1ac1db4195d77b4bc40e0a96 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Fri, 10 Jan 2025 14:57:39 +0000
Subject: [PATCH 02/34] remove unneeded argument

---
 R/simulate_infections.R | 2 --
 1 file changed, 2 deletions(-)

diff --git a/R/simulate_infections.R b/R/simulate_infections.R
index 6998a8171..8524bf820 100644
--- a/R/simulate_infections.R
+++ b/R/simulate_infections.R
@@ -113,7 +113,6 @@ simulate_infections <- function(estimates, R, initial_infections,
     log_initial_infections <- log(initial_infections) -
       (seeding_time - 1) * initial_growth
   } else {
-    initial_growth <- numeric(0)
     log_initial_infections <- log(initial_infections)
   }
 
@@ -123,7 +122,6 @@ simulate_infections <- function(estimates, R, initial_infections,
     seeding_time = seeding_time,
     future_time = 0,
     initial_infections = array(log_initial_infections, dim = c(1, 1)),
-    initial_growth = array(initial_growth, dim = c(1, length(initial_growth))),
     R = array(R$R, dim = c(1, nrow(R))),
     pop = pop
   )

From 06bd8092345ff308b65379a8538b7edd1694ee5a Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Fri, 10 Jan 2025 14:57:47 +0000
Subject: [PATCH 03/34] update tests

---
 tests/testthat/test-stan-infections.R | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/tests/testthat/test-stan-infections.R b/tests/testthat/test-stan-infections.R
index d08e06c01..fc6633be5 100644
--- a/tests/testthat/test-stan-infections.R
+++ b/tests/testthat/test-stan-infections.R
@@ -25,35 +25,35 @@ gt_rev_pmf <- get_delay_rev_pmf(
 # test generate infections
 test_that("generate_infections works as expected", {
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 0, 0, 0), 0),
     c(rep(1000, 10), 995, 996, rep(997, 8))
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(20), 0.03, 0, 0, 0, 0), 0),
-    c(20, 21, 21, 22, 23, 23, 24, 25, 25, 26, 24, 27, 28, 29, 30, 30, 31, 32, 33, 34)
+    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(20), 0, 0, 0, 0), 0),
+    c(rep(20, 11), 22, 22, 23, 24, 24, 25, 26, 27, 28)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(100), 0, 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(100), 0, 0, 0, 0), 0),
     c(rep(100, 10), 99, 110, 112, 115, 119, 122, 126, 130, 134, 138)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 4, gt_rev_pmf, log(500), -0.02, 0, 0, 0, 0), 0),
-    c(500, 490, 480, 471, 382, 403, 408, rep(409, 7))
+    round(generate_infections(c(1, rep(1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0), 0),
+    c(rep(500, 4), 394, 418, 424, rep(425, 7))
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1.1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0), 0),
     c(rep(500, 4), 394, 460, 475, 489, 505, 520, 536, 553, 570, 588)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 1, gt_rev_pmf, log(40), numeric(0), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1, 9)), 1, gt_rev_pmf, log(40), 0, 0, 0, 0), 0),
     c(40, 8, 11, 12, 12, rep(13, 6))
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 1, gt_rev_pmf, log(100), 0.01, 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1.1, 9)), 1, gt_rev_pmf, log(100), 0, 0, 0, 0), 0),
     c(100, 20, 31, 35, 36, 37, 38, 39, 41, 42, 43)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 100000, 4, 0, 0), 0),
+    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 100000, 4, 0, 0), 0),
     c(rep(1000, 10), 995, 996, rep(997, 4), 980, 965, 947, 926)
   )
 })

From 5239e6308d055717ed5d56ff15ffa928a801a9fe Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Fri, 10 Jan 2025 15:01:16 +0000
Subject: [PATCH 04/34] update naming to be more consistent

---
 R/simulate_infections.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/simulate_infections.R b/R/simulate_infections.R
index 8524bf820..26da448e4 100644
--- a/R/simulate_infections.R
+++ b/R/simulate_infections.R
@@ -108,10 +108,10 @@ simulate_infections <- function(estimates, R, initial_infections,
   if (seeding_time > 1) {
     ## estimate initial growth from initial reproduction number if seeding time
     ## is greater than 1
-    initial_growth <- (R$R[1] - 1) / mean(generation_time)
+    initial_growth_estimate <- (R$R[1] - 1) / mean(generation_time)
     ## adjust initial infections for initial exponential growth
     log_initial_infections <- log(initial_infections) -
-      (seeding_time - 1) * initial_growth
+      (seeding_time - 1) * initial_growth_estimate
   } else {
     log_initial_infections <- log(initial_infections)
   }

From a77ce56923a8b1756bf50859583fe13fdd54b265 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Fri, 10 Jan 2025 15:21:42 +0000
Subject: [PATCH 05/34] don't use sundials solver

---
 inst/stan/functions/rt.stan | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 5bfa42eb3..1f9e4965f 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -106,7 +106,7 @@ real R_to_r(real R, vector gt_rev_pmf) {
     linspaced_vector(gt_len, 0, gt_len - 1)
   );
   vector[1] r_approx = [ (R - 1) / (R * mean_gt) ]';
-  vector[1] r = solve_newton(eq, r_approx, gt_pmf, R);
+  vector[1] r = solve_powell(eq, r_approx, gt_pmf, R);
 
   return(r[1]);
 }

From b4f18b4bf5ab18686b96f5b903d7873c94726c8f Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Fri, 10 Jan 2025 15:21:50 +0000
Subject: [PATCH 06/34] update snapshot

---
 tests/testthat/_snaps/simulate-infections.md | 56 ++++++++++----------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 9250ab91e..6e7ac747d 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -36,33 +36,33 @@
 
               variable       date     value
                 <char>     <Date>     <num>
-     1:     infections 2023-01-01 101.60082
-     2:     infections 2023-01-02 107.29011
-     3:     infections 2023-01-03 113.03527
-     4:     infections 2023-01-04 119.01464
-     5:     infections 2023-01-05 125.28049
-     6:     infections 2023-01-06 131.63823
-     7:     infections 2023-01-07 138.33602
-     8:     infections 2023-01-08  96.91916
-     9:     infections 2023-01-09  94.05215
-    10:     infections 2023-01-10  89.94060
-    11:     infections 2023-01-11  85.64594
-    12:     infections 2023-01-12  81.40206
-    13:     infections 2023-01-13  77.28891
-    14:     infections 2023-01-14  73.33458
-    15: reported_cases 2023-01-01  75.00000
-    16: reported_cases 2023-01-02  67.00000
-    17: reported_cases 2023-01-03  73.00000
-    18: reported_cases 2023-01-04  63.00000
-    19: reported_cases 2023-01-05 276.00000
-    20: reported_cases 2023-01-06 121.00000
-    21: reported_cases 2023-01-07 215.00000
-    22: reported_cases 2023-01-08  84.00000
-    23: reported_cases 2023-01-09 116.00000
-    24: reported_cases 2023-01-10 132.00000
-    25: reported_cases 2023-01-11  91.00000
-    26: reported_cases 2023-01-12  77.00000
-    27: reported_cases 2023-01-13 120.00000
-    28: reported_cases 2023-01-14  41.00000
+     1:     infections 2023-01-01  98.02473
+     2:     infections 2023-01-02 103.42527
+     3:     infections 2023-01-03 108.93812
+     4:     infections 2023-01-04 114.69047
+     5:     infections 2023-01-05 120.72404
+     6:     infections 2023-01-06 126.84017
+     7:     infections 2023-01-07 133.28626
+     8:     infections 2023-01-08  93.37732
+     9:     infections 2023-01-09  90.61201
+    10:     infections 2023-01-10  86.64853
+    11:     infections 2023-01-11  82.50955
+    12:     infections 2023-01-12  78.42043
+    13:     infections 2023-01-13  74.45823
+    14:     infections 2023-01-14  70.65007
+    15: reported_cases 2023-01-01  71.00000
+    16: reported_cases 2023-01-02  64.00000
+    17: reported_cases 2023-01-03  70.00000
+    18: reported_cases 2023-01-04  60.00000
+    19: reported_cases 2023-01-05 266.00000
+    20: reported_cases 2023-01-06 117.00000
+    21: reported_cases 2023-01-07 208.00000
+    22: reported_cases 2023-01-08  81.00000
+    23: reported_cases 2023-01-09  97.00000
+    24: reported_cases 2023-01-10   5.00000
+    25: reported_cases 2023-01-11 109.00000
+    26: reported_cases 2023-01-12  75.00000
+    27: reported_cases 2023-01-13 116.00000
+    28: reported_cases 2023-01-14  40.00000
               variable       date     value
 

From 8c9750f2960b67f98fe67987c46cba65ee3cbc2e Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 09:50:48 +0000
Subject: [PATCH 07/34] estimate initial infections within model

---
 R/create.R                                    | 57 +------------------
 inst/stan/data/rt.stan                        |  2 -
 inst/stan/estimate_infections.stan            |  2 +-
 inst/stan/functions/infections.stan           |  7 ++-
 inst/stan/functions/rt.stan                   |  8 +--
 tests/testthat/test-estimate-early-dynamics.R | 49 ----------------
 6 files changed, 9 insertions(+), 116 deletions(-)
 delete mode 100644 tests/testthat/test-estimate-early-dynamics.R

diff --git a/R/create.R b/R/create.R
index 740000fea..5eac4a5dc 100644
--- a/R/create.R
+++ b/R/create.R
@@ -440,54 +440,6 @@ create_forecast_data <- function(forecast = forecast_opts(), data) {
   return(data)
 }
 
-#' Calculate prior infections and fit early growth
-#'
-#' @description Calculates the prior infections and growth rate based on the
-#' first week's data.
-#'
-#' @param cases Numeric vector; the case counts from the input data.
-#' @inheritParams create_stan_data
-#' @return A list containing `initial_infections_estimate` and
-#'   `initial_growth_estimate`.
-#' @keywords internal
-estimate_early_dynamics <- function(cases, seeding_time) {
-  initial_period <- data.table::data.table(
-    confirm = cases[seq_len(min(7, seeding_time, length(cases)))],
-    t = seq_len(min(7, seeding_time, length(cases))) - 1
-  )[!is.na(confirm)]
-
-  # Calculate initial infections and growth estimate
-  if (seeding_time > 1 && nrow(initial_period) > 1) {
-    safe_lm <- purrr::safely(stats::lm)
-    log_linear_estimate <- safe_lm(log(confirm) ~ t, data = initial_period)[[1]]
-    initial_infections_estimate <- ifelse(
-      is.null(log_linear_estimate), 0, log_linear_estimate$coefficients[1]
-    )
-    initial_growth_estimate <- ifelse(
-      is.null(log_linear_estimate), 0, log_linear_estimate$coefficients[2]
-    )
-  } else {
-    initial_infections_estimate <- 0
-    initial_growth_estimate <- 0
-  }
-
-  # Calculate prior infections
-  if (initial_infections_estimate == 0) {
-    initial_infections_estimate <- log(
-      mean(initial_period$confirm, na.rm = TRUE)
-    )
-    if (is.na(initial_infections_estimate) ||
-        is.null(initial_infections_estimate)) {
-      initial_infections_estimate <- 0
-    }
-  }
-
-  return(list(
-    initial_infections_estimate = initial_infections_estimate,
-    initial_growth_estimate = initial_growth_estimate
-  ))
-}
-
 #' Create Stan Data Required for estimate_infections
 #'
 #' @description`r lifecycle::badge("stable")`
@@ -553,11 +505,6 @@ create_stan_data <- function(data, seeding_time, rt, gp, obs, backcalc,
       delay = stan_data$seeding_time, horizon = stan_data$horizon
     )
   )
-  # calculate prior infections and fit early growth
-  stan_data <- c(
-    stan_data,
-    estimate_early_dynamics(confirmed_cases, seeding_time)
-  )
   # backcalculation settings
   stan_data <- c(stan_data, create_backcalc_data(backcalc))
   # gaussian process data
@@ -639,9 +586,7 @@ create_initial_conditions <- function(data) {
       out$eta <- array(numeric(0))
     }
     if (data$estimate_r == 1) {
-      out$initial_infections <- array(
-        rnorm(1, data$initial_infections_estimate, 0.2)
-      )
+      out$initial_infections <- array(rnorm(1))
     }
 
     if (data$bp_n > 0) {
diff --git a/inst/stan/data/rt.stan b/inst/stan/data/rt.stan
index 357330c19..71890accc 100644
--- a/inst/stan/data/rt.stan
+++ b/inst/stan/data/rt.stan
@@ -1,6 +1,4 @@
   int estimate_r;                    // should the reproduction no be estimated (1 = yes)
-  real initial_infections_estimate;  // point estimate of initial infections
-  real initial_growth_estimate;      // point estimate of initial growth rate
   int bp_n;                          // no of breakpoints (0 = no breakpoints)
   array[t - seeding_time] int breakpoints; // when do breakpoints occur
   int future_fixed;                  // is underlying future Rt assumed to be fixed
diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
index db1746d75..0f6f7be24 100644
--- a/inst/stan/estimate_infections.stan
+++ b/inst/stan/estimate_infections.stan
@@ -204,7 +204,7 @@ model {
     // priors on Rt
     profile("rt lp") {
       rt_lp(
-        initial_infections, bp_effects, bp_sd, bp_n, initial_infections_estimate
+        initial_infections, bp_effects, bp_sd, bp_n
       );
     }
   }
diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan
index 5a2dc09bf..c074c6e32 100644
--- a/inst/stan/functions/infections.stan
+++ b/inst/stan/functions/infections.stan
@@ -29,15 +29,16 @@ vector generate_infections(vector R, int uot, vector gt_rev_pmf,
   vector[t] infections = rep_vector(0, t);
   vector[ot] cum_infections;
   vector[ot] infectiousness;
+  real growth = R_to_r(R[1], gt_rev_pmf);
   // Initialise infections using daily growth
-  infections[1] = exp(initial_infections[1]);
+  infections[1] = exp(initial_infections[1] - growth * uot);
   if (obs_scale) {
     infections[1] = infections[1] / frac_obs;
   }
   if (uot > 1) {
-    real growth = exp(R_to_r(R[1], gt_rev_pmf));
+    real exp_growth = exp(growth);
     for (s in 2:uot) {
-      infections[s] = infections[s - 1] * growth;
+      infections[s] = infections[s - 1] * exp_growth;
     }
   }
   // calculate cumulative infections
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 1f9e4965f..b0366c7e3 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -52,18 +52,16 @@ vector update_Rt(int t, real R0, vector noise, array[] int bps,
  * @param bp_effects Vector of breakpoint effects
  * @param bp_sd Array of breakpoint standard deviations
  * @param bp_n Number of breakpoints
- * @param prior_infections Prior mean for initial infections
  */
 void rt_lp(array[] real initial_infections, vector bp_effects,
-           array[] real bp_sd, int bp_n, real prior_infections) {
+           array[] real bp_sd, int bp_n) {
   //breakpoint effects on Rt
   if (bp_n > 0) {
     bp_sd[1] ~ normal(0, 0.1) T[0,];
     bp_effects ~ normal(0, bp_sd[1]);
   }
-  // initial infections
-  initial_infections ~ normal(prior_infections, sqrt(prior_infections));
-  
+  // initial infections scaling (on the log scale)
+  initial_infections ~ std_normal();
 }
 
 /**
diff --git a/tests/testthat/test-estimate-early-dynamics.R b/tests/testthat/test-estimate-early-dynamics.R
deleted file mode 100644
index 513f1b9df..000000000
--- a/tests/testthat/test-estimate-early-dynamics.R
+++ /dev/null
@@ -1,49 +0,0 @@
-test_that("estimate_early_dynamics works", {
-  cases <- EpiNow2::example_confirmed[1:30]
-  early_estimates <- estimate_early_dynamics(cases$confirm, 7)
-  # Check dimensions
-  expect_identical(
-    names(early_estimates),
-    c("initial_infections_estimate", "initial_growth_estimate")
-  )
-  expect_identical(length(early_estimates), 2L)
-  # Check values
-  expect_identical(
-    round(early_estimates$initial_infections_estimate, 2),
-    3.21
-  )
-  expect_identical(
-    round(early_estimates$initial_growth_estimate, 2),
-    0.35
-  )
-})
-
-test_that("estimate_early_dynamics handles NA values correctly", {
-  cases <- c(10, 20, NA, 40, 50, NA, 70)
-  early_estimates <- estimate_early_dynamics(cases, 7)
-  expect_identical(
-    round(early_estimates$initial_infections_estimate, 2),
-    2.55
-  )
-  expect_true(!is.na(early_estimates$initial_growth_estimate))
-})
-
-test_that("estimate_early_dynamics handles exponential growth", {
-  cases <- 2^(c(0:6)) # Exponential growth
-  early_estimates <- estimate_early_dynamics(cases, 7)
-  expect_equal(early_estimates$initial_infections_estimate, log(2^0))
-  expect_true(early_estimates$initial_growth_estimate > 0) # Growth should be positive
-})
-
-test_that("estimate_early_dynamics handles exponential decline", {
-  cases <- rev(2^(c(0:6))) # Exponential decline
-  early_estimates <- estimate_early_dynamics(cases, 7)
-  expect_equal(early_estimates$initial_infections_estimate, log(2^6))
-  expect_true(early_estimates$initial_growth_estimate < 0) # Growth should be negative
-})
-
-test_that("estimate_early_dynamics correctly handles seeding time less than 2", {
-  cases <- c(5, 10, 20) # Less than 7 days of data
-  early_estimates <- estimate_early_dynamics(cases, 1)
-  expect_equal(early_estimates$initial_growth_estimate, 0) # Growth should be 0 if seeding time is <= 1
-})

From 72340b9fd7f9d61dfc6de1e9ea2d3bc917f94937 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 10:21:37 +0000
Subject: [PATCH 08/34] update simulation snapshot

---
 tests/testthat/_snaps/simulate-infections.md | 112 +++++++++----------
 1 file changed, 56 insertions(+), 56 deletions(-)

diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 6e7ac747d..10583bd23 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -2,67 +2,67 @@
 
               variable       date     value
                 <char>     <Date>     <num>
-     1:     infections 2023-01-01 120.00000
-     2:     infections 2023-01-02 144.00000
-     3:     infections 2023-01-03 172.80000
-     4:     infections 2023-01-04 207.36000
-     5:     infections 2023-01-05 248.83200
-     6:     infections 2023-01-06 298.59840
-     7:     infections 2023-01-07 358.31808
-     8:     infections 2023-01-08 286.65446
-     9:     infections 2023-01-09 229.32357
-    10:     infections 2023-01-10 183.45886
-    11:     infections 2023-01-11 146.76709
-    12:     infections 2023-01-12 117.41367
-    13:     infections 2023-01-13  93.93093
-    14:     infections 2023-01-14  75.14475
-    15: reported_cases 2023-01-01 125.00000
-    16: reported_cases 2023-01-02 135.00000
-    17: reported_cases 2023-01-03 195.00000
-    18: reported_cases 2023-01-04 224.00000
-    19: reported_cases 2023-01-05 253.00000
-    20: reported_cases 2023-01-06 328.00000
-    21: reported_cases 2023-01-07 364.00000
-    22: reported_cases 2023-01-08 278.00000
-    23: reported_cases 2023-01-09 206.00000
-    24: reported_cases 2023-01-10 169.00000
-    25: reported_cases 2023-01-11 144.00000
-    26: reported_cases 2023-01-12 109.00000
-    27: reported_cases 2023-01-13  80.00000
-    28: reported_cases 2023-01-14  95.00000
+     1:     infections 2023-01-01 100.00000
+     2:     infections 2023-01-02 120.00000
+     3:     infections 2023-01-03 144.00000
+     4:     infections 2023-01-04 172.80000
+     5:     infections 2023-01-05 207.36000
+     6:     infections 2023-01-06 248.83200
+     7:     infections 2023-01-07 298.59840
+     8:     infections 2023-01-08 238.87872
+     9:     infections 2023-01-09 191.10298
+    10:     infections 2023-01-10 152.88238
+    11:     infections 2023-01-11 122.30590
+    12:     infections 2023-01-12  97.84472
+    13:     infections 2023-01-13  78.27578
+    14:     infections 2023-01-14  62.62062
+    15: reported_cases 2023-01-01 105.00000
+    16: reported_cases 2023-01-02 112.00000
+    17: reported_cases 2023-01-03 166.00000
+    18: reported_cases 2023-01-04 188.00000
+    19: reported_cases 2023-01-05 211.00000
+    20: reported_cases 2023-01-06 277.00000
+    21: reported_cases 2023-01-07 304.00000
+    22: reported_cases 2023-01-08 231.00000
+    23: reported_cases 2023-01-09 170.00000
+    24: reported_cases 2023-01-10 140.00000
+    25: reported_cases 2023-01-11 120.00000
+    26: reported_cases 2023-01-12  90.00000
+    27: reported_cases 2023-01-13  66.00000
+    28: reported_cases 2023-01-14  77.00000
               variable       date     value
 
 # simulate_infections works as expected with additional parameters
 
               variable       date     value
                 <char>     <Date>     <num>
-     1:     infections 2023-01-01  98.02473
-     2:     infections 2023-01-02 103.42527
-     3:     infections 2023-01-03 108.93812
-     4:     infections 2023-01-04 114.69047
-     5:     infections 2023-01-05 120.72404
-     6:     infections 2023-01-06 126.84017
-     7:     infections 2023-01-07 133.28626
-     8:     infections 2023-01-08  93.37732
-     9:     infections 2023-01-09  90.61201
-    10:     infections 2023-01-10  86.64853
-    11:     infections 2023-01-11  82.50955
-    12:     infections 2023-01-12  78.42043
-    13:     infections 2023-01-13  74.45823
-    14:     infections 2023-01-14  70.65007
-    15: reported_cases 2023-01-01  71.00000
-    16: reported_cases 2023-01-02  64.00000
-    17: reported_cases 2023-01-03  70.00000
-    18: reported_cases 2023-01-04  60.00000
-    19: reported_cases 2023-01-05 266.00000
-    20: reported_cases 2023-01-06 117.00000
-    21: reported_cases 2023-01-07 208.00000
-    22: reported_cases 2023-01-08  81.00000
-    23: reported_cases 2023-01-09  97.00000
-    24: reported_cases 2023-01-10   5.00000
-    25: reported_cases 2023-01-11 109.00000
-    26: reported_cases 2023-01-12  75.00000
-    27: reported_cases 2023-01-13 116.00000
-    28: reported_cases 2023-01-14  40.00000
+     1:     infections 2023-01-01  59.65590
+     2:     infections 2023-01-02  62.94256
+     3:     infections 2023-01-03  66.29757
+     4:     infections 2023-01-04  69.79834
+     5:     infections 2023-01-05  73.47025
+     6:     infections 2023-01-06  77.19240
+     7:     infections 2023-01-07  81.11536
+     8:     infections 2023-01-08  56.82758
+     9:     infections 2023-01-09  55.14467
+    10:     infections 2023-01-10  52.73257
+    11:     infections 2023-01-11  50.21367
+    12:     infections 2023-01-12  47.72511
+    13:     infections 2023-01-13  45.31379
+    14:     infections 2023-01-14  42.99623
+    15: reported_cases 2023-01-01  57.00000
+    16: reported_cases 2023-01-02 128.00000
+    17: reported_cases 2023-01-03  35.00000
+    18: reported_cases 2023-01-04 152.00000
+    19: reported_cases 2023-01-05  67.00000
+    20: reported_cases 2023-01-06 106.00000
+    21: reported_cases 2023-01-07  61.00000
+    22: reported_cases 2023-01-08   3.00000
+    23: reported_cases 2023-01-09  82.00000
+    24: reported_cases 2023-01-10  64.00000
+    25: reported_cases 2023-01-11  84.00000
+    26: reported_cases 2023-01-12  28.00000
+    27: reported_cases 2023-01-13  70.00000
+    28: reported_cases 2023-01-14  42.00000
               variable       date     value
 

From 33692c6e6a6b845c3e359cf756a8b4507dcfedd2 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 11:24:04 +0000
Subject: [PATCH 09/34] another snapshot update

---
 tests/testthat/_snaps/simulate-infections.md | 56 ++++++++++----------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 10583bd23..662ac64e5 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -36,33 +36,33 @@
 
               variable       date     value
                 <char>     <Date>     <num>
-     1:     infections 2023-01-01  59.65590
-     2:     infections 2023-01-02  62.94256
-     3:     infections 2023-01-03  66.29757
-     4:     infections 2023-01-04  69.79834
-     5:     infections 2023-01-05  73.47025
-     6:     infections 2023-01-06  77.19240
-     7:     infections 2023-01-07  81.11536
-     8:     infections 2023-01-08  56.82758
-     9:     infections 2023-01-09  55.14467
-    10:     infections 2023-01-10  52.73257
-    11:     infections 2023-01-11  50.21367
-    12:     infections 2023-01-12  47.72511
-    13:     infections 2023-01-13  45.31379
-    14:     infections 2023-01-14  42.99623
-    15: reported_cases 2023-01-01  57.00000
-    16: reported_cases 2023-01-02 128.00000
-    17: reported_cases 2023-01-03  35.00000
-    18: reported_cases 2023-01-04 152.00000
-    19: reported_cases 2023-01-05  67.00000
-    20: reported_cases 2023-01-06 106.00000
-    21: reported_cases 2023-01-07  61.00000
-    22: reported_cases 2023-01-08   3.00000
-    23: reported_cases 2023-01-09  82.00000
-    24: reported_cases 2023-01-10  64.00000
-    25: reported_cases 2023-01-11  84.00000
-    26: reported_cases 2023-01-12  28.00000
-    27: reported_cases 2023-01-13  70.00000
-    28: reported_cases 2023-01-14  42.00000
+     1:     infections 2023-01-01  97.87994
+     2:     infections 2023-01-02 103.27251
+     3:     infections 2023-01-03 108.77721
+     4:     infections 2023-01-04 114.52107
+     5:     infections 2023-01-05 120.54573
+     6:     infections 2023-01-06 126.65283
+     7:     infections 2023-01-07 133.08940
+     8:     infections 2023-01-08  93.23940
+     9:     infections 2023-01-09  90.47818
+    10:     infections 2023-01-10  86.52055
+    11:     infections 2023-01-11  82.38768
+    12:     infections 2023-01-12  78.30460
+    13:     infections 2023-01-13  74.34825
+    14:     infections 2023-01-14  70.54572
+    15: reported_cases 2023-01-01  71.00000
+    16: reported_cases 2023-01-02  64.00000
+    17: reported_cases 2023-01-03  70.00000
+    18: reported_cases 2023-01-04  60.00000
+    19: reported_cases 2023-01-05 265.00000
+    20: reported_cases 2023-01-06 116.00000
+    21: reported_cases 2023-01-07 208.00000
+    22: reported_cases 2023-01-08  81.00000
+    23: reported_cases 2023-01-09  97.00000
+    24: reported_cases 2023-01-10   5.00000
+    25: reported_cases 2023-01-11 108.00000
+    26: reported_cases 2023-01-12  75.00000
+    27: reported_cases 2023-01-13 116.00000
+    28: reported_cases 2023-01-14  40.00000
               variable       date     value
 

From 693175132328eba886947e3aa40c5c00492cbcbb Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 11:27:21 +0000
Subject: [PATCH 10/34] doc update

---
 man/estimate_early_dynamics.Rd | 23 -----------------------
 1 file changed, 23 deletions(-)
 delete mode 100644 man/estimate_early_dynamics.Rd

diff --git a/man/estimate_early_dynamics.Rd b/man/estimate_early_dynamics.Rd
deleted file mode 100644
index 9f6d19ad9..000000000
--- a/man/estimate_early_dynamics.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/create.R
-\name{estimate_early_dynamics}
-\alias{estimate_early_dynamics}
-\title{Calculate prior infections and fit early growth}
-\usage{
-estimate_early_dynamics(cases, seeding_time)
-}
-\arguments{
-\item{cases}{Numeric vector; the case counts from the input data.}
-
-\item{seeding_time}{Integer; seeding time, usually obtained using
-\code{\link[=get_seeding_time]{get_seeding_time()}}.}
-}
-\value{
-A list containing \code{initial_infections_estimate} and
-\code{initial_growth_estimate}.
-}
-\description{
-Calculates the prior infections and growth rate based on the
-first week's data.
-}
-\keyword{internal}

From 1f8289d99f327abb97a4d22ca858486428266b43 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 13:17:53 +0000
Subject: [PATCH 11/34] update snapshot again

---
 tests/testthat/_snaps/simulate-infections.md | 56 ++++++++++----------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 662ac64e5..10583bd23 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -36,33 +36,33 @@
 
               variable       date     value
                 <char>     <Date>     <num>
-     1:     infections 2023-01-01  97.87994
-     2:     infections 2023-01-02 103.27251
-     3:     infections 2023-01-03 108.77721
-     4:     infections 2023-01-04 114.52107
-     5:     infections 2023-01-05 120.54573
-     6:     infections 2023-01-06 126.65283
-     7:     infections 2023-01-07 133.08940
-     8:     infections 2023-01-08  93.23940
-     9:     infections 2023-01-09  90.47818
-    10:     infections 2023-01-10  86.52055
-    11:     infections 2023-01-11  82.38768
-    12:     infections 2023-01-12  78.30460
-    13:     infections 2023-01-13  74.34825
-    14:     infections 2023-01-14  70.54572
-    15: reported_cases 2023-01-01  71.00000
-    16: reported_cases 2023-01-02  64.00000
-    17: reported_cases 2023-01-03  70.00000
-    18: reported_cases 2023-01-04  60.00000
-    19: reported_cases 2023-01-05 265.00000
-    20: reported_cases 2023-01-06 116.00000
-    21: reported_cases 2023-01-07 208.00000
-    22: reported_cases 2023-01-08  81.00000
-    23: reported_cases 2023-01-09  97.00000
-    24: reported_cases 2023-01-10   5.00000
-    25: reported_cases 2023-01-11 108.00000
-    26: reported_cases 2023-01-12  75.00000
-    27: reported_cases 2023-01-13 116.00000
-    28: reported_cases 2023-01-14  40.00000
+     1:     infections 2023-01-01  59.65590
+     2:     infections 2023-01-02  62.94256
+     3:     infections 2023-01-03  66.29757
+     4:     infections 2023-01-04  69.79834
+     5:     infections 2023-01-05  73.47025
+     6:     infections 2023-01-06  77.19240
+     7:     infections 2023-01-07  81.11536
+     8:     infections 2023-01-08  56.82758
+     9:     infections 2023-01-09  55.14467
+    10:     infections 2023-01-10  52.73257
+    11:     infections 2023-01-11  50.21367
+    12:     infections 2023-01-12  47.72511
+    13:     infections 2023-01-13  45.31379
+    14:     infections 2023-01-14  42.99623
+    15: reported_cases 2023-01-01  57.00000
+    16: reported_cases 2023-01-02 128.00000
+    17: reported_cases 2023-01-03  35.00000
+    18: reported_cases 2023-01-04 152.00000
+    19: reported_cases 2023-01-05  67.00000
+    20: reported_cases 2023-01-06 106.00000
+    21: reported_cases 2023-01-07  61.00000
+    22: reported_cases 2023-01-08   3.00000
+    23: reported_cases 2023-01-09  82.00000
+    24: reported_cases 2023-01-10  64.00000
+    25: reported_cases 2023-01-11  84.00000
+    26: reported_cases 2023-01-12  28.00000
+    27: reported_cases 2023-01-13  70.00000
+    28: reported_cases 2023-01-14  42.00000
               variable       date     value
 

From 998ef5c1026c58a9abd934c6ed291451ff7bad5e Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 15:10:19 +0000
Subject: [PATCH 12/34] try manual approach

---
 inst/stan/functions/infections.stan |  2 +-
 inst/stan/functions/rt.stan         | 43 +++++++++++++++++++++++++++--
 2 files changed, 41 insertions(+), 4 deletions(-)

diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan
index c074c6e32..77eb35fc0 100644
--- a/inst/stan/functions/infections.stan
+++ b/inst/stan/functions/infections.stan
@@ -29,7 +29,7 @@ vector generate_infections(vector R, int uot, vector gt_rev_pmf,
   vector[t] infections = rep_vector(0, t);
   vector[ot] cum_infections;
   vector[ot] infectiousness;
-  real growth = R_to_r(R[1], gt_rev_pmf);
+  real growth = R_to_r_manual(R[1], gt_rev_pmf, 1e-3);
   // Initialise infections using daily growth
   infections[1] = exp(initial_infections[1] - growth * uot);
   if (obs_scale) {
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index b0366c7e3..f5580d4b7 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -77,6 +77,21 @@ void rt_lp(array[] real initial_infections, vector bp_effects,
     return(dot_product(pmf, exp_r));
   }
 
+/**
+ * Helper function for calculating r from R using Newton's method
+ *
+ * @param R Reproduction number
+ * @param r growth rate
+ * @param pmf generation time probability mass function (first index: 0)
+ */
+real R_to_r_newton_step(real R, real r, vector pmf) {
+  int len = num_elements(pmf);
+  vector[len] zero_series = linspaced_vector(len, 0, len - 1);
+  vector[len] exp_r = exp(-r * zero_series);
+  real ret = (R * dot_product(pmf, exp_r) - 1) /
+    (- R * dot_product(pmf .* zero_series, exp_r));
+  return(ret);
+}
 
 /**
  * Helper function used by the [R_to_r()] function to estimate r from R
@@ -90,13 +105,13 @@ vector eq(vector r, vector gt_pmf, real R) {
 }
 
 /**
- * Estiamte the growth rate r from reproduction number R. Used in the model to
- * estimate the initial growth rate
+ * Estimate the growth rate r from reproduction number R. Used in the model to
+ * estimate the initial growth rate using stan's algebraic solver
  *
  * @param R reproduction number
  * @param gt_rev_pmf reverse probability mass function of the generation time
  */
-real R_to_r(real R, vector gt_rev_pmf) {
+real R_to_r(real R, vector gt_rev_pmf, int newton_steps) {
   int gt_len = num_elements(gt_rev_pmf);
   vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
   real mean_gt = dot_product(
@@ -108,3 +123,25 @@ real R_to_r(real R, vector gt_rev_pmf) {
 
   return(r[1]);
 }
+
+/**
+ * Estimate the growth rate r from reproduction number R. Used in the model to
+ * estimate the initial growth rate using Newton's method.
+ *
+ * @param R reproduction number
+ * @param gt_rev_pmf reverse probability mass function of the generation time
+ * @param abs_tol absolute tolerance of the solver
+ */
+real R_to_r_manual(real R, vector gt_rev_pmf, real abs_tol) {
+  int gt_len = num_elements(gt_rev_pmf);
+  vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
+  real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));
+  real r = (R - 1) / (R * mean_gt);
+  real step = abs_tol + 1;
+  while (fabs(step) > abs_tol) {
+    step = R_to_r_newton_step(R, r, gt_pmf);
+    r -= step;
+  }
+
+  return(r);
+}

From 2a006a465f45321a70e7dad41998741030afd186 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 15:41:51 +0000
Subject: [PATCH 13/34] relax prior

---
 inst/stan/functions/rt.stan | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index f5580d4b7..b1ee6b10d 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -61,7 +61,7 @@ void rt_lp(array[] real initial_infections, vector bp_effects,
     bp_effects ~ normal(0, bp_sd[1]);
   }
   // initial infections scaling (on the log scale)
-  initial_infections ~ std_normal();
+  initial_infections ~ normal(0, 5);
 }
 
 /**

From 618f8bce41a086538d09cadc891f3895aa90e7c9 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 16:06:15 +0000
Subject: [PATCH 14/34] scale with early cases

---
 inst/stan/estimate_infections.stan | 3 ++-
 inst/stan/functions/rt.stan        | 4 ++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
index 0f6f7be24..4cd3dc703 100644
--- a/inst/stan/estimate_infections.stan
+++ b/inst/stan/estimate_infections.stan
@@ -204,7 +204,8 @@ model {
     // priors on Rt
     profile("rt lp") {
       rt_lp(
-        initial_infections, bp_effects, bp_sd, bp_n
+        initial_infections, bp_effects, bp_sd, bp_n,
+        fmax(1, mean(head(cases, 7)))
       );
     }
   }
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index b1ee6b10d..064b3eb6f 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -54,14 +54,14 @@ vector update_Rt(int t, real R0, vector noise, array[] int bps,
  * @param bp_n Number of breakpoints
  */
 void rt_lp(array[] real initial_infections, vector bp_effects,
-           array[] real bp_sd, int bp_n) {
+           array[] real bp_sd, int bp_n, real prior_infections) {
   //breakpoint effects on Rt
   if (bp_n > 0) {
     bp_sd[1] ~ normal(0, 0.1) T[0,];
     bp_effects ~ normal(0, bp_sd[1]);
   }
   // initial infections scaling (on the log scale)
-  initial_infections ~ normal(0, 5);
+  initial_infections ~ normal(log(prior_infections), 2);
 }
 
 /**

From 487f26037bfac88610a0f016e079ab37524a8560 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 16:59:06 +0000
Subject: [PATCH 15/34] update snapshots

---
 tests/testthat/_snaps/simulate-infections.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 10583bd23..9fd28840d 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -12,7 +12,7 @@
      8:     infections 2023-01-08 238.87872
      9:     infections 2023-01-09 191.10298
     10:     infections 2023-01-10 152.88238
-    11:     infections 2023-01-11 122.30590
+    11:     infections 2023-01-11 122.30591
     12:     infections 2023-01-12  97.84472
     13:     infections 2023-01-13  78.27578
     14:     infections 2023-01-14  62.62062
@@ -37,18 +37,18 @@
               variable       date     value
                 <char>     <Date>     <num>
      1:     infections 2023-01-01  59.65590
-     2:     infections 2023-01-02  62.94256
+     2:     infections 2023-01-02  62.94257
      3:     infections 2023-01-03  66.29757
      4:     infections 2023-01-04  69.79834
      5:     infections 2023-01-05  73.47025
-     6:     infections 2023-01-06  77.19240
-     7:     infections 2023-01-07  81.11536
+     6:     infections 2023-01-06  77.19241
+     7:     infections 2023-01-07  81.11537
      8:     infections 2023-01-08  56.82758
      9:     infections 2023-01-09  55.14467
     10:     infections 2023-01-10  52.73257
     11:     infections 2023-01-11  50.21367
     12:     infections 2023-01-12  47.72511
-    13:     infections 2023-01-13  45.31379
+    13:     infections 2023-01-13  45.31380
     14:     infections 2023-01-14  42.99623
     15: reported_cases 2023-01-01  57.00000
     16: reported_cases 2023-01-02 128.00000

From 8986fb102042ffa40b2a056ab339e2675bbf3a45 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 16:59:15 +0000
Subject: [PATCH 16/34] fabs -> abs

---
 inst/stan/functions/rt.stan | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 064b3eb6f..fca2c3f3c 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -138,7 +138,7 @@ real R_to_r_manual(real R, vector gt_rev_pmf, real abs_tol) {
   real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));
   real r = (R - 1) / (R * mean_gt);
   real step = abs_tol + 1;
-  while (fabs(step) > abs_tol) {
+  while (abs(step) > abs_tol) {
     step = R_to_r_newton_step(R, r, gt_pmf);
     r -= step;
   }

From 60d7cd5ce42a59a5d4a11926f2d29d46043dd5d3 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 13 Jan 2025 17:33:03 +0000
Subject: [PATCH 17/34] pass cases

---
 inst/stan/estimate_infections.stan | 2 +-
 inst/stan/functions/rt.stan        | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
index 4cd3dc703..5d9010e4c 100644
--- a/inst/stan/estimate_infections.stan
+++ b/inst/stan/estimate_infections.stan
@@ -205,7 +205,7 @@ model {
     profile("rt lp") {
       rt_lp(
         initial_infections, bp_effects, bp_sd, bp_n,
-        fmax(1, mean(head(cases, 7)))
+        cases
       );
     }
   }
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index fca2c3f3c..22934adde 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -54,14 +54,16 @@ vector update_Rt(int t, real R0, vector noise, array[] int bps,
  * @param bp_n Number of breakpoints
  */
 void rt_lp(array[] real initial_infections, vector bp_effects,
-           array[] real bp_sd, int bp_n, real prior_infections) {
+           array[] real bp_sd, int bp_n, array[] int cases) {
   //breakpoint effects on Rt
   if (bp_n > 0) {
     bp_sd[1] ~ normal(0, 0.1) T[0,];
     bp_effects ~ normal(0, bp_sd[1]);
   }
   // initial infections scaling (on the log scale)
-  initial_infections ~ normal(log(prior_infections), 2);
+  int head_count = num_elements(cases) > 7 ? 7 : num_elements(cases);
+  real infections_guess = mean(head(cases, head_count));
+  initial_infections ~ normal(log(infections_guess), 2);
 }
 
 /**

From 9ed36e3ab62a653aa596903b55eabd26e94328ce Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Tue, 14 Jan 2025 10:56:01 +0000
Subject: [PATCH 18/34] ensure mean init is >=1

---
 inst/stan/functions/rt.stan | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 22934adde..aff8bc505 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -62,7 +62,7 @@ void rt_lp(array[] real initial_infections, vector bp_effects,
   }
   // initial infections scaling (on the log scale)
   int head_count = num_elements(cases) > 7 ? 7 : num_elements(cases);
-  real infections_guess = mean(head(cases, head_count));
+  real infections_guess = fmax(1, mean(head(cases, head_count)));
   initial_infections ~ normal(log(infections_guess), 2);
 }
 

From 9beca9ab0300a6a18297e844807b8e80a88a322b Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 20 Jan 2025 07:48:45 +0000
Subject: [PATCH 19/34] move to transformed data

---
 inst/stan/estimate_infections.stan | 8 +++++++-
 inst/stan/functions/rt.stan        | 6 ++----
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
index 5d9010e4c..9a8e4112b 100644
--- a/inst/stan/estimate_infections.stan
+++ b/inst/stan/estimate_infections.stan
@@ -38,6 +38,12 @@ transformed data {
       delay_types_groups, delay_max, delay_np_pmf_groups
     );
   }
+
+  // initial infections scaling (on the log scale)
+  real initial_infections_guess = fmax(
+    1,
+    mean(head(cases, num_elements(cases) > 7 ? 7 : num_elements(cases)))
+  );
 }
 
 parameters {
@@ -205,7 +211,7 @@ model {
     profile("rt lp") {
       rt_lp(
         initial_infections, bp_effects, bp_sd, bp_n,
-        cases
+        cases, initial_infections_guess
       );
     }
   }
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index aff8bc505..798fc6e24 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -54,15 +54,13 @@ vector update_Rt(int t, real R0, vector noise, array[] int bps,
  * @param bp_n Number of breakpoints
  */
 void rt_lp(array[] real initial_infections, vector bp_effects,
-           array[] real bp_sd, int bp_n, array[] int cases) {
+           array[] real bp_sd, int bp_n, array[] int cases,
+           real infections_guess) {
   //breakpoint effects on Rt
   if (bp_n > 0) {
     bp_sd[1] ~ normal(0, 0.1) T[0,];
     bp_effects ~ normal(0, bp_sd[1]);
   }
-  // initial infections scaling (on the log scale)
-  int head_count = num_elements(cases) > 7 ? 7 : num_elements(cases);
-  real infections_guess = fmax(1, mean(head(cases, head_count)));
   initial_infections ~ normal(log(infections_guess), 2);
 }
 

From 20e4f6f2b1a0609bba99522018a7ccbf58fde4b6 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 20 Jan 2025 07:49:05 +0000
Subject: [PATCH 20/34] add source

---
 inst/stan/functions/rt.stan | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 798fc6e24..9bfa2d6bc 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -68,6 +68,10 @@ void rt_lp(array[] real initial_infections, vector bp_effects,
  * Calculate the negative moment generating function (nMGF) of a probability
  * distribution with given probability mass function.
  *
+ * Function code is based on Julia code from
+ * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
+ * under Apache license 2.0.
+ *
  * @param r function argument of the nMGF
  * @param pmf probability mass function as vector (first index: 0)
  */
@@ -80,6 +84,10 @@ void rt_lp(array[] real initial_infections, vector bp_effects,
 /**
  * Helper function for calculating r from R using Newton's method
  *
+ * Function code is based on Julia code from
+ * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
+ * under Apache license 2.0.
+ *
  * @param R Reproduction number
  * @param r growth rate
  * @param pmf generation time probability mass function (first index: 0)
@@ -96,6 +104,10 @@ real R_to_r_newton_step(real R, real r, vector pmf) {
 /**
  * Helper function used by the [R_to_r()] function to estimate r from R
  *
+ * Function code is based on Julia code from
+ * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
+ * under Apache license 2.0.
+ *
  * @param r A vector of length 1, the growth rate
  * @param gt_pmf probability mass function of the generation time
  * @param R reproduction number
@@ -108,6 +120,10 @@ vector eq(vector r, vector gt_pmf, real R) {
  * Estimate the growth rate r from reproduction number R. Used in the model to
  * estimate the initial growth rate using stan's algebraic solver
  *
+ * Function code is based on Julia code from
+ * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
+ * under Apache license 2.0.
+ *
  * @param R reproduction number
  * @param gt_rev_pmf reverse probability mass function of the generation time
  */
@@ -128,6 +144,10 @@ real R_to_r(real R, vector gt_rev_pmf, int newton_steps) {
  * Estimate the growth rate r from reproduction number R. Used in the model to
  * estimate the initial growth rate using Newton's method.
  *
+ * Function code is based on Julia code from
+ * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
+ * under Apache license 2.0.
+ *
  * @param R reproduction number
  * @param gt_rev_pmf reverse probability mass function of the generation time
  * @param abs_tol absolute tolerance of the solver

From 27fa43de3e204f6652ffe8585c23daf98847a31b Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 20 Jan 2025 07:49:13 +0000
Subject: [PATCH 21/34] update news item

---
 NEWS.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 02fbc8cd6..4d5e1206d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -20,7 +20,7 @@
 - All parameters have been changed to the new parameter interface. By @sbfnk in #871 and #890 and reviewed by @seabbs.
 - The Gaussian Process lengthscale is now scaled internally by half the length of the time series. By @sbfnk in #890 and reviewed by #seabbs.
 - A bug was fixed where `plot.dist_spec()` wasn't throwing an informative error due to an incomplete check for the max of the specified delay. By @jamesmbaazam in #858 and reviewed by @.
-- Updated the early dynamics calculation to use the full linear model if available. Also changed the prior for initial infections to be approximately Poisson and the initial growth rate to the point estimate of the initial growth rate scaled linearly by the estimated initial infections term. By @sbfnk in #903 and reviewed by @seabbs and @SamuelBrand1 
+- Updated the early dynamics calculation to estimate growth from the initial reproduction number instead of a separate linear model. Also changed the prior calculation for initial infections to be a scaling factor of early case numbers adjusted by the growth estimate, instead a true number of initial infections. By @sbfnk in #923 (with initial exploration in #903) and reviewed by @seabbs and @SamuelBrand1. 
 
 ## Package changes
 

From d3b63f390ad73325e3cc34faef3ba9555420a61e Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 20 Jan 2025 08:01:36 +0000
Subject: [PATCH 22/34] fix simulations

---
 R/simulate_infections.R             | 17 +++++------------
 inst/stan/data/simulation_rt.stan   |  1 +
 inst/stan/estimate_infections.stan  |  6 +++---
 inst/stan/functions/infections.stan | 13 +++++++++----
 inst/stan/functions/rt.stan         |  8 ++++----
 inst/stan/simulate_infections.stan  |  2 +-
 6 files changed, 23 insertions(+), 24 deletions(-)

diff --git a/R/simulate_infections.R b/R/simulate_infections.R
index 26da448e4..71c00dfcc 100644
--- a/R/simulate_infections.R
+++ b/R/simulate_infections.R
@@ -105,23 +105,14 @@ simulate_infections <- function(estimates, R, initial_infections,
   if (missing(seeding_time)) {
     seeding_time <- sum(max(generation_time))
   }
-  if (seeding_time > 1) {
-    ## estimate initial growth from initial reproduction number if seeding time
-    ## is greater than 1
-    initial_growth_estimate <- (R$R[1] - 1) / mean(generation_time)
-    ## adjust initial infections for initial exponential growth
-    log_initial_infections <- log(initial_infections) -
-      (seeding_time - 1) * initial_growth_estimate
-  } else {
-    log_initial_infections <- log(initial_infections)
-  }
 
   data <- list(
     n = 1,
     t = nrow(R) + seeding_time,
     seeding_time = seeding_time,
     future_time = 0,
-    initial_infections = array(log_initial_infections, dim = c(1, 1)),
+    initial_infections = array(log(initial_infections), dim = c(1, 1)),
+    initial_as_scale = 0,
     R = array(R$R, dim = c(1, nrow(R))),
     pop = pop
   )
@@ -431,7 +422,9 @@ forecast_infections <- function(estimates,
     draws <- map(draws, ~ as.matrix(.[nstart:nend, ]))
 
     ## prepare data for stan command
-    data <- c(list(n = dim(draws$R)[1]), draws, estimates$args)
+    data <- c(
+      list(n = dim(draws$R)[1], initial_as_scale = 1), draws, estimates$args
+    )
 
     ## allocate empty parameters
     data <- allocate_empty(
diff --git a/inst/stan/data/simulation_rt.stan b/inst/stan/data/simulation_rt.stan
index 3bad2c61a..ff2b1c933 100644
--- a/inst/stan/data/simulation_rt.stan
+++ b/inst/stan/data/simulation_rt.stan
@@ -1,4 +1,5 @@
   array[n, 1] real initial_infections; // initial logged infections
+  int initial_as_scale; // whether to interpret initial infections as scaling
 
   matrix[n, t - seeding_time] R; // reproduction number
   int pop;                       // susceptible population
diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
index 9a8e4112b..bceabfd8e 100644
--- a/inst/stan/estimate_infections.stan
+++ b/inst/stan/estimate_infections.stan
@@ -41,8 +41,8 @@ transformed data {
 
   // initial infections scaling (on the log scale)
   real initial_infections_guess = fmax(
-    1,
-    mean(head(cases, num_elements(cases) > 7 ? 7 : num_elements(cases)))
+    0,
+    log(mean(head(cases, num_elements(cases) > 7 ? 7 : num_elements(cases))))
   );
 }
 
@@ -109,7 +109,7 @@ transformed parameters {
       );
       infections = generate_infections(
         R, seeding_time, gt_rev_pmf, initial_infections, pop,
-        future_time, obs_scale, frac_obs
+        future_time, obs_scale, frac_obs, 1
       );
     }
   } else {
diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan
index 77eb35fc0..f9de031eb 100644
--- a/inst/stan/functions/infections.stan
+++ b/inst/stan/functions/infections.stan
@@ -17,10 +17,11 @@ real update_infectiousness(vector infections, vector gt_rev_pmf,
   );
   return(new_inf);
 }
+
 // generate infections by using Rt = Rt-1 * sum(reversed generation time pmf * infections)
 vector generate_infections(vector R, int uot, vector gt_rev_pmf,
                            array[] real initial_infections, int pop, int ht,
-                           int obs_scale, real frac_obs) {
+                           int obs_scale, real frac_obs, int initial_as_scale) {
   // time indices and storage
   int ot = num_elements(R);
   int nht = ot - ht;
@@ -31,9 +32,13 @@ vector generate_infections(vector R, int uot, vector gt_rev_pmf,
   vector[ot] infectiousness;
   real growth = R_to_r_manual(R[1], gt_rev_pmf, 1e-3);
   // Initialise infections using daily growth
-  infections[1] = exp(initial_infections[1] - growth * uot);
-  if (obs_scale) {
-    infections[1] = infections[1] / frac_obs;
+  if (initial_as_scale) {
+    infections[1] = exp(initial_infections[1] - growth * uot);
+    if (obs_scale) {
+      infections[1] = infections[1] / frac_obs;
+    }
+  } else {
+    infections[1] = exp(initial_infections[1]);
   }
   if (uot > 1) {
     real exp_growth = exp(growth);
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 9bfa2d6bc..516992b65 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -48,20 +48,20 @@ vector update_Rt(int t, real R0, vector noise, array[] int bps,
 /**
  * Calculate the log-probability of the reproduction number (Rt) priors
  *
- * @param initial_infections Array of initial infection values
+ * @param initial_infections_scale Array of initial infection values
  * @param bp_effects Vector of breakpoint effects
  * @param bp_sd Array of breakpoint standard deviations
  * @param bp_n Number of breakpoints
  */
-void rt_lp(array[] real initial_infections, vector bp_effects,
+void rt_lp(array[] real initial_infections_scale, vector bp_effects,
            array[] real bp_sd, int bp_n, array[] int cases,
-           real infections_guess) {
+           real initial_infections_guess) {
   //breakpoint effects on Rt
   if (bp_n > 0) {
     bp_sd[1] ~ normal(0, 0.1) T[0,];
     bp_effects ~ normal(0, bp_sd[1]);
   }
-  initial_infections ~ normal(log(infections_guess), 2);
+  initial_infections_scale ~ normal(initial_infections_guess, 2);
 }
 
 /**
diff --git a/inst/stan/simulate_infections.stan b/inst/stan/simulate_infections.stan
index 8e3c40aa9..f93c46357 100644
--- a/inst/stan/simulate_infections.stan
+++ b/inst/stan/simulate_infections.stan
@@ -62,7 +62,7 @@ generated quantities {
 
       infections[i] = to_row_vector(generate_infections(
         to_vector(R[i]), seeding_time, gt_rev_pmf, initial_infections[i],
-        pop, future_time, obs_scale, frac_obs[i]
+        pop, future_time, obs_scale, frac_obs[i], initial_as_scale
       ));
 
       if (delay_id) {

From 31ebdc9b84b986c1d9dfd9a9961521c461fc85b8 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 20 Jan 2025 09:08:29 +0000
Subject: [PATCH 23/34] update tests

---
 tests/testthat/test-stan-infections.R | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/tests/testthat/test-stan-infections.R b/tests/testthat/test-stan-infections.R
index fc6633be5..2d14885ed 100644
--- a/tests/testthat/test-stan-infections.R
+++ b/tests/testthat/test-stan-infections.R
@@ -25,35 +25,35 @@ gt_rev_pmf <- get_delay_rev_pmf(
 # test generate infections
 test_that("generate_infections works as expected", {
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 0, 0, 0, 1), 0),
     c(rep(1000, 10), 995, 996, rep(997, 8))
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(20), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(20), 0, 0, 0, 0, 1), 0),
     c(rep(20, 11), 22, 22, 23, 24, 24, 25, 26, 27, 28)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(100), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(100), 0, 0, 0, 0, 1), 0),
     c(rep(100, 10), 99, 110, 112, 115, 119, 122, 126, 130, 134, 138)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0, 1), 0),
     c(rep(500, 4), 394, 418, 424, rep(425, 7))
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1.1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0, 1), 0),
     c(rep(500, 4), 394, 460, 475, 489, 505, 520, 536, 553, 570, 588)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 1, gt_rev_pmf, log(40), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1, 9)), 1, gt_rev_pmf, log(40), 0, 0, 0, 0, 1), 0),
     c(40, 8, 11, 12, 12, rep(13, 6))
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1.1, 9)), 1, gt_rev_pmf, log(100), 0, 0, 0, 0), 0),
+    round(generate_infections(c(1, rep(1.1, 9)), 1, gt_rev_pmf, log(100), 0, 0, 0, 0, 1), 0),
     c(100, 20, 31, 35, 36, 37, 38, 39, 41, 42, 43)
   )
   expect_equal(
-    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 100000, 4, 0, 0), 0),
+    round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 100000, 4, 0, 0, 1), 0),
     c(rep(1000, 10), 995, 996, rep(997, 4), 980, 965, 947, 926)
   )
 })

From 11534c7a199c5b11c774477012c97ad79de0560e Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 20 Jan 2025 09:08:35 +0000
Subject: [PATCH 24/34] update snapshots

---
 tests/testthat/_snaps/simulate-infections.md | 118 +++++++++----------
 1 file changed, 59 insertions(+), 59 deletions(-)

diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 9fd28840d..8d44802eb 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -2,67 +2,67 @@
 
               variable       date     value
                 <char>     <Date>     <num>
-     1:     infections 2023-01-01 100.00000
-     2:     infections 2023-01-02 120.00000
-     3:     infections 2023-01-03 144.00000
-     4:     infections 2023-01-04 172.80000
-     5:     infections 2023-01-05 207.36000
-     6:     infections 2023-01-06 248.83200
-     7:     infections 2023-01-07 298.59840
-     8:     infections 2023-01-08 238.87872
-     9:     infections 2023-01-09 191.10298
-    10:     infections 2023-01-10 152.88238
-    11:     infections 2023-01-11 122.30591
-    12:     infections 2023-01-12  97.84472
-    13:     infections 2023-01-13  78.27578
-    14:     infections 2023-01-14  62.62062
-    15: reported_cases 2023-01-01 105.00000
-    16: reported_cases 2023-01-02 112.00000
-    17: reported_cases 2023-01-03 166.00000
-    18: reported_cases 2023-01-04 188.00000
-    19: reported_cases 2023-01-05 211.00000
-    20: reported_cases 2023-01-06 277.00000
-    21: reported_cases 2023-01-07 304.00000
-    22: reported_cases 2023-01-08 231.00000
-    23: reported_cases 2023-01-09 170.00000
-    24: reported_cases 2023-01-10 140.00000
-    25: reported_cases 2023-01-11 120.00000
-    26: reported_cases 2023-01-12  90.00000
-    27: reported_cases 2023-01-13  66.00000
-    28: reported_cases 2023-01-14  77.00000
+     1:     infections 2023-01-01 120.00000
+     2:     infections 2023-01-02 144.00000
+     3:     infections 2023-01-03 172.80000
+     4:     infections 2023-01-04 207.36000
+     5:     infections 2023-01-05 248.83200
+     6:     infections 2023-01-06 298.59840
+     7:     infections 2023-01-07 358.31808
+     8:     infections 2023-01-08 286.65446
+     9:     infections 2023-01-09 229.32357
+    10:     infections 2023-01-10 183.45886
+    11:     infections 2023-01-11 146.76709
+    12:     infections 2023-01-12 117.41367
+    13:     infections 2023-01-13  93.93093
+    14:     infections 2023-01-14  75.14475
+    15: reported_cases 2023-01-01 125.00000
+    16: reported_cases 2023-01-02 135.00000
+    17: reported_cases 2023-01-03 195.00000
+    18: reported_cases 2023-01-04 224.00000
+    19: reported_cases 2023-01-05 253.00000
+    20: reported_cases 2023-01-06 328.00000
+    21: reported_cases 2023-01-07 364.00000
+    22: reported_cases 2023-01-08 278.00000
+    23: reported_cases 2023-01-09 206.00000
+    24: reported_cases 2023-01-10 169.00000
+    25: reported_cases 2023-01-11 144.00000
+    26: reported_cases 2023-01-12 109.00000
+    27: reported_cases 2023-01-13  80.00000
+    28: reported_cases 2023-01-14  95.00000
               variable       date     value
 
 # simulate_infections works as expected with additional parameters
 
-              variable       date     value
-                <char>     <Date>     <num>
-     1:     infections 2023-01-01  59.65590
-     2:     infections 2023-01-02  62.94257
-     3:     infections 2023-01-03  66.29757
-     4:     infections 2023-01-04  69.79834
-     5:     infections 2023-01-05  73.47025
-     6:     infections 2023-01-06  77.19241
-     7:     infections 2023-01-07  81.11537
-     8:     infections 2023-01-08  56.82758
-     9:     infections 2023-01-09  55.14467
-    10:     infections 2023-01-10  52.73257
-    11:     infections 2023-01-11  50.21367
-    12:     infections 2023-01-12  47.72511
-    13:     infections 2023-01-13  45.31380
-    14:     infections 2023-01-14  42.99623
-    15: reported_cases 2023-01-01  57.00000
-    16: reported_cases 2023-01-02 128.00000
-    17: reported_cases 2023-01-03  35.00000
-    18: reported_cases 2023-01-04 152.00000
-    19: reported_cases 2023-01-05  67.00000
-    20: reported_cases 2023-01-06 106.00000
-    21: reported_cases 2023-01-07  61.00000
-    22: reported_cases 2023-01-08   3.00000
-    23: reported_cases 2023-01-09  82.00000
-    24: reported_cases 2023-01-10  64.00000
-    25: reported_cases 2023-01-11  84.00000
-    26: reported_cases 2023-01-12  28.00000
-    27: reported_cases 2023-01-13  70.00000
-    28: reported_cases 2023-01-14  42.00000
-              variable       date     value
+              variable       date    value
+                <char>     <Date>    <num>
+     1:     infections 2023-01-01 160.8333
+     2:     infections 2023-01-02 169.6942
+     3:     infections 2023-01-03 178.7393
+     4:     infections 2023-01-04 188.1775
+     5:     infections 2023-01-05 198.0770
+     6:     infections 2023-01-06 208.1120
+     7:     infections 2023-01-07 218.6883
+     8:     infections 2023-01-08 153.2081
+     9:     infections 2023-01-09 148.6709
+    10:     infections 2023-01-10 142.1679
+    11:     infections 2023-01-11 135.3769
+    12:     infections 2023-01-12 128.6677
+    13:     infections 2023-01-13 122.1667
+    14:     infections 2023-01-14 115.9185
+    15: reported_cases 2023-01-01 117.0000
+    16: reported_cases 2023-01-02 107.0000
+    17: reported_cases 2023-01-03 119.0000
+    18: reported_cases 2023-01-04  99.0000
+    19: reported_cases 2023-01-05 439.0000
+    20: reported_cases 2023-01-06 192.0000
+    21: reported_cases 2023-01-07 328.0000
+    22: reported_cases 2023-01-08 130.0000
+    23: reported_cases 2023-01-09 161.0000
+    24: reported_cases 2023-01-10  12.0000
+    25: reported_cases 2023-01-11 179.0000
+    26: reported_cases 2023-01-12  84.0000
+    27: reported_cases 2023-01-13 183.0000
+    28: reported_cases 2023-01-14  64.0000
+              variable       date    value
 

From c62be58028a138f7b58b734fa58675f745707f11 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Tue, 21 Jan 2025 08:29:12 +0100
Subject: [PATCH 25/34] temporarily remove additional repos

---
 DESCRIPTION | 2 --
 1 file changed, 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 396d0c3a1..739dedf18 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -146,8 +146,6 @@ LinkingTo:
     RcppParallel (>= 5.0.1),
     rstan (>= 2.26.0),
     StanHeaders (>= 2.26.0)
-Additional_repositories:
-    https://stan-dev.r-universe.dev
 Biarch: true
 Config/testthat/edition: 3
 Encoding: UTF-8

From 9dd97d508af41e52b785a5ff0052c902df11f6db Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Tue, 21 Jan 2025 09:08:35 +0100
Subject: [PATCH 26/34] Revert "temporarily remove additional repos"

This reverts commit c62be58028a138f7b58b734fa58675f745707f11.
---
 DESCRIPTION | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/DESCRIPTION b/DESCRIPTION
index 739dedf18..396d0c3a1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -146,6 +146,8 @@ LinkingTo:
     RcppParallel (>= 5.0.1),
     rstan (>= 2.26.0),
     StanHeaders (>= 2.26.0)
+Additional_repositories:
+    https://stan-dev.r-universe.dev
 Biarch: true
 Config/testthat/edition: 3
 Encoding: UTF-8

From d70e4e5c3354dd1747a2261c496368f473e54237 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Tue, 21 Jan 2025 20:59:41 +0100
Subject: [PATCH 27/34] touchstone: don't upgrade

---
 .github/workflows/touchstone-receive.yaml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/touchstone-receive.yaml b/.github/workflows/touchstone-receive.yaml
index 60b8b6634..40a0a42d0 100644
--- a/.github/workflows/touchstone-receive.yaml
+++ b/.github/workflows/touchstone-receive.yaml
@@ -52,7 +52,8 @@ jobs:
         with:
           cache-version: 1
           extra-packages: stan-dev/cmdstanr
-          touchstone_ref: '@289cfc7'
+          touchstone_package: 'github::sbfnk/touchstone'
+          touchstone_ref: '@dont-upgrade'
           benchmarking_repo: ${{ matrix.config.benchmarking_repo }}
           benchmarking_ref: ${{ matrix.config.benchmarking_ref }}
           benchmarking_path: ${{ matrix.config.benchmarking_path }}

From 7f48f9ee8dccc5986e831673c6349e6557aaf498 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 22 Jan 2025 09:33:41 +0100
Subject: [PATCH 28/34] stabilise initial guess

---
 inst/stan/functions/rt.stan | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 516992b65..6d03feae8 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -156,7 +156,7 @@ real R_to_r_manual(real R, vector gt_rev_pmf, real abs_tol) {
   int gt_len = num_elements(gt_rev_pmf);
   vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
   real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));
-  real r = (R - 1) / (R * mean_gt);
+  real r = (R - 1) / (R * mean_gt + 1);
   real step = abs_tol + 1;
   while (abs(step) > abs_tol) {
     step = R_to_r_newton_step(R, r, gt_pmf);

From 684f139af1d82d8c9427b586a4f6123affa3c3a1 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 22 Jan 2025 09:34:10 +0100
Subject: [PATCH 29/34] Revert "touchstone: don't upgrade"

This reverts commit d70e4e5c3354dd1747a2261c496368f473e54237.
---
 .github/workflows/touchstone-receive.yaml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/.github/workflows/touchstone-receive.yaml b/.github/workflows/touchstone-receive.yaml
index 40a0a42d0..60b8b6634 100644
--- a/.github/workflows/touchstone-receive.yaml
+++ b/.github/workflows/touchstone-receive.yaml
@@ -52,8 +52,7 @@ jobs:
         with:
           cache-version: 1
           extra-packages: stan-dev/cmdstanr
-          touchstone_package: 'github::sbfnk/touchstone'
-          touchstone_ref: '@dont-upgrade'
+          touchstone_ref: '@289cfc7'
           benchmarking_repo: ${{ matrix.config.benchmarking_repo }}
           benchmarking_ref: ${{ matrix.config.benchmarking_ref }}
           benchmarking_path: ${{ matrix.config.benchmarking_path }}

From 75628ff0ac5c6690288ca21bda79ce17d2fc6123 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 22 Jan 2025 10:41:30 +0100
Subject: [PATCH 30/34] update sim snapshot

---
 tests/testthat/_snaps/simulate-infections.md | 28 ++++++++++----------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 8d44802eb..8bf3150f6 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -36,20 +36,20 @@
 
               variable       date    value
                 <char>     <Date>    <num>
-     1:     infections 2023-01-01 160.8333
-     2:     infections 2023-01-02 169.6942
-     3:     infections 2023-01-03 178.7393
-     4:     infections 2023-01-04 188.1775
-     5:     infections 2023-01-05 198.0770
-     6:     infections 2023-01-06 208.1120
-     7:     infections 2023-01-07 218.6883
-     8:     infections 2023-01-08 153.2081
-     9:     infections 2023-01-09 148.6709
-    10:     infections 2023-01-10 142.1679
-    11:     infections 2023-01-11 135.3769
-    12:     infections 2023-01-12 128.6677
-    13:     infections 2023-01-13 122.1667
-    14:     infections 2023-01-14 115.9185
+     1:     infections 2023-01-01 160.8325
+     2:     infections 2023-01-02 169.6934
+     3:     infections 2023-01-03 178.7385
+     4:     infections 2023-01-04 188.1766
+     5:     infections 2023-01-05 198.0760
+     6:     infections 2023-01-06 208.1110
+     7:     infections 2023-01-07 218.6873
+     8:     infections 2023-01-08 153.2073
+     9:     infections 2023-01-09 148.6702
+    10:     infections 2023-01-10 142.1672
+    11:     infections 2023-01-11 135.3762
+    12:     infections 2023-01-12 128.6671
+    13:     infections 2023-01-13 122.1661
+    14:     infections 2023-01-14 115.9180
     15: reported_cases 2023-01-01 117.0000
     16: reported_cases 2023-01-02 107.0000
     17: reported_cases 2023-01-03 119.0000

From 0449171a47d048b56ec79250e8d39909027d0b2c Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 22 Jan 2025 12:40:20 +0100
Subject: [PATCH 31/34] try max instead of correction

---
 inst/stan/functions/rt.stan                  |  2 +-
 tests/testthat/_snaps/simulate-infections.md | 28 ++++++++++----------
 2 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 6d03feae8..cf4dfd83d 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -156,7 +156,7 @@ real R_to_r_manual(real R, vector gt_rev_pmf, real abs_tol) {
   int gt_len = num_elements(gt_rev_pmf);
   vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
   real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));
-  real r = (R - 1) / (R * mean_gt + 1);
+  real r = fmax((R - 1) / (R * mean_gt), -1);
   real step = abs_tol + 1;
   while (abs(step) > abs_tol) {
     step = R_to_r_newton_step(R, r, gt_pmf);
diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md
index 8bf3150f6..8d44802eb 100644
--- a/tests/testthat/_snaps/simulate-infections.md
+++ b/tests/testthat/_snaps/simulate-infections.md
@@ -36,20 +36,20 @@
 
               variable       date    value
                 <char>     <Date>    <num>
-     1:     infections 2023-01-01 160.8325
-     2:     infections 2023-01-02 169.6934
-     3:     infections 2023-01-03 178.7385
-     4:     infections 2023-01-04 188.1766
-     5:     infections 2023-01-05 198.0760
-     6:     infections 2023-01-06 208.1110
-     7:     infections 2023-01-07 218.6873
-     8:     infections 2023-01-08 153.2073
-     9:     infections 2023-01-09 148.6702
-    10:     infections 2023-01-10 142.1672
-    11:     infections 2023-01-11 135.3762
-    12:     infections 2023-01-12 128.6671
-    13:     infections 2023-01-13 122.1661
-    14:     infections 2023-01-14 115.9180
+     1:     infections 2023-01-01 160.8333
+     2:     infections 2023-01-02 169.6942
+     3:     infections 2023-01-03 178.7393
+     4:     infections 2023-01-04 188.1775
+     5:     infections 2023-01-05 198.0770
+     6:     infections 2023-01-06 208.1120
+     7:     infections 2023-01-07 218.6883
+     8:     infections 2023-01-08 153.2081
+     9:     infections 2023-01-09 148.6709
+    10:     infections 2023-01-10 142.1679
+    11:     infections 2023-01-11 135.3769
+    12:     infections 2023-01-12 128.6677
+    13:     infections 2023-01-13 122.1667
+    14:     infections 2023-01-14 115.9185
     15: reported_cases 2023-01-01 117.0000
     16: reported_cases 2023-01-02 107.0000
     17: reported_cases 2023-01-03 119.0000

From dba0c2ef34412f88d5ea527818dbe41bb7e9ec61 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 22 Jan 2025 14:18:35 +0100
Subject: [PATCH 32/34] remove/rename

---
 inst/stan/functions/infections.stan |  2 +-
 inst/stan/functions/rt.stan         | 62 ++---------------------------
 2 files changed, 4 insertions(+), 60 deletions(-)

diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan
index f9de031eb..38a0f4b17 100644
--- a/inst/stan/functions/infections.stan
+++ b/inst/stan/functions/infections.stan
@@ -30,7 +30,7 @@ vector generate_infections(vector R, int uot, vector gt_rev_pmf,
   vector[t] infections = rep_vector(0, t);
   vector[ot] cum_infections;
   vector[ot] infectiousness;
-  real growth = R_to_r_manual(R[1], gt_rev_pmf, 1e-3);
+  real growth = R_to_r(R[1], gt_rev_pmf, 1e-3);
   // Initialise infections using daily growth
   if (initial_as_scale) {
     infections[1] = exp(initial_infections[1] - growth * uot);
diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index cf4dfd83d..53ab53283 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -64,27 +64,10 @@ void rt_lp(array[] real initial_infections_scale, vector bp_effects,
   initial_infections_scale ~ normal(initial_infections_guess, 2);
 }
 
-/**
- * Calculate the negative moment generating function (nMGF) of a probability
- * distribution with given probability mass function.
- *
- * Function code is based on Julia code from
- * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
- * under Apache license 2.0.
- *
- * @param r function argument of the nMGF
- * @param pmf probability mass function as vector (first index: 0)
- */
-  real neg_MGF(real r, vector pmf) {
-    int len = num_elements(pmf);
-    vector[len] exp_r = exp(-r * linspaced_vector(len, 0, len - 1));
-    return(dot_product(pmf, exp_r));
-  }
-
 /**
  * Helper function for calculating r from R using Newton's method
  *
- * Function code is based on Julia code from
+ * Code is based on Julia code from
  * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
  * under Apache license 2.0.
  *
@@ -101,50 +84,11 @@ real R_to_r_newton_step(real R, real r, vector pmf) {
   return(ret);
 }
 
-/**
- * Helper function used by the [R_to_r()] function to estimate r from R
- *
- * Function code is based on Julia code from
- * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
- * under Apache license 2.0.
- *
- * @param r A vector of length 1, the growth rate
- * @param gt_pmf probability mass function of the generation time
- * @param R reproduction number
- */
-vector eq(vector r, vector gt_pmf, real R) {
-  return([ neg_MGF(r[1], gt_pmf) - 1 / R ]');
-}
-
-/**
- * Estimate the growth rate r from reproduction number R. Used in the model to
- * estimate the initial growth rate using stan's algebraic solver
- *
- * Function code is based on Julia code from
- * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
- * under Apache license 2.0.
- *
- * @param R reproduction number
- * @param gt_rev_pmf reverse probability mass function of the generation time
- */
-real R_to_r(real R, vector gt_rev_pmf, int newton_steps) {
-  int gt_len = num_elements(gt_rev_pmf);
-  vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
-  real mean_gt = dot_product(
-    gt_pmf,
-    linspaced_vector(gt_len, 0, gt_len - 1)
-  );
-  vector[1] r_approx = [ (R - 1) / (R * mean_gt) ]';
-  vector[1] r = solve_powell(eq, r_approx, gt_pmf, R);
-
-  return(r[1]);
-}
-
 /**
  * Estimate the growth rate r from reproduction number R. Used in the model to
  * estimate the initial growth rate using Newton's method.
  *
- * Function code is based on Julia code from
+ * Code is based on Julia code from
  * https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl
  * under Apache license 2.0.
  *
@@ -152,7 +96,7 @@ real R_to_r(real R, vector gt_rev_pmf, int newton_steps) {
  * @param gt_rev_pmf reverse probability mass function of the generation time
  * @param abs_tol absolute tolerance of the solver
  */
-real R_to_r_manual(real R, vector gt_rev_pmf, real abs_tol) {
+real R_to_r(real R, vector gt_rev_pmf, real abs_tol) {
   int gt_len = num_elements(gt_rev_pmf);
   vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
   real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));

From 6a360589232286967a9651e4f918527d1c58926d Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Thu, 23 Jan 2025 14:36:09 +0000
Subject: [PATCH 33/34] version 2: testing for speed

---
 inst/stan/functions/rt.stan | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 53ab53283..2326e5cd6 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -100,7 +100,7 @@ real R_to_r(real R, vector gt_rev_pmf, real abs_tol) {
   int gt_len = num_elements(gt_rev_pmf);
   vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
   real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));
-  real r = fmax((R - 1) / (R * mean_gt), -1);
+  real r = (R - 1) / (R * mean_gt + 1);
   real step = abs_tol + 1;
   while (abs(step) > abs_tol) {
     step = R_to_r_newton_step(R, r, gt_pmf);

From 35740e58c5deb8965d5010777db7e9e7514bc843 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Thu, 23 Jan 2025 15:52:07 +0000
Subject: [PATCH 34/34] Revert "version 2: testing for speed"

This reverts commit 6a360589232286967a9651e4f918527d1c58926d.
---
 inst/stan/functions/rt.stan | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
index 2326e5cd6..53ab53283 100644
--- a/inst/stan/functions/rt.stan
+++ b/inst/stan/functions/rt.stan
@@ -100,7 +100,7 @@ real R_to_r(real R, vector gt_rev_pmf, real abs_tol) {
   int gt_len = num_elements(gt_rev_pmf);
   vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
   real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));
-  real r = (R - 1) / (R * mean_gt + 1);
+  real r = fmax((R - 1) / (R * mean_gt), -1);
   real step = abs_tol + 1;
   while (abs(step) > abs_tol) {
     step = R_to_r_newton_step(R, r, gt_pmf);