Merge branch 'main' into accumulate-design

NEWS.md

@@ -3,6 +3,7 @@
 ## Model changes
 
 - A bug in the spectral densities of Matern kernels of order other than 3/2 has been fixed and the vignette updated accordingly. By @sbfnk in #835 and reviewed by @seabbs.
+- Changed the prior on the (square root of the) magnitude alpha of the Gaussian Process back to HalfNormal(0, 0.05) based on user feedback of unexpected results. By @sbfnk in #840 and reviewed by @jamesmbaazam.
 
 ## Package changes
 
@@ -25,7 +26,7 @@ A release that introduces model improvements to the Gaussian Process models, alo
 - When defining probability distributions these can now be truncated using the `tolerance` argument
 - Ornstein-Uhlenbeck and 5 / 2 Matérn kernels have been added. By @sbfnk in #741 and reviewed by @seabbs.
 - Gaussian processes have been vectorised, leading to some speed gains 🚀 , and the `gp_opts()` function has gained three more options, "periodic", "ou",  and "se", to specify periodic and linear kernels respectively. By @seabbs in #742 and reviewed by @jamesmbaazam.
-- Prior predictive checks have been used to update the following priors: the prior on the magnitude of the Gaussian process (from HalfNormal(0, 1) to HalfNormal(0, 0.1)), and the prior on the overdispersion (from 1 / HalfNormal(0, 1)^2 to 1 / HalfNormal(0, 0.25)). In the user-facing API, this is a change in default values of the `sd` of `phi` in `obs_opts()` from 1 to 0.25. By @seabbs in #742 and reviewed by @jamesmbaazam.
+- Prior predictive checks have been used to update the following priors: the prior on the magnitude of the Gaussian process (from HalfNormal(0, 0.05)^2 to HalfNormal(0, 0.01)^2), and the prior on the overdispersion (from 1 / HalfNormal(0, 1)^2 to 1 / HalfNormal(0, 0.25)). In the user-facing API, this is a change in default values of the `sd` of `phi` in `obs_opts()` from 1 to 0.25. By @seabbs in #742 and reviewed by @jamesmbaazam.
 - The default stan control options have been updated from `list(adapt_delta = 0.95, max_treedepth = 15)` to `list(adapt_delta = 0.9, max_treedepth = 12)` due to improved performance and to reduce the runtime of the default parameterisations. By @seabbs in #742 and reviewed by @jamesmbaazam.
 - Initialisation has been simplified by sampling directly from the priors, where possible, rather than from a constrained space. By @seabbs in #742 and reviewed by @jamesmbaazam.
 - Unnecessary normalisation of delay priors has been removed. By @seabbs in #742 and reviewed by @jamesmbaazam.

R/opts.R

@@ -458,7 +458,7 @@ backcalc_opts <- function(prior = c("reports", "none", "infections"),
 #' deviation of the Gaussian process (logged Rt in case of the renewal model,
 #' logged infections in case of the nonmechanistic model).
 #'
-#' @param alpha_sd Numeric, defaults to 0.1. The standard deviation of the
+#' @param alpha_sd Numeric, defaults to 0.05. The standard deviation of the
 #' magnitude parameter of the Gaussian process kernel. Can be tuned to adjust
 #' how far alpha is allowed to deviate form its prior mean (`alpha_mean`).
 #'
@@ -509,7 +509,7 @@ gp_opts <- function(basis_prop = 0.2,
                     ls_min = 0,
                     ls_max = 60,
                     alpha_mean = 0,
-                    alpha_sd = 0.01,
+                    alpha_sd = 0.05,
                     kernel = c("matern", "se", "ou", "periodic"),
                     matern_order = 3 / 2,
                     matern_type,

inst/stan/functions/gaussian_process.stan

@@ -33,7 +33,7 @@ vector diagSPD_EQ(real alpha, real rho, real L, int M) {
 vector diagSPD_Matern12(real alpha, real rho, real L, int M) {
   vector[M] indices = linspaced_vector(M, 1, M);
   real factor = 2;
-  vector[M] denom = rho * ((1 / rho)^2 + (pi() / 2 / L) * indices);
+  vector[M] denom = rho * ((1 / rho)^2 + pow(pi() / 2 / L * indices, 2));
   return alpha * sqrt(factor * inv(denom));
 }
 
@@ -65,7 +65,7 @@ vector diagSPD_Matern32(real alpha, real rho, real L, int M) {
 vector diagSPD_Matern52(real alpha, real rho, real L, int M) {
   vector[M] indices = linspaced_vector(M, 1, M);
   real factor = 3 * pow(sqrt(5) / rho, 5);
-  vector[M] denom = 2 * (sqrt(5) / rho)^2 + pow((pi() / 2 / L) * indices, 3);
+  vector[M] denom = 2 * pow((sqrt(5) / rho)^2 + pow((pi() / 2 / L) * indices, 2), 3);
   return alpha * sqrt(factor * inv(denom));
 }
 

tests/testthat/test-create_gp_data.R

@@ -11,7 +11,7 @@ test_that("create_gp_data returns correct default values when GP is disabled", {
   expect_equal(gp_data$ls_sdlog, convert_to_logsd(21, 7))
   expect_equal(gp_data$ls_min, 0)
   expect_equal(gp_data$ls_max, 3.54, tolerance = 0.01)
-  expect_equal(gp_data$alpha_sd, 0.01)
+  expect_equal(gp_data$alpha_sd, 0.05)
   expect_equal(gp_data$gp_type, 2)  # Default to Matern
   expect_equal(gp_data$nu, 3 / 2)
   expect_equal(gp_data$w0, 1.0)

tests/testthat/test-gp_opts.R

@@ -6,7 +6,7 @@ test_that("gp_opts returns correct default values", {
   expect_equal(gp$ls_sd, 7)
   expect_equal(gp$ls_min, 0)
   expect_equal(gp$ls_max, 60)
-  expect_equal(gp$alpha_sd, 0.01)
+  expect_equal(gp$alpha_sd, 0.05)
   expect_equal(gp$kernel, "matern")
   expect_equal(gp$matern_order, 3 / 2)
   expect_equal(gp$w0, 1.0)

tests/testthat/test-stan-guassian-process.R

@@ -36,7 +36,7 @@ test_that("diagSPD_Matern functions return correct dimensions and values", {
   # Check specific values for known inputs
   indices <- linspaced_vector(M, 1, M)
   factor12 <- 2
-  denom12 <- rho * ((1 / rho)^2 + (pi / 2 / L) * indices)
+  denom12 <- rho * ((1 / rho)^2 + (pi / 2 / L * indices)^2)
   expected_result12 <- alpha * sqrt(factor12 / denom12)
   expect_equal(result12, expected_result12, tolerance = 1e-8)
 
@@ -58,7 +58,7 @@ test_that("diagSPD_Matern functions return correct dimensions and values", {
   # Check specific values for known inputs
   indices <- linspaced_vector(M, 1, M)
   factor52 <- 3 * (sqrt(5) / rho)^5
-  denom52 <- 2 * (sqrt(5) / rho)^2 + ((pi / 2 / L) * indices)^3
+  denom52 <- 2 * ((sqrt(5) / rho)^2 + (pi / 2 / L * indices)^2)^3
   expected_result52 <- alpha * sqrt(factor52 / denom52)
   expect_equal(result52, expected_result52, tolerance = 1e-8)
 })

vignettes/gaussian_process_implementation_details.Rmd

@@ -104,28 +104,28 @@ S_k(x) = \alpha^2 \frac{2 \sqrt{\pi} \Gamma(\nu + 1/2) (2\nu)^\nu}{\Gamma(\nu) \
 For $\nu = 3 / 2$ (the default in `EpiNow2`) this simplifies to
 
 \begin{equation}
-    S_k(x) =
+    S_k(\omega) =
     \alpha^2 \frac{4 \left(\sqrt{3} / \rho \right)^3}{\left(\left(\sqrt{3} / \rho\right)^2 + \omega^2\right)^2} = 
     \left(\frac{2 \alpha \left(\sqrt{3} / \rho \right)^{\frac{3}{2}}}{\left(\sqrt{3} / \rho\right)^2 + \omega^2} \right)^2
 \end{equation}
 
 For $\nu = 1 / 2$ it is
 
 \begin{equation}
-S_k(x) = \alpha^2 \frac{2}{\rho \left(1 / \rho^2 + \omega^2\right)}
+S_k(\omega) = \alpha^2 \frac{2}{\rho \left(1 / \rho^2 + \omega^2\right)}
 \end{equation}
 
 and for $\nu = 5 / 2$ it is
 
 \begin{equation}
-S_k(x) =
+S_k(\omega) =
     \alpha^2 \frac{3 \left(\sqrt{5} / \rho \right)^5}{2 \left(\left(\sqrt{5} / \rho\right)^2 + \omega^2\right)^3}
 \end{equation}
 
 In the case of a squared exponential the spectral kernel density is given by
 
 \begin{equation}
-S_k(x) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 w^2 \right)
+S_k(\omega) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 \omega^2 \right)
 \end{equation}
 
 The functions $\phi_{j}(x)$ are the eigenfunctions of the Laplace operator,